search

JuliaLinearAlgebra / SkewLinearAlgebra.jl

Julia Package for skew-symmetric matrices
MIT License
17 stars 4 forks source link

eigen fails to diagonalize a matrix #130

Closed cvsvensson closed 9 months ago

cvsvensson commented 9 months ago

I encountered a matrix that eigen fails to diagonalize. Both the values and vectors are wrong. However, eigvals returns the correct eigenvalues.

using SkewLinearAlgebra
#sp is the sparse matrix in the code block below
A = SkewHermitian(Matrix(sp))
eigvals(A) #These are the correct, mostly non-zero, eigenvalues. Agrees with eigvals(Matrix(A)).
vals, vecs = eigen(A)#vals are all zero, which is incorrect.
vecs'*A*vecs #Is not diagonal. It's tridiagonal.
svd(A).S #Also all zeros.
using SparseArrays
sp = sparse([26, 50, 51, 52, 27, 51, 52, 53, 28, 52, 53, 54, 29, 53, 54, 55, 30, 54, 55, 56, 31, 55, 56, 32, 33, 57, 58, 32, 33, 34, 57, 58, 59, 33, 34, 35, 58, 59, 60, 34, 35, 36, 59, 60, 61, 35, 36, 37, 60, 61, 62, 36, 37, 38, 61, 62, 63, 37, 38, 39, 62, 63, 64, 38, 39, 40, 63, 64, 65, 39, 40, 41, 64, 65, 66, 40, 41, 42, 65, 66, 41, 42, 43, 66, 42, 43, 44, 43, 44, 45, 44, 45, 46, 45, 46, 47, 46, 47, 48, 47, 48, 49, 48, 49, 50, 49, 50, 51, 1, 50, 51, 52, 2, 51, 52, 53, 3, 52, 53, 54, 4, 53, 54, 55, 5, 54, 55, 56, 6, 55, 56, 7, 8, 7, 8, 9, 58, 8, 9, 10, 59, 9, 10, 11, 60, 10, 11, 12, 61, 11, 12, 13, 62, 12, 13, 14, 63, 13, 14, 15, 64, 14, 15, 16, 65, 15, 16, 17, 66, 16, 17, 18, 17, 18, 19, 18, 19, 20, 19, 20, 21, 20, 21, 22, 21, 22, 23, 22, 23, 24, 23, 24, 25, 1, 24, 25, 26, 1, 2, 25, 26, 27, 1, 2, 3, 26, 27, 28, 2, 3, 4, 27, 28, 29, 3, 4, 5, 28, 29, 30, 4, 5, 6, 29, 30, 31, 5, 6, 30, 31, 7, 8, 7, 8, 9, 33, 8, 9, 10, 34, 9, 10, 11, 35, 10, 11, 12, 36, 11, 12, 13, 37, 12, 13, 14, 38, 13, 14, 15, 39, 14, 15, 16, 40, 15, 16, 17, 41], [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 33, 33, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 37, 38, 38, 38, 38, 39, 39, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60, 60, 60, 61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 64, 64, 64, 64, 65, 65, 65, 65, 66, 66, 66, 66], [-1.176e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690911, -4.9e-14, -0.0980580675690911, 0.987744983947536, -0.0980580675690927, 1.96e-14, -0.0980580675690927, 0.987744983947536, -0.098058067569091, -1.96e-13, -0.098058067569091, 0.987744983947536, -0.0980580675690928, -1.274e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690928, -5.88e-14, -0.0980580675690928, 0.987744983947536, -10.0, -0.4902903378454601, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.0980580675690921, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690921, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.49029033784546, -0.0980580675690919, -100.97143868952186, -0.0980580675690923, 0.49029033784546, -10.0, -0.4902903378454601, -0.0980580675690923, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454602, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 1.176e-13, 0.4902903378454599, -10.0, -0.4902903378454602, 4.9e-14, 0.4902903378454602, -10.0, -0.4902903378454599, -1.96e-14, 0.4902903378454599, -10.0, -0.4902903378454603, 1.96e-13, 0.4902903378454603, -10.0, -0.4902903378454599, 1.274e-13, 0.4902903378454599, -10.0, -0.4902903378454599, 5.88e-14, 0.4902903378454599, -10.0, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454601, 2.4e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 4.9e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 7.3e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 9.8e-15, 0.4902903378454601, 10.0, -0.49029033784546, 1.22e-14, 0.49029033784546, 10.0, -0.4902903378454601, 1.47e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.71e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.96e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 2.2e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454603, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454602, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, -0.987744983947536, 0.0980580675690911, 0.4902903378454599, 10.0, -0.4902903378454602, 0.0980580675690911, -0.987744983947536, 0.0980580675690927, 0.4902903378454602, 10.0, -0.4902903378454599, 0.0980580675690927, -0.987744983947536, 0.098058067569091, 0.4902903378454599, 10.0, -0.4902903378454603, 0.098058067569091, -0.987744983947536, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.0980580675690928, 0.4902903378454599, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.4902903378454599, 10.0, 100.97143868952186, 0.098058067569092, 0.098058067569092, 100.97143868952186, 0.098058067569092, -2.4e-15, 0.098058067569092, 100.97143868952186, 0.0980580675690921, -4.9e-15, 0.0980580675690921, 100.97143868952186, 0.0980580675690919, -7.3e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -9.8e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690923, -1.22e-14, 0.0980580675690923, 100.97143868952186, 0.0980580675690919, -1.47e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.71e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.96e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -2.2e-14], 66, 66)
cvsvensson commented 9 months ago

Found a smaller 10x10 matrix with the same problem.

sp = sparse([2, 7, 1, 3, 6, 8, 2, 4, 7, 9, 3, 5, 8, 10, 4, 9, 2, 7, 1, 3, 6, 8, 2, 4, 7, 9, 3, 5, 8, 10, 4, 9], [1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10], [-0.8414709848, 1.5403023059, 0.8414709848, -0.8414709848, 0.4596976941, 1.5403023059, 0.8414709848, -0.8414709848, 0.4596976941, 1.5403023059, 0.8414709848, -0.8414709848, 0.4596976941, 1.5403023059, 0.8414709848, 0.4596976941, -0.4596976941, 0.8414709848, -1.5403023059, -0.4596976941, -0.8414709848, 0.8414709848, -1.5403023059, -0.4596976941, -0.8414709848, 0.8414709848, -1.5403023059, -0.4596976941, -0.8414709848, 0.8414709848, -1.5403023059, -0.8414709848], 10, 10)
smataigne commented 9 months ago

I encountered a matrix that eigen fails to diagonalize. Both the values and vectors are wrong. However, eigvals returns the correct eigenvalues.

using SkewLinearAlgebra
#sp is the sparse matrix in the code block below
A = SkewHermitian(Matrix(sp))
eigvals(A) #These are the correct, mostly non-zero, eigenvalues. Agrees with eigvals(Matrix(A)).
vals, vecs = eigen(A)#vals are all zero, which is incorrect.
vecs'*A*vecs #Is not diagonal. It's tridiagonal.
svd(A).S #Also all zeros.
using SparseArrays
sp = sparse([26, 50, 51, 52, 27, 51, 52, 53, 28, 52, 53, 54, 29, 53, 54, 55, 30, 54, 55, 56, 31, 55, 56, 32, 33, 57, 58, 32, 33, 34, 57, 58, 59, 33, 34, 35, 58, 59, 60, 34, 35, 36, 59, 60, 61, 35, 36, 37, 60, 61, 62, 36, 37, 38, 61, 62, 63, 37, 38, 39, 62, 63, 64, 38, 39, 40, 63, 64, 65, 39, 40, 41, 64, 65, 66, 40, 41, 42, 65, 66, 41, 42, 43, 66, 42, 43, 44, 43, 44, 45, 44, 45, 46, 45, 46, 47, 46, 47, 48, 47, 48, 49, 48, 49, 50, 49, 50, 51, 1, 50, 51, 52, 2, 51, 52, 53, 3, 52, 53, 54, 4, 53, 54, 55, 5, 54, 55, 56, 6, 55, 56, 7, 8, 7, 8, 9, 58, 8, 9, 10, 59, 9, 10, 11, 60, 10, 11, 12, 61, 11, 12, 13, 62, 12, 13, 14, 63, 13, 14, 15, 64, 14, 15, 16, 65, 15, 16, 17, 66, 16, 17, 18, 17, 18, 19, 18, 19, 20, 19, 20, 21, 20, 21, 22, 21, 22, 23, 22, 23, 24, 23, 24, 25, 1, 24, 25, 26, 1, 2, 25, 26, 27, 1, 2, 3, 26, 27, 28, 2, 3, 4, 27, 28, 29, 3, 4, 5, 28, 29, 30, 4, 5, 6, 29, 30, 31, 5, 6, 30, 31, 7, 8, 7, 8, 9, 33, 8, 9, 10, 34, 9, 10, 11, 35, 10, 11, 12, 36, 11, 12, 13, 37, 12, 13, 14, 38, 13, 14, 15, 39, 14, 15, 16, 40, 15, 16, 17, 41], [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 33, 33, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 37, 38, 38, 38, 38, 39, 39, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60, 60, 60, 61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 64, 64, 64, 64, 65, 65, 65, 65, 66, 66, 66, 66], [-1.176e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690911, -4.9e-14, -0.0980580675690911, 0.987744983947536, -0.0980580675690927, 1.96e-14, -0.0980580675690927, 0.987744983947536, -0.098058067569091, -1.96e-13, -0.098058067569091, 0.987744983947536, -0.0980580675690928, -1.274e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690928, -5.88e-14, -0.0980580675690928, 0.987744983947536, -10.0, -0.4902903378454601, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.0980580675690921, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690921, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.49029033784546, -0.0980580675690919, -100.97143868952186, -0.0980580675690923, 0.49029033784546, -10.0, -0.4902903378454601, -0.0980580675690923, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454602, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 1.176e-13, 0.4902903378454599, -10.0, -0.4902903378454602, 4.9e-14, 0.4902903378454602, -10.0, -0.4902903378454599, -1.96e-14, 0.4902903378454599, -10.0, -0.4902903378454603, 1.96e-13, 0.4902903378454603, -10.0, -0.4902903378454599, 1.274e-13, 0.4902903378454599, -10.0, -0.4902903378454599, 5.88e-14, 0.4902903378454599, -10.0, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454601, 2.4e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 4.9e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 7.3e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 9.8e-15, 0.4902903378454601, 10.0, -0.49029033784546, 1.22e-14, 0.49029033784546, 10.0, -0.4902903378454601, 1.47e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.71e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.96e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 2.2e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454603, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454602, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, -0.987744983947536, 0.0980580675690911, 0.4902903378454599, 10.0, -0.4902903378454602, 0.0980580675690911, -0.987744983947536, 0.0980580675690927, 0.4902903378454602, 10.0, -0.4902903378454599, 0.0980580675690927, -0.987744983947536, 0.098058067569091, 0.4902903378454599, 10.0, -0.4902903378454603, 0.098058067569091, -0.987744983947536, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.0980580675690928, 0.4902903378454599, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.4902903378454599, 10.0, 100.97143868952186, 0.098058067569092, 0.098058067569092, 100.97143868952186, 0.098058067569092, -2.4e-15, 0.098058067569092, 100.97143868952186, 0.0980580675690921, -4.9e-15, 0.0980580675690921, 100.97143868952186, 0.0980580675690919, -7.3e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -9.8e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690923, -1.22e-14, 0.0980580675690923, 100.97143868952186, 0.0980580675690919, -1.47e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.71e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.96e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -2.2e-14], 66, 66)

Thank you very much for reporting the issue! eigvals and eigen should output the same result since they use the same QR algorithm. Your example shows that there must be a bug in the routine skewtrieigen_divided! (which underlies eigen) that is not present in skewtrieigvals (which underlies eigvals). I will try to find where the difference comes from.

smataigne commented 9 months ago

I encountered a matrix that eigen fails to diagonalize. Both the values and vectors are wrong. However, eigvals returns the correct eigenvalues.

using SkewLinearAlgebra
#sp is the sparse matrix in the code block below
A = SkewHermitian(Matrix(sp))
eigvals(A) #These are the correct, mostly non-zero, eigenvalues. Agrees with eigvals(Matrix(A)).
vals, vecs = eigen(A)#vals are all zero, which is incorrect.
vecs'*A*vecs #Is not diagonal. It's tridiagonal.
svd(A).S #Also all zeros.
using SparseArrays
sp = sparse([26, 50, 51, 52, 27, 51, 52, 53, 28, 52, 53, 54, 29, 53, 54, 55, 30, 54, 55, 56, 31, 55, 56, 32, 33, 57, 58, 32, 33, 34, 57, 58, 59, 33, 34, 35, 58, 59, 60, 34, 35, 36, 59, 60, 61, 35, 36, 37, 60, 61, 62, 36, 37, 38, 61, 62, 63, 37, 38, 39, 62, 63, 64, 38, 39, 40, 63, 64, 65, 39, 40, 41, 64, 65, 66, 40, 41, 42, 65, 66, 41, 42, 43, 66, 42, 43, 44, 43, 44, 45, 44, 45, 46, 45, 46, 47, 46, 47, 48, 47, 48, 49, 48, 49, 50, 49, 50, 51, 1, 50, 51, 52, 2, 51, 52, 53, 3, 52, 53, 54, 4, 53, 54, 55, 5, 54, 55, 56, 6, 55, 56, 7, 8, 7, 8, 9, 58, 8, 9, 10, 59, 9, 10, 11, 60, 10, 11, 12, 61, 11, 12, 13, 62, 12, 13, 14, 63, 13, 14, 15, 64, 14, 15, 16, 65, 15, 16, 17, 66, 16, 17, 18, 17, 18, 19, 18, 19, 20, 19, 20, 21, 20, 21, 22, 21, 22, 23, 22, 23, 24, 23, 24, 25, 1, 24, 25, 26, 1, 2, 25, 26, 27, 1, 2, 3, 26, 27, 28, 2, 3, 4, 27, 28, 29, 3, 4, 5, 28, 29, 30, 4, 5, 6, 29, 30, 31, 5, 6, 30, 31, 7, 8, 7, 8, 9, 33, 8, 9, 10, 34, 9, 10, 11, 35, 10, 11, 12, 36, 11, 12, 13, 37, 12, 13, 14, 38, 13, 14, 15, 39, 14, 15, 16, 40, 15, 16, 17, 41], [1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 11, 11, 11, 11, 11, 11, 12, 12, 12, 12, 12, 12, 13, 13, 13, 13, 13, 13, 14, 14, 14, 14, 14, 14, 15, 15, 15, 15, 15, 15, 16, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 18, 19, 19, 19, 20, 20, 20, 21, 21, 21, 22, 22, 22, 23, 23, 23, 24, 24, 24, 25, 25, 25, 26, 26, 26, 26, 27, 27, 27, 27, 28, 28, 28, 28, 29, 29, 29, 29, 30, 30, 30, 30, 31, 31, 31, 32, 32, 33, 33, 33, 33, 34, 34, 34, 34, 35, 35, 35, 35, 36, 36, 36, 36, 37, 37, 37, 37, 38, 38, 38, 38, 39, 39, 39, 39, 40, 40, 40, 40, 41, 41, 41, 41, 42, 42, 42, 43, 43, 43, 44, 44, 44, 45, 45, 45, 46, 46, 46, 47, 47, 47, 48, 48, 48, 49, 49, 49, 50, 50, 50, 50, 51, 51, 51, 51, 51, 52, 52, 52, 52, 52, 52, 53, 53, 53, 53, 53, 53, 54, 54, 54, 54, 54, 54, 55, 55, 55, 55, 55, 55, 56, 56, 56, 56, 57, 57, 58, 58, 58, 58, 59, 59, 59, 59, 60, 60, 60, 60, 61, 61, 61, 61, 62, 62, 62, 62, 63, 63, 63, 63, 64, 64, 64, 64, 65, 65, 65, 65, 66, 66, 66, 66], [-1.176e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690911, -4.9e-14, -0.0980580675690911, 0.987744983947536, -0.0980580675690927, 1.96e-14, -0.0980580675690927, 0.987744983947536, -0.098058067569091, -1.96e-13, -0.098058067569091, 0.987744983947536, -0.0980580675690928, -1.274e-13, -0.0980580675690928, 0.987744983947536, -0.0980580675690928, -5.88e-14, -0.0980580675690928, 0.987744983947536, -10.0, -0.4902903378454601, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.098058067569092, 0.4902903378454601, -10.0, -0.4902903378454601, -0.098058067569092, -100.97143868952186, -0.0980580675690921, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690921, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.49029033784546, -0.0980580675690919, -100.97143868952186, -0.0980580675690923, 0.49029033784546, -10.0, -0.4902903378454601, -0.0980580675690923, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, -100.97143868952186, 0.4902903378454601, -10.0, -0.4902903378454601, -0.0980580675690919, 0.4902903378454601, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454602, 0.4902903378454603, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454599, 0.4902903378454599, -10.0, -0.4902903378454603, 0.4902903378454603, -10.0, -0.4902903378454599, 1.176e-13, 0.4902903378454599, -10.0, -0.4902903378454602, 4.9e-14, 0.4902903378454602, -10.0, -0.4902903378454599, -1.96e-14, 0.4902903378454599, -10.0, -0.4902903378454603, 1.96e-13, 0.4902903378454603, -10.0, -0.4902903378454599, 1.274e-13, 0.4902903378454599, -10.0, -0.4902903378454599, 5.88e-14, 0.4902903378454599, -10.0, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454601, 2.4e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 4.9e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 7.3e-15, 0.4902903378454601, 10.0, -0.4902903378454601, 9.8e-15, 0.4902903378454601, 10.0, -0.49029033784546, 1.22e-14, 0.49029033784546, 10.0, -0.4902903378454601, 1.47e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.71e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 1.96e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 2.2e-14, 0.4902903378454601, 10.0, -0.4902903378454601, 0.4902903378454601, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454603, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.4902903378454602, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454599, 0.4902903378454599, 10.0, -0.4902903378454603, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, -0.987744983947536, 0.0980580675690911, 0.4902903378454599, 10.0, -0.4902903378454602, 0.0980580675690911, -0.987744983947536, 0.0980580675690927, 0.4902903378454602, 10.0, -0.4902903378454599, 0.0980580675690927, -0.987744983947536, 0.098058067569091, 0.4902903378454599, 10.0, -0.4902903378454603, 0.098058067569091, -0.987744983947536, 0.0980580675690928, 0.4902903378454603, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.0980580675690928, 0.4902903378454599, 10.0, -0.4902903378454599, 0.0980580675690928, -0.987744983947536, 0.4902903378454599, 10.0, 100.97143868952186, 0.098058067569092, 0.098058067569092, 100.97143868952186, 0.098058067569092, -2.4e-15, 0.098058067569092, 100.97143868952186, 0.0980580675690921, -4.9e-15, 0.0980580675690921, 100.97143868952186, 0.0980580675690919, -7.3e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -9.8e-15, 0.0980580675690919, 100.97143868952186, 0.0980580675690923, -1.22e-14, 0.0980580675690923, 100.97143868952186, 0.0980580675690919, -1.47e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.71e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -1.96e-14, 0.0980580675690919, 100.97143868952186, 0.0980580675690919, -2.2e-14], 66, 66)

I entered your first example on my machine and I obtain the right answer. It means that the problem must be related to the version of Julia installed, which will make it more difficult to solve. I have Julia 1.6.3.

However, I have also a problem for the second problem. It is an accuracy issue that will try to solve.

stevengj commented 9 months ago

Thanks for looking into this, @smataigne!

smataigne commented 9 months ago

I found why the second example fails. It is the exact same reason as in #118 . We reach a corner case of Wilkinson shifts which leaves the matrix unchanged. This is exactly the case described in the file corner_case.pdf attached. A solution is provided in the file strategy.pdf attached. We have to choose special shifts in this case (QR without shifts does not solve the problem). corner_case.pdf strategy.pdf

smataigne commented 9 months ago

Concretely, for the second example, (just like in #118) there is a moment where the tail of the tridiagonal form becomes

julia> v  =[2.0000000000122116, 1.004623745816162e-11, -2.000000000008927]
julia> A = SkewHermTridiagonal(v)
4×4 SkewHermTridiagonal{Float64, Vector{Float64}, Nothing}:
  ⋅   -2.0            ⋅            ⋅
 2.0    ⋅           -1.00462e-11   ⋅
  ⋅    1.00462e-11    ⋅           2.0
  ⋅     ⋅           -2.0           ⋅

With the Wilkinson shifts, this matrix converges so slowly to the real Schur form that the accuracy of the problem will be lost much before reaching the convergence.

cvsvensson commented 8 months ago

Thanks for the fix! However, now one can encounter LAPACKException(22) as in https://github.com/JuliaLinearAlgebra/SkewLinearAlgebra.jl/issues/64. In the attached txt file, there's a 92x92 matrix with this issue. The standard eigen works, but the SkewHermitian version uses a SymTridiagonal version that throws the exception.

eigen(Matrix(sp)) #Works
eigen(SkewHermitian(Matrix(sp))) #LAPACKException(22)

sp.txt