LOBPCG unstable, PosDefExceptions

[jagot]

May be related to #223.

Some background info: I'm trying the solve the Hartree–Fock equations, which constitute a system of coupled integro-differential eigenproblems, which are solved iteratively:

where P is the eigenvector to be optimized to minimize an overall energy, ε its associated eigenvalue, Υ is a potential (diagonal matrix) formed by the other coupled vectors, and X is an integral operator (full matrix), also formed by the other coupled vectors. The algorithm is to

1. Given a set of initial vectors, form the operators Y and X,
2. Solve all the eigenproblems, giving new estimates,
3. Goto 1. until converged

So far, I have been using ArnoldiMethod.jl which mostly works, but has problems with condition numbers when decreasing the grid step size.

Using the following script

using Pkg
Pkg.activate(".")

using IterativeSolvers
using ArnoldiMethod
using LinearAlgebra
using BSON

function test_diagonalization(filename, orbital, method, tol=1e-8)
    println("Optimizing $(orbital) using $(method)")
    data = BSON.load(filename)[orbital]
    A,P,x = data[:A],data[:P],data[:x]
    if method == :arnoldi
        schurQR,history = partialschur(A, nev=1, tol=tol, which=SR())
        println(history)
        λ = schurQR.R[1]
        println("λ = $λ")
        λ
    else
        try
            res = lobpcg(A, false, x, 1, P=factorize(Symmetric(P)), tol=tol, maxiter=1000)
            show(res)
            res.λ[1]
        catch PosDefException
            "PosDefExc"
        end
    end
end

filenames = [
    "helium.bson", # 2 equations
    "beryllium.bson", # 4 equations
    "neon.bson" # 10 equations
]

orbitals = [
    "1s₀α", "1s₀β", # Equations 1-2
    "2s₀α", "2s₀β", # Equations 3-4
    "2p₋₁α", "2p₋₁β", "2p₀α", "2p₀β", "2p₁α", "2p₁β" # Equations 5-10
]

methods = [
    :arnoldi,
    :lobpcg
]

test_diagonalization(filenames[1], orbitals[1], methods[1])

and the attached file lobpcg-data.zip, which contains the operators formed in the first step above, for the first iteration, I get the following matrix:

┌───────────┬─────────┬───────────────────────┬───────────────────────┬───────────────────────┬───────────────────────┬───────────────────┬───────────────────┬───────────────────┬───────────────────┬───────────────────┬───────────────────┐
│   Element │  Method │                  1s₀α │                  1s₀β │                  2s₀α │                  2s₀β │             2p₋₁α │             2p₋₁β │              2p₀α │              2p₀β │              2p₁α │              2p₁β │
├───────────┼─────────┼───────────────────────┼───────────────────────┼───────────────────────┼───────────────────────┼───────────────────┼───────────────────┼───────────────────┼───────────────────┼───────────────────┼───────────────────┤
│    helium │ arnoldi │   -0.8207003606948218 │   -0.8207003606952347 │                       │                       │                   │                   │                   │                   │                   │                   │
│    helium │  lobpcg │   -0.8207003606950175 │   -0.8207003606950175 │                       │                       │                   │                   │                   │                   │                   │                   │
│ beryllium │ arnoldi │   -3.9520775411349813 │   -3.9520775411552345 │   -0.2454902061332506 │  -0.24549020613648817 │                   │                   │                   │                   │                   │                   │
│ beryllium │  lobpcg │ -4.364707164462907e14 │             PosDefExc │             PosDefExc │             PosDefExc │                   │                   │                   │                   │                   │                   │
│      neon │ arnoldi │    -102.3242323345127 │   -102.32423233446154 │     -463.004359839325 │    -463.0043598393252 │ 4.885303298697446 │ 4.885303298697389 │ 4.884194503733136 │ 4.884194503733244 │ 4.885303298697248 │ 4.885303298697376 │
│      neon │  lobpcg │ -3.817247259165677e16 │ -4.091724928996284e16 │ -5.881061779982185e13 │ -5.881061779982185e13 │ 4.885303298697289 │ 4.885303298697285 │ 4.884194503733101 │ 4.884194503733094 │ 4.885303298697291 │ 4.885303298697287 │
└───────────┴─────────┴───────────────────────┴───────────────────────┴───────────────────────┴───────────────────────┴───────────────────┴───────────────────┴───────────────────┴───────────────────┴───────────────────┴───────────────────┘

Here we can see that Arnoldi gives reasonable eigenvalues for all elements/equations, whereas LOBPCG struggles for anything but helium.

Am I using LOBPCG wrongly, or is it simply not suited for this problem?

[jagot]

I discovered that my matrices are, in fact, not symmetric, which they should be. I need to investigate why this is the case. Maybe this issue can be closed.

[mohamed82008]

If they are not symmetric up to some minor error, try wrapping it with the Symmetric wrapper. Check norm(A - A').

[jagot]

I figured out (had forgotten) why the matrices are not symmetric; to maintain orthogonality between the different vectors, every time I apply the operator to a vector, I project out the component parallel to the other vectors (this is equivalent to Lagrange multipliers). This forces e.g. 2s₀α to be different from 1s₀α, however, this ruins the symmetry of the matrix.

[jagot]

Indeed, applying the projection operators symmetrically on each side of the linear operator makes the overall operator symmetric again, and LOBPCG does not crash anymore. However, when the self-consistent iteration approaches the correct eigenpair, LOBPCG takes a big jump elsewhere and the self-consistent iteration either starts converging again or diverges.

[mohamed82008]

You can also use the constraint functionality in LOBPCG designed for this purpose. Use the C keyword argument, where C is a matrix whose columns are going to be B-orthogonal to the solution eigenvector.

[mohamed82008]

In your case B is the identity.