implement sparse solver

[williamesamuelson]

Added sparse RQ solver from #1.

However, I am a bit confused about the RQ_EIG solver. Do we allow $b$ to be a matrix?

[cvsvensson]

b can be a matrix, as in the section Extension in https://www.cis.upenn.edu/~jshi/papers/supplement_nips2006.pdf. We return one eigenvector for each column of b.

For weak majoranas, that means we can solve for both majoranas at once (the two columns of b are the two different right hand sides), while reusing the factorization of CC'.

[williamesamuelson]

I see, thanks!

[cvsvensson]

Perhaps https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl is preferable over Arpack.

[cvsvensson]

Perhaps https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl is preferable over Arpack.

Actually, our problem is hermitian, so the Lanczos algorithm is appropriate. Arpack.jl and KrylovKit.jl has this. I'd suggest trying the latter.

[williamesamuelson]

we decided that we should return a vector from solve and that $b$ can be given as a single vector or as a Span struct.

todo

• [x] add a struct Span that can be put in $b$ in ConstrainedRQProblem
• [x] replace _select_vecs(..., b::AbstractMatrix) with _select_vecs(..., b::Span). No, one _select_vecs should work.
• [x] replace eigen with eigsol (?) from KrylovKit for sparse solver
• [x] add tests with random matrices

[williamesamuelson]

There seems to be something going on in select_vectors causing RQ_EIG and RQ_SPARSE return different vectors.

function _select_vectors(eigvecs, P, eps=DEFAULT_SELECTVEC_EPS)
    vecs = eachcol(eigvecs)
    v = similar(first(vecs))
    f = x -> any(>(eps) ∘ abs, mul!(v, P, x))
    i1 = findfirst(f, vecs)
    return vecs[i1]
end

Could we instead return the vector corresponding to the first eigenvalue>epsilon?

Related to this, it seems like it would be easier to solve the maximization problem instead. In that case, with KrylovKit, we just have to calculate one eigenvector without comparing a bunch of eigenvalues to epsilon.

I guess that would be fine in the WeakMajorana problem if we adjust the penalty matrix.

Actually, it would be nice to implement both maximize and minimize

[codecov[bot]]

Welcome to Codecov :tada:

Once you merge this PR into your default branch, you're all set! Codecov will compare coverage reports and display results in all future pull requests.

Thanks for integrating Codecov - We've got you covered :open_umbrella:

[williamesamuelson]

A few questions/issues:

• How should we handle sparseness? Should we depend on SparseArrays and convert where we need to keep sparsity? or should the user supply sparse matrices ($C$ and $b$) in the problem? Now I convert into sparse in RQ_SPARSE.
• ldiv! is not supported for Sparse Cholesky factorizations, see e.g., https://github.com/JuliaSparse/SparseArrays.jl/issues/123. Should we change the factorization of $CC'$ or the solver to something else than ldiv! in the InverseMap? Now I changed the factorization to LU.
• If solving just one eigenvalue didn't provide the correct solution, we need to request more solutions in Krylov eigsolve. But now we request too many compared to krylovdim. I guess one solution would be to have some default values/let the user decide and if it is not enough, we let the user know that a solution wasn't found.

[cvsvensson]

• User should supply sparse arrays. We may need SparseArrays dependency for dispatch and for selecting default solver.
• I guess options are lu+ldiv! or chol+. Pick whatever works for now. One can benchmark the two choices later.
• It might make sense to expand krylovdim if the user didn't set it manually and nev is larger.