Extend KrylovKit.eigsolve by allowing the initial guess x₀ to be a Tensor

[jofrevalles]

Summary

This PR extends the eigsolve function introduced in PR #157 by allowing the initial guess eigenvector x₀ to be a Tensor. Previously, this variable had to be an AbstractVector. We also updated the tests to cover this extension.

Example

julia> using KrylovKit

julia> using Tenet; using LinearAlgebra; using Test

julia> A = rand(ComplexF64, 4, 4)
4×4 Matrix{ComplexF64}: ...

julia> data = (A + A') / 2 # Make it Hermitian

julia> tensor = Tensor(data, (:i, :j))
4×4 Tensor{ComplexF64, 2, Matrix{ComplexF64}}: ..

julia> vals, vecs, info = eigsolve(tensor, Tensor(rand(ComplexF64, 4), (:i,)), 4, :SR, Lanczos(; krylovdim=4, tol=1e-16); left_inds=[:i], right_inds=[:j]) # Perform eigensolve
...

julia> V_matrix = hcat([reshape(parent(vec), :) for vec in vecs]...) # Convert vecs to matrix form for reconstruction
4×4 Matrix{ComplexF64}: ...

julia> D_matrix = Diagonal(vals)
4×4 Diagonal{Float64, Vector{Float64}}:
 -0.573245   ⋅          ⋅         ⋅ 
   ⋅        0.0709279   ⋅         ⋅ 
   ⋅         ⋅         0.430396   ⋅ 
   ⋅         ⋅          ⋅        2.32354

julia> reconstructed_matrix = V_matrix * D_matrix * inv(V_matrix)
4×4 Matrix{ComplexF64}: ...

julia> @test isapprox(reconstructed_matrix, parent(tensor))
Test Passed

[mofeing]

Do you gain something by making x0 a Tensor?

[mofeing]

Why have you removed this?

[jofrevalles]

See the comment from above.

[jofrevalles]

Do you gain something by making x0 a Tensor?

In some iterative algorithms (e.g., DMRG), we may have an initial tensor θ, and we update it with some eigenvector from eisolve. In this case we usually can pass θ to Lanczos, so it makes sense that we have that in a Tensor form.

[mofeing]

In some iterative algorithms (e.g., DMRG), we may have an initial tensor θ, and we update it with some eigenvector from eisolve. In this case we usually can pass θ to Lanczos, so it makes sense that we have that in a Tensor form.

Yes, but as far as I've seen, the result of calling eigsolve now it's not a Tensor right. Maybe it would nice to return a Tensor?

[jofrevalles]

Yes, but as far as I've seen, the result of calling eigsolve now it's not a Tensor right. Maybe it would nice to return a Tensor?

No, the result of eigvecs is a list of Tensors, so this is already working as expected.

[mofeing]

Ok 👍

You just need to deduplicate the code and the PR should be ready.

[jofrevalles]

Ok 👍

You just need to deduplicate the code and the PR should be ready.

Done!