RalphAS / GenericSchur.jl

Julia package for Schur decomposition of matrices with generic element types
Other
27 stars 4 forks source link

should these eigenvecs differ this much? #1

Closed JeffreySarnoff closed 4 years ago

JeffreySarnoff commented 5 years ago
julia> using GenericSchur, LinearAlgebra

julia> setprecision(BigFloat,53)
53

julia> m = reshape([1.0, 12.0, 8.0, 5.0],2,2)
2×2 Array{Float64,2}:
  1.0  8.0
 12.0  5.0

julia> b=BigFloat.(m)
2×2 Array{BigFloat,2}:
  1.0  8.0
 12.0  5.0

julia> eigen(m)
Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}
eigenvalues:
2-element Array{Float64,1}:
 -7.0
 13.0
eigenvectors:
2×2 Array{Float64,2}:
 -0.707107  -0.5547
  0.707107  -0.83205

julia> eigen(b)
Eigen{Complex{BigFloat},Complex{BigFloat},Array{Complex{BigFloat},2},Array{Complex{BigFloat},1}}
eigenvalues:
2-element Array{Complex{BigFloat},1}:
 -7.0 + 0.0im
 13.0 + 0.0im
eigenvectors:
2×2 Array{Complex{BigFloat},2}:
 -1.0+0.0im  0.666667+0.0im
  1.0+0.0im       1.0+0.0im
RalphAS commented 5 years ago

In principle there is an arbitrary scale factor for eigenvectors, and the LinearAlgebra developers didn't document the selection criterion, but consistency would be preferable. I should have followed the LAPACK convention for normalization (unit Euclidean norm), since that is inherited by Julia for BLAS types. I'll plan on doing that in an upcoming release. Note that an ambiguity of phase factor will likely perdure. Thanks for bringing it to my attention.

kagalenko-m-b commented 4 years ago

Not normalizing eigenvectors results in breakdown of the relation

E,Q = svd(A)
A \approx E*Q*E'

It is a serious bug.

RalphAS commented 4 years ago

Can you clarify? I'm not aware of a version of svd which would invoke methods in this package. But

E,Q = eigen(A)
A \approx Q*Diagonal(E)*inv(Q)

should work with the existing (lack of) normalization, as long as Q is reasonably well-conditioned.

(I don't mean to be dismissive; I do intend to standardize the normalization.)

kagalenko-m-b commented 4 years ago

Sorry, svd() was a typo.

Note, that in your version you use inv(Q), because the Q is not orthogonal. It is generally expected that the eigenmatrx is orthogonal. Also, if that expression used in program, inv() is going to call eigen() or svd().

That's not a "choice of normalization", it's a serious bug unfixed for a year that broke my program and cost a good part of day to track down.

JeffreySarnoff commented 4 years ago

Leading with ignorance and hoping to provide a coherent way forward ... Here is an answer to how matlab normalizes eigenvectors.

kagalenko-m-b commented 4 years ago

Matlab relies on LAPACK for most of its linear algebra, so it is not surprising that it normalizes the eigenvectors the same way. So does scipy/numpy. Mpmath (python arbitrary precision library) does not use LAPACK, but also normalizes the eigenvectors to unity.

JeffreySarnoff commented 4 years ago

Normalize usually implies to unity for this class of numerical algorithms.

kagalenko-m-b commented 4 years ago

And that's what Matlab's eig() does. In case of generalized eigenvalue problem it uses non-euclidean norm.

JeffreySarnoff commented 4 years ago

can you make a PR to help, so we resolve this?

kagalenko-m-b commented 4 years ago

I am insufficiently familiar with GenericSchur internals, sorry. In my own program I just normalize them myself

for k = 1:size(Q)[2]
  Q[:, k] /= norm(Q[:, k])
end

That will grow superfluous once the eigenvalue normalization is fixed.

RalphAS commented 4 years ago

The requested normalization convention is in v0.4.0, and is now documented and tested.

JeffreySarnoff commented 4 years ago

(thanks Ralph -- its a hit, a long ball .. and the crowd goes wild!)

kagalenko-m-b commented 4 years ago

Thank you, Ralph, that is very helpful!