τ = [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?