Function to return any eigenvector

[ParadaCarleton]

Sometimes it's useful to compute just a single eigenvector (e.g. by the power method), and any eigenvector will do, e.g. I have an irreducible Markov chain, where the eigenvector is unique. I couldn't find any functions for this; should we add one to LinearAlgebra?

[oscardssmith]

see https://github.com/JuliaLang/julia/pull/49487 which was my attempt to write something similar to this for making opnorm2 faster.

[cafaxo]

Packages that offer partial eigendecompositions:

• Arpack.jl
• ArnoldiMethod.jl
• KrylovKit.jl
• IterativeSolvers.jl

I am not sure that such iterative methods fit into LinearAlgebra when we already have specialized packages for this.

[ParadaCarleton]

Packages that offer partial eigendecompositions:

• Arpack.jl
• ArnoldiMethod.jl
• KrylovKit.jl
• IterativeSolvers.jl

I am not sure that such iterative methods fit into LinearAlgebra when we already have specialized packages for this.

Definitely true for some of these! Fitting the entirety of KrylovKit.jl into LinearAlgebra (including the GPUArrays dependency) would definitely not be the best idea :sweat_smile:

But in some cases, I think they'd fit in perfectly. The power method seems basic, small, and well-established enough that I think it'd be a good fit for LinearAlgebra; I was actually very surprised to learn it isn't included.

[stevengj]

The power method seems basic, small, and well-established

Also mostly obsolete. Krylov methods are generally strictly better, even with restarting.

[ParadaCarleton]

The power method seems basic, small, and well-established

Also mostly obsolete. Krylov methods are generally strictly better, even with restarting.

Huh, that's surprising, I thought the power method was still standard for the largest eigenvalue problem.

[aravindh-krishnamoorthy]

Sometimes it's useful to compute just a single eigenvector (e.g. by the power method), and any eigenvector will do, e.g. I have an irreducible Markov chain, where the eigenvector is unique. I couldn't find any functions for this; should we add one to LinearAlgebra?

If your matrix is symmetric, kindly take a look at eigen(A::RealHermSymComplexHerm, irange::UnitRange) in symmetriceigen.jl, which can optionally return a single eigenvalue and associated eigenvector. It uses LAPACK's XSYEVR, which is quite efficient.

Huh, that's surprising, I thought the power method was still standard for the largest eigenvalue problem.

If deemed sufficient, a basic version of the power algorithm can be easily implemented by users themselves. Please take a look at the code in https://github.com/JuliaLang/julia/pull/49487 (mentioned in the earlier comments) for an example implementation. I feel it might be a great idea to add this algorithm to the Julia documentation.

[ParadaCarleton]

If deemed sufficient, a basic version of the power algorithm can be easily implemented by users themselves.

Definitely can be, but didn't @stevengj just say it's obsolete? The largest eigenvalue problem is common enough that I'd think LinearAlgebra should include it, but I'm not sure what algorithm it should use.

[aravindh-krishnamoorthy]

If deemed sufficient, a basic version of the power algorithm can be easily implemented by users themselves.

Definitely can be, but didn't @stevengj just say it's obsolete?

It's a first algorithm taught in numerical linear algebra classes. The advantages and disadvantages are well known, please see, e.g., the Wiki page and references therein. The method works well when the eigenvalues are unique and sufficiently separated. However, since it is only beneficial for a relatively narrow class of matrices and non-convergent (or very slowly convergent) for others, it may not be suitable for stdlib, especially when several better alternatives targeting specialised applications are known. Hardly ever used in serious business, in my experience so far.

Edit: R has a slightly modified version with scaling in the matlib package aimed at teaching. So, may be worth looking at if there's sufficient interest. But, perhaps better done in a specialised package rather than stdlib.

[ParadaCarleton]

Yeah, if it's only really useful for educational reasons, it should go in a package. But the largest eigenvalue problem is common enough that I'd guess there should be a good default option in LinearAlgebra, or people might just reach for the first algorithm they know (usually, power iteration).

[stevengj]

We used to have an eigs function in base (to return extremal eigenvalues/vectors), but it was moved to a package: https://github.com/JuliaLang/julia/pull/27616

I don't really see us reversing this decision and putting Krylov iterative methods back into Base. There are now several good packages for this. The algorithms are qualitatively rather different from the dense linear algebra in LinearAlgebra, or even from the sparse-direct methods in SparseArrays.