Generalize 2-arg eigen towards AbstractMatrix

[dkarrasch]

Fixes JuliaLang/LinearAlgebra.jl#953, and cleans up symmetriceigen.jl code a bit (without functional changes).

[pcjentsch]

As I mentioned in my reply in JuliaLang/LinearAlgebra.jl#953, with this implementation,

eigen(A::Diagonal,B::Diagonal) != eigen(Matrix(A),B::Diagonal) != eigen(Matrix(A),Matrix(B))

despite these all representing the same matrices, but I think this would hold for the other matrix types. They're all correct of course, just represented differently. Does that matter?

[dkarrasch]

As long as you get a correct answer, I'm not sure a user would bother whether they would get the same result if they chose a different (i.e. dense) representation. Provided correctness, the next goal is performance, so simply densifying doesn't sound like an appealing option to me.

[jishnub]

While this is being addressed, I wonder if specializing cases like eigen(A::AbstractMatrix, B::Symmetric) are worth it, where B is positive-definite? Such cases arise at times (I am working with one such case in a fluid dynamics application), and this may be transformed to a standard eigenvalue problem by performing a Cholesky decomposition of B, which would lead to a significant speedup over evaluating eigen(A, Matrix(B)).

An example for 1000x1000 matrices:

julia> A = rand(1000,1000);

julia> X = rand(1000,1000);

julia> S2 = Symmetric(X'*X);

julia> @time λ1, v1 = eigen(A, Matrix(S2));
  4.999689 seconds (25 allocations: 45.953 MiB)

julia> @time λ2, v2 = begin
           C = cholesky(S2)
           λ, v = eigen(C.L\(A/C.U))
           λ, C.U\v
       end;
  0.446539 seconds (41 allocations: 94.568 MiB)

julia> real(λ1) ≈ real(λ2)
true

julia> λ2s = λ2[[argmin(abs.(λ1 .- λ)) for λ in λ2]];

julia> λ1 ≈ λ2s
true

Currently, this method errors,

julia> eigen(A, S2);
ERROR: MethodError: no method matching eigen!(::Matrix{Float64}, ::Symmetric{Float64, Matrix{Float64}})

so limiting it to the special case of positive definite matrices might not be too bad, I guess?

I can make a separate PR for this, if this seems reasonable. One concern with this is that the numerical precision of eigen(A, B::Symmetric) would differ from that of eigen(A, Matrix(B::Symmetric)) (I'm not sure which one is better).

[dkarrasch]

We have

julia/stdlib/LinearAlgebra/src/symmetriceigen.jl

Lines 170 to 175 in 276af84

	function eigen!(A::RealHermSymComplexHerm{T,S}, B::AbstractMatrix{T}; sortby::Union{Function,Nothing}=nothing) where {T<:Number,S<:StridedMatrix}
	U = cholesky(B).U
	vals, w = eigen!(UtiAUi!(A, U))
	vecs = U \ w
	GeneralizedEigen(sorteig!(vals, vecs, sortby)...)
	end

which does exactly what you propose, but the method signature don't lead there for your example. I'll try to come up with some "promotion mechanism" that would help your case.

[dkarrasch]

Actually, can you help me benchmark the UtiAUi thing and compare to rdiv!(ldiv!(U', parent(A)), U)? I must have made a mistake, but on my machine UtiAUi! is not allocation-free, and performance is just horrible.

[dkarrasch]

Tests pass locally. It would be great if people could try this on their favorite use cases, and look at precision and performance.

[antoine-levitt]

I'm incompetent to review the dispatch logic (which feels a bit convoluted, but then most of LinearAlgebra is...), but otherwise looks good to me. I've been wanting eigen for more general types for a while so this is great

[antoine-levitt]

In particular, it looks like it does eigen of tridiagonal matrices, which is awesome (this comes up all the time for quick & dirty 1D PDE examples)

[dkarrasch]

In particular, it looks like it does eigen of tridiagonal matrices, which is awesome (this comes up all the time for quick & dirty 1D PDE examples)

Well, but it's not specialized. It densifies the Tridiagonal lhs, but keeps diagonal rhs.

[antoine-levitt]

Yeah I don't mind the fact it's not specialized, it's nice that it works at all!