Fusing LOBPCG with Laplacians.KMPLapSolver() preconditioner from Laplacians.jl

[bodhi91]

This may be a novice question that I am putting forward. I am fairly new to Julia. I am solving a generalized eigenvalue problem in Julia leveraging the LOBPCG solver from the IterativeSolvers.jl package for Laplacian matrices of size ranging between 40Kx40K to 500Kx500K. The package provides a way to add a preconditioner and the documentation mentions that 'P' the preconditioner must overload the ldiv! function. The preconditioner of my choice for this problem is taken from the Laplacians.jl package. I intend to use the Laplacians.KMPLapSolver, based on the paper "Approaching optimality for solving SDD systems" by Koutis et-al as my preconditioner. The KMPLapSolver returns a function which needs to be translated to a format compatible with LOBPCG.

Can someone provide some guidance as to how to use a different preconditioner than ldiv! for the LOBPCG?

Random.seed!(1)
n = 1000
a = wtedChimera(n, 1) #Adjacency Matrix
la = lap(a) #Laplacian Matrix
sddm = copy(la) 
sddm = sddm + spdiagm(0=>rand(n)/n) #SDDM Matrix
largest = false
nev = 1

# Without Preconditioner
lobpcg(sddm, largest, nev) #works

# With Preconditioner from Laplacians.jl
lapSolver = KMPLapSolver(a)

lobpcg(sddm, P=lapSolver, largest, nev) #Does not work

ERROR: MethodError: no method matching ldiv!(::SubArray{Float64,2,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true}, ::getfield(Laplacians, Symbol("##200#202")){Float64,Int64,Float64,Bool,Array{Int64,1},getfield(Laplacians, Symbol("###200#201#203")){Int64,Array{Int64,1},SparseMatrixCSC{Float64,Int64},getfield(Laplacians, Symbol("##133#135")){Float64,Float64,Float64,Bool,Array{Int64,1},getfield(Laplacians, Symbol("###133#134#136")){Array{Int64,1},getfield(Laplacians, Symbol("##121#123")){getfield(Laplacians, Symbol("###121#122#124")){Array{Int64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SuiteSparse.CHOLMOD.Factor{Float64}}}}}}}, ::SubArray{Float64,2,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true})

Is there any way I can use the preconditioners from the Laplacian package into the LOBPCG? Is it possible to add some documentation regarding the use of preconditioners for LOBPCG with examples?

[mohamed82008]

Try this

using LinearAlgebra
function LinearAlgebra.ldiv!(c::AbstractArray{T}, P::KMPLapSolver, b::AbstractArray{T}) where {T}
    c .= P(b)
    return c
end

before calling lobpcg.

[bodhi91]

Thanks a lot. The following code works:

using LinearAlgebra
using IterativeSolvers
using Laplacians

Random.seed!(1)
n = 1000
a = wtedChimera(n, 1)
b = randn(n)
b = b .- mean(b)
la = lap(a)
sddm = copy(la)
sddm = sddm + spdiagm(0=>rand(n)/n)
largest = false
nev = 1

function LinearAlgebra.ldiv!(c::AbstractArray{T}, P::getfield(Laplacians, Symbol("##200#202")){Float64,Int64,Float64,Bool,Array{Int64,1},getfield(Laplacians, Symbol("###200#201#203")){Int64,Array{Int64,1},SparseMatrixCSC{Float64,Int64},getfield(Laplacians, Symbol("##133#135")){Float64,Float64,Float64,Bool,Array{Int64,1},getfield(Laplacians, Symbol("###133#134#136")){Array{Int64,1},getfield(Laplacians, Symbol("##121#123")){getfield(Laplacians, Symbol("###121#122#124")){Array{Int64,1},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},SuiteSparse.CHOLMOD.Factor{Float64}}}}}}}, b::AbstractArray{T}) where {T}
       c .= P(b)
       return c
end

lobpcg(sddm, largest, nev) #Without Preconditioner
lobpcg(sddm, largest, nev, P=KMPLapSolver(a)) #With Preconditioner

# Results without Preconditioner #

#= 
Results of LOBPCG Algorithm
 * Algorithm: LOBPCG - CholQR
 * λ: [0.0005096320997139823]
 * Residual norm(s): [1.852468775228979e-5]
 * Convergence
   * Iterations: 82
   * Converged: true
   * Iterations limit: 200
=#

# Results with Preconditioner #

#=
Results of LOBPCG Algorithm
 * Algorithm: LOBPCG - CholQR
 * λ: [0.0005096315407290366]
 * Residual norm(s): [2.4995233997278383e-7]
 * Convergence
   * Iterations: 6
   * Converged: true
   * Iterations limit: 200
=#

Can we add this to the document as a future reference for other users?

[dkarrasch]

This is a horrible example. Besides the fact that the code is not working as is, more importantly, KMPLapSolver(a) does not return a concrete type. That's why you have to use this super-specific signature, that does not generalize to another example. You would end up defining ldiv!s for each example! I tried to modify your example slightly, but then I seem to be unable to reproduce your results...

using LinearAlgebra
using SparseArrays
using Random
using Statistics
using IterativeSolvers
using Laplacians

struct LapPreconditioner{T,F}
    f::F
    M::Int
    function LapPreconditioner(A::SparseMatrixCSC{T}) where {T}
        f = KMPLapSolver(A)
        return new{T,typeof(f)}(f, size(A, 1))
    end
end

function LinearAlgebra.ldiv!(c::AbstractVecOrMat{T}, P::LapPreconditioner{T}, b::AbstractVecOrMat{T}) where {T}
    copyto!(c, P.f(b))
end

Random.seed!(1)
n = 1000
a = wted_chimera(n, 1; verbose=false)
b = randn(n)
b = b .- mean(b)
la = lap(a)
sddm = copy(a)
sddm = sddm + spdiagm(0=>rand(n)/n)
largest = false
nev = 1

@time lobpcg(sddm, largest, 1) #Without Preconditioner
@time lobpcg(sddm, largest, 1, P=LapPreconditioner(a)) #With Preconditioner

The results are

Results of LOBPCG Algorithm # without preconditioner
 * Algorithm: LOBPCG - CholQR
 * λ: [-4.586226043892103]
 * Residual norm(s): [1.992905543714895e-5]
 * Convergence
   * Iterations: 56
   * Converged: true
   * Iterations limit: 200

Results of LOBPCG Algorithm # with preconditioner
 * Algorithm: LOBPCG - CholQR
 * λ: [-4.486180608427051]
 * Residual norm(s): [0.41045530978237926]
 * Convergence
   * Iterations: 201
   * Converged: false
   * Iterations limit: 200

In any case, I believe one should make an effort to wrap the result of KMPLapSolver(A) into a nice type, and then dispatch ldiv! on it.

[mohamed82008]

I can provide a PreconditionerFunction in Preconditioners.jl to do this automatically.

[bodhi91]

Thank you Daniel. Even wrapping the result of KMPLapSolver(A) into a nice type does not ensure it will work for all matrices. For example, the example you provided fails for n = 10000 or more. Also, there is a typo in the original code. sddm = copy(la). I will, for now, leverage the AMG.jl package for a suitable preconditioner since that is directly compatible with LOBPCG. Having a PreconditionerFunction in Preconditioners.jl would be great.