Hi all, I’ve release ArnoldiMethod v0.3.5 and v0.4.0 recently, after a relatively long break from writing Julia code.

See GitHub - JuliaLinearAlgebra/ArnoldiMethod.jl: The Arnoldi Method with Krylov-Schur restart, natively in Julia.

ArnoldiMethod.jl can be used to compute eigenvalues and eigenvectors (or Schur vectors) of large and sparse matrices. If you’ve used MATLAB, you may know the `eigs`

function, which does something similar.

The package implements the Krylov-Schur algorithm, which is a slight improvement over the Implicitly Restarted Arnoldi Method you may know from ARPACK and MATLAB’s `eigs`

function.

The ArnoldiMethod.jl package should be much more stable than Arpack.jl, and slightly more stable than KrylovKit.jl’s `eigsolve`

/ `schursolve`

functions. It should also have better performance than KrylovKit.jl (mostly because it pre-allocates and works in-place as much as possible). If not, open an issue.

Finally, ArnoldiMethod.jl is I believe the only implementation of a partial Schur decomposition where *everything* is implemented natively in Julia, meaning that it works on arbitrary number types (I’ve only tested DoubeFloats.jl, BigFloat).

There are no algorithmic differences between v0.4 and v0.3. It’s just that v0.4 exports fewer names.

What’s new in both v0.3 and v0.4:

- You can provide an initial guess
`partialschur(A; v1=...)`

, which is the vector that induces the Krylov subspace.
- You can pre-allocate large arrays and temporaries using the
`ArnoldiWorkspace`

struct, and run the algorithm in-place using `partialschur!`

- You can reuse an existing partial Schur decomposition in an
`ArnoldiWorkspace`

and continue to compute additional eigenvalues and Schur vectors.
- Stability of restart is improved: converged eigenvectors with eigenvalues far away from the target are properly purged (= removed from the Krylov subspace) so that it’s unlikely the algorithm has to do computations on an array with relatively very big and really tiny numbers, which otherwise leads to loss of precision. Also, in one crucial bit where large and tiny numbers are expected, Given’s rotations are now used instead of Householder reflections.
- Fixed a few bugs, mostly related to locking of converged vectors

What’s still *not* in v0.3 and v0.4 after all these years:

A utility function like `eigs`

that directly returns eigenvalues and eigenvectors. Right now ArnoldiMethod.jl is still a two-step method: first compute a Schur decomposition, then transform it into an eigendecomposition. I’m open for PRs.

Give it a try!

28 Likes

Concerning the desirable extensions: for the “vibration” problem (symmetric generalized stencil, looking for smallest eigenvalues), the package GitHub - PetrKryslUCSD/VibrationGEPHelpers.jl: Generalized Eigenvalue Problem helper functions wraps ArnoldiMethod.jl (and other similar packages), and so provides nearly uniform interface similar to eigs.

4 Likes

I’ve been using this package with great success to compute the stationary distribution for markov transition matrices. As promised, it delivers a modest performance gain relative to KrylovKit and a large gain relative to Arpack.

@stabbles, could you explain the extent to which it takes advantage of multithreading? Does it use the number of threads specified by Julia (e.g.,`threads.nthreads()`

)? …or does it use the the number of BLAS threads?

I find that it is not achieving very high rates of CPU utilization when doing the partial schur decomposition. My transition matrix that I’m working with is a sparse matrix of about 25,000 x 25,000 with about 60% fill.

I’m not @stabbles, but I took a look at the code and didn’t see any julia native use of threading (no use of the `@threads`

or `@spawn`

macros, and no dependence on other packages providing multithreading APIs). However, it uses `LinearAlgebra.mul!`

and friends, which are multithreaded when used on array types that dispatch to BLAS (provided `BLAS.get_num_threads() > 1`

).

That’s not a very sparse matrix. You’d probably have better luck using a regular dense `Matrix`

.

3 Likes

Thanks. Very insightful. Unfortunately memory constraints preclude me from working with the transition matrix in dense form.

Any chance you can access a machine with more memory? Or use Float32? We’re talking about less than a factor of 2 here; sparse matrices are typically used when they cut down on memory by orders of magnitude and can be brutally slow for practically dense matrices like this.

Actually, I misstated the size and fill rate. Apollogies for that. Size is actually about 580,000 x 580,000. And the fill is less than 1%.

My eyeballing the fill from the `spy()`

plot was apparently really inaccurate!

2 Likes

Most of the time should be taken by matrix-vector products for a problem of this size (profile to check). So the challenge is to make sure that the full memory bandwidth is being used; one often needs multiple threads to get that. Try the MKL sparse routines.

A GPU sparse matrix library might provide dramatic speed improvements.

For the size you mention, you might need a system with multiple GPUs and a GPU wizard to help design a partitioned product method for your matrix structure. It would be interesting to verify that ArnoldiMethod can take advantage of such things.

Does this package work with arrays on GPU?