LOBPCG constraint requirements

[doddgray]

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:

IterativeSolvers.jl/src/lobpcg.jl

Lines 144 to 224 in ae01dfe

	mutable struct Constraint{T, TVorM<:Union{AbstractVector{T}, AbstractMatrix{T}}, TM<:AbstractMatrix{T}, TC}
	Y::TVorM
	BY::TVorM
	gram_chol::TC
	gramYBV::TM # to be used in view
	tmp::TM # to be used in view
	end
	struct BWrapper end
	struct NotBWrapper end
	Constraint(::Nothing, B, X) = Constraint(nothing, B, X, BWrapper())
	Constraint(::Nothing, B, X, ::NotBWrapper) = Constraint(nothing, B, X, BWrapper())
	function Constraint(::Nothing, B, X, ::BWrapper)
	return Constraint{Nothing, Matrix{Nothing}, Matrix{Nothing}, Nothing}(Matrix{Nothing}(undef,0,0), Matrix{Nothing}(undef,0,0), nothing, Matrix{Nothing}(undef,0,0), Matrix{Nothing}(undef,0,0))
	end
	Constraint(Y, B, X) = Constraint(Y, B, X, BWrapper())
	function Constraint(Y, B, X, ::BWrapper)
	T = eltype(X)
	if B isa Nothing
	BY = Y
	else
	BY = similar(Y)
	mul!(BY, B, Y)
	end
	return Constraint(Y, BY, X, NotBWrapper())
	end
	function Constraint(Y, BY, X, ::NotBWrapper)
	T = eltype(X)
	if Y isa SubArray
	gramYBY = @view Matrix{T}(I, size(Y.parent, 2), size(Y.parent, 2))[1:size(Y, 2), 1:size(Y, 2)]
	mul!(gramYBY, adjoint(Y), BY)
	gramYBV = @view zeros(T, size(Y.parent, 2), size(X, 2))[1:size(Y, 2), :]
	else
	gramYBY = adjoint(Y)*BY
	gramYBV = zeros(T, size(Y, 2), size(X, 2))
	end
	realdiag!(gramYBY)
	gramYBY_chol = cholesky!(Hermitian(gramYBY))
	tmp = deepcopy(gramYBV)
	return Constraint{eltype(Y), typeof(Y), typeof(gramYBV), typeof(gramYBY_chol)}(Y, BY, gramYBY_chol, gramYBV, tmp)
	end
	function update!(c::Constraint, X, BX)
	sizeY = size(c.Y, 2)
	sizeX = size(X, 2)
	c.Y.parent[:, sizeY+1:sizeY+sizeX] .= X
	if X !== BX
	c.BY.parent[:, sizeY+1:sizeY+sizeX] .= BX
	end
	sizeY += sizeX
	Y = @view c.Y.parent[:, 1:sizeY]
	BY = @view c.BY.parent[:, 1:sizeY]
	c.Y = Y
	c.BY = BY
	gram_chol = c.gram_chol
	new_factors = @view gram_chol.factors.parent[1:sizeY, 1:sizeY]
	c.gram_chol = typeof(gram_chol)(new_factors, gram_chol.uplo, convert(LinearAlgebra.BlasInt, 0))
	c.gramYBV = @view c.gramYBV.parent[1:sizeY, :]
	c.tmp = @view c.tmp.parent[1:sizeY, :]
	return c
	end
	function (constr!::Constraint{Nothing})(X, X_temp)
	nothing
	end
	function (constr!::Constraint)(X, X_temp)
	@views if size(constr!.Y, 2) > 0
	sizeX = size(X, 2)
	sizeY = size(constr!.Y, 2)
	gramYBV_view = constr!.gramYBV[1:sizeY, 1:sizeX]
	mul!(gramYBV_view, adjoint(constr!.BY), X)
	tmp_view = constr!.tmp[1:sizeY, 1:sizeX]
	ldiv!(tmp_view, constr!.gram_chol, gramYBV_view)
	mul!(X_temp, constr!.Y, tmp_view)
	@inbounds X .= X .- X_temp
	end
	nothing
	end

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]

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]

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

[mohamed82008]

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.