Jutho / KrylovKit.jl

Krylov methods for linear problems, eigenvalues, singular values and matrix functions
Other
284 stars 37 forks source link

eigsolve with x0 as all eigenvectors #92

Closed GabrielPonte closed 2 weeks ago

GabrielPonte commented 2 weeks ago

Hi, I would like to know if it is possible to input all the initial eigenvectors as an input.

For example: n=5; A = rand(n,n); A = A*A'; vals,vecs,info = eigsolve(A); Then, I would like to get the new eigendecomposition of a modified problem giving initial eigenvectors to compute the solution faster B = rand(n,n); B = 1e-2*B*B'; vals2,vecs2,info2 = eigsolve(A+B,vecs); but vecs2 is not equal to n vectors of length n.

Is it possible to get vecs2 to be the same size as vecs1?

Jutho commented 2 weeks ago

The starting vector, i.e. the second argument to eigsolve, can only be one vector. What you are observing here, is the fact that if vecs is a Vector of vectors, then A * vecs actually still works and behaves as a linear operator on vecs, which should be thought of as a single vector with a nested / block structure. KrylovKit.jl is written to support very general vector types, and thus it supports such nested vectors, provided the linear operator can deal with them.

What you get in this case is that A * vecs produces a new nested vector (i.e. a Vector{Vector{...}} instance) with subvector entries given by

outvecs = A vecs => outvecs[i] = sum( A[i,j] vecs[j] for j=...)

This is a perfect linear operator (it is equivalent to kron(A, I) with I the identity operator acting on the vectors in vecs itself), and is being diagonalised.

Krylov methods can only take one starting vector. If you want to include the information of multiple vectors, the best way to do so is to start from a linear combination of those vectors. So what you could do is

vals2, vecs2, info2 = eigsolve(A + B, rand(length(vecs))' * vecs)

GabrielPonte commented 2 weeks ago

Thanks!