Factorization failing?

[SimonCichy]

Hi!
I need a Takagi factorization for a project I am working on, so I tried to use your package. However, in some cases it fails and I don't know why. Here is an example:
I start with a symmetric matrix τ and then as in the Readme I compute the decomposition to get the diagonal matrix D and the unitary W. However, it seems to fail as τ ≠ transpose(W) * D * W. I also don't get any convergence error

τ = ComplexF64[0.925 + 0.0im 0.0 + 0.0im 0.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im -0.02399982992272 + 0.0im -0.00489937871047 + 0.0im -0.00500042517513 + 0.0im; 0.0 + 0.0im -0.00489937871047 + 0.0im -0.00100017007728 + 0.0im 0.02449481063548 + 0.0im; 0.0 + 0.0im -0.00500042517513 + 0.0im 0.02449481063548 + 0.0im 0.0 + 0.0im]
issymmetric(τ)                             # true
D, W = takagi_factor(τ; sort=-1)
τ ≈ transpose(W) * D * W                   # false 

I am probably doing something wrong, but since I can't find what I thought to share in case it is some edge case.

Also, as a side note, I found this paper which seems to give a closed-form formula for the Takagi decomposition based on the SVD (Section IV, in particular Box 1). It may be interesting to complement the other issue, but in my case there are some instances where it fails dramatically and I don't know why.

[JLTastet]

Hi Simon,

This package is a fairly dumb Julia reimplementation of this routine (possibly an earlier version) that is used in a popular high-energy physics software. I am not sure if it has been broadly validated, or if there are any guarantees on its accuracy (probably not, given the issue you are encountering).

It would be amazing if you could check if the issue also appears in the original (C or Fortran) version. If not, then it’s a bug on my side, and I can try to fix it.

If the original version also features the same issue, then we’ll have to find a better algorithm.

I am open to contributions, and even to adding new maintainers, since I don’t use this package that much nowadays :)

[SimonCichy]

Hi Jean-Loup,

Thanks a lot for the fast reply, even though you don't actively work on this currently. I have now written an implementation based on the singular value decomposition from this paper and their corresponding julia implementation SymplecticDecompositions.jl, combined with an approach based on the eigenvalue decomposition for real matrices, taken from Strawbery Fields, a python package for simulation and programming quantum devices.

If you are interested in my implementation (which is not iterative) I could share it, but even if it works for the cases I tried it is not thouroughly tested, and none of it are my ideas.

[JLTastet]

Hi Simon,

I would be glad if you could contribute your implementation!

If you are confident that it covers more corner cases than the current one, then we could make it the default. Otherwise, we can have both implementations coexist side by side. And in either case, I should probably emphasise in the README that this package is still experimental and that user discretion is advised.

So feel free to open a pull request!