polyquantique / symplecticdecompositions.jl Goto Github PK

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 I compute the decomposition to get the elements of the diagonal matrix S and the unitary W. However, it seems to fail as τ ≠ W * Diagonal(S) * transpose(W).

τ = [0.925 0.0 0.0 0.0; 0.0 -0.02399982992272 -0.00489937871047 -0.00500042517513; 0.0 -0.00489937871047 -0.00100017007728 0.02449481063548; 0.0 -0.00500042517513 0.02449481063548 0.0]
issymmetric(τ)                                       # true
S, W = takagi(τ);
τ ≈ W * Diagonal(S) * transpose(W)                   # false 
norm(τ - W * Diagonal(S) * transpose(W))             # 0.14488271314030668

For most other examples I ran it works nicely and I get a norm error of the order of $10^{-17}$, but this example just dramatically fails. Here is one very close example where it works:

τ = [0.9325 0.0 0.0 -0.0; 0.0 -0.0135 0.018 0.0; 0.0 0.018 0.0135 0.0; -0.0 0.0 0.0 -0.0225]
issymmetric(τ)                                       # true
S, W = takagi(τ);
τ ≈ W * Diagonal(S) * transpose(W)                   # true
norm(τ - W * Diagonal(S) * transpose(W))             # 1.758776770179577e-17

I am probably doing something wrong, but since I can't find what I thought to share in case it is some edge case. Do you know if I am doing something wrong, or there could be some numerical instability somewhere?