search

JuliaLinearAlgebra / IterativeSolvers.jl

Iterative algorithms for solving linear systems, eigensystems, and singular value problems
MIT License
394 stars 106 forks source link

LOBPCG constraint requirements #294

Open doddgray opened 3 years ago

doddgray commented 3 years ago

IterativeSolvers.jl's LOBPCG methods for partial eigen-decomposition provide a keyword argument interface for constraint operators/matrices. The documentation for this argument currently says

C: constraint to deflate the residual and solution vectors orthogonal to a subspace; must overload mul!

The requirements for constraint maps in the current implementation are actually much more strict than this (and much more strict than for the target operator and preconditioner).

I ran into this issue while trying to use constraints of type LinearMaps.FunctionMap with IterativeSolvers.lobpcg. During construction of callable Constraint structs, constraint operators pass through various pre-processing steps that are incompatible with general LinearMaps types (similar and Hermitian). More broadly, the algorithm for constraints seems to require setindex! via view and slicing, as well as ldiv!, see update! and methods of Constraint:

https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/ae01dfe228138d0db9164eb343252150f6c0bfc6/src/lobpcg.jl#L144-L224

Could somebody more knowledgeable explain the actual requirements and suggested types for constraint operators? I think the main appeal of IterativeSolvers.lobpcg is compatibility with "matrix-free"/implicit linear operator types, so the currently undocumented constraint interface requirements could lead to confusion.

It would be great to allow more general constraint functions & operators if it wouldn't be too much work. IIUC, these extra requirements for constraint operators arise because nev (# eigenvalues/vectors solved) extra columns are glued onto the constraint matrix and mutated during each iteration. Could the same processing be achieved without constraining the constraints (hehe), just by creating a separate array? I would be glad to take a crack at it if a PR is desired, although this solver seems to be stuck mid-overhaul at the moment (#247).

doddgray commented 3 years ago

apologies, I think I misunderstood the intended input style for constraints here. Comparing with the docs for scipy.sparse.linalg.lobpcg (here) it looks like constraints is expecting an N x n_constraints array (just a few concatenated column vectors), where the target operator is N x N and n_constraints is determined by the number of independent constraints you want to apply. Is that right?

It would be helpful to see an example of this being using correctly.

mohamed82008 commented 3 years ago

Here is an example https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/master/test/lobpcg.jl#L220.

mohamed82008 commented 3 years ago

Comparing with the docs for scipy.sparse.linalg.lobpcg (here) it looks like constraints is expecting an N x n_constraints array (just a few concatenated column vectors), where the target operator is N x N and n_constraints is determined by the number of independent constraints you want to apply. Is that right?

Yes that's correct.