Notebooks · lib/lobpcg

LOBPCG, the lowest eigenpairs of an operator you can only apply

lobpcg(apply_h, precondition, vectors, values, max_steps, nconv, tol) finds the lowest few eigenvalues and eigenvectors of a Hermitian operator H that is given as a function: apply_h(x) returns H applied to each row of x. The matrix is never formed, so H can be a stencil, an FFT sandwich, or anything else that is cheaper to apply than to store.

LOBPCG is short for locally optimal block preconditioned conjugate gradient. Each step searches a small subspace for the vectors with the lowest Rayleigh quotient $x^{H} H x / x^{H} x$ , and the subspace is three blocks of rows: the current approximations, their preconditioned residuals (the direction in which the quotient falls fastest), and the step each approximation took last time. The lowest eigenvalues of H restricted to that subspace (Ritz values) and the vectors that achieve them (Ritz vectors) become the next approximations. The post lobpcg-finds-the-lowest-states explains why this works, on a particle in a box.

The arguments:

vectors is [nb, n] complex, one starting guess per row; it comes back holding the eigenvectors, orthonormal, lowest first. A real problem is a complex one whose imaginary parts are zero.
values is [nb]; it comes back holding the eigenvalues.
precondition(r, x) returns an approximation to $(H - λ)^{- 1} r$ for the residual r of the row x, or r itself when there is nothing better: no_preconditioner is that.
Only the lowest nconv rows are held to tol, a bound on the norm of $H x - λ x$ for a unit x. The rows above them are a buffer: the top of a block converges slowest, because nothing above it in the block keeps it apart from the rest of the spectrum.
It returns the largest residual norm among those nconv rows.

Ask for a tol the arithmetic can reach. On the test problems below the residual stops falling near $1 0^{- 10}$ , and an iteration held at that floor with nothing left to find can come apart: on the 90 by 90 test matrix, asked for $1 0^{- 13}$ , the largest residual was back above $1 0^{- 4}$ by step 120 and above 1 by step 130. The DFT engine asks for $1 0^{- 6}$ .

lobpcg_start and lobpcg_step are the same iteration a step at a time, for a caller that wants to watch it.

Used with import lobpcg. The plane-wave DFT engine dft/pwdft_it uses it for its bands.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 49–472

#: LOBPCG_X, LOBPCG_HX, LOBPCG_P, LOBPCG_HP  where each block sits in a LOBPCG
#:     state, a [4, nb, n] complex array: the current approximations x, H x,
#:     the previous step p, and H p
LOBPCG_X = 0
LOBPCG_HX = 1
LOBPCG_P = 2
LOBPCG_HP = 3

#: LOBPCG_BLOCKS  how many blocks a state holds, the four above
LOBPCG_BLOCKS = 4

#: LOBPCG_LOCK_FLOOR  a row whose residual is under max(tol, this) stops
#:     offering search directions, however small tol is
LOBPCG_LOCK_FLOOR = 1.0e-9

#| How many entries each row of a block has.
def row_length(a) {
    dims = shape(a)
    return dims[1]
}

#| The Euclidean norm of each row of `a`: one number per vector in a block.
def row_norms(a) {
    norms = zeros(len(a))
    for b to len(a) {
        v = slab(a, b)
        norms[b] = sqrt(sum(abs(v) * abs(v)))
    }
    return norms
}

#| Divide each row of `a` by its entry in `norms`, in place. A zero row stays
#| zero instead of becoming a row of NaN.
def divide_rows(a, norms) {
    for b to len(a) {
        if norms[b] > 1.0e-300 {
            setslab(a, b, slab(a, b) / norms[b])
        }
    }
}

#| The rows of `top` followed by the rows of `bottom`, as one block.
def stacked(top, bottom) {
    ntop = len(top)
    both = czeros(ntop + len(bottom), row_length(top))
    for b to ntop {
        setslab(both, b, slab(top, b))
    }
    for b to len(bottom) {
        setslab(both, ntop + b, slab(bottom, b))
    }
    return both
}

#| The Rayleigh quotient of each row, $\lambda_b = x_b^H H x_b$ for a unit row
#| $x_b$: the best estimate of an eigenvalue that the vector alone can give.
def rayleigh_quotients(x, hx) {
    lam = zeros(len(x))
    for b to len(x) {
        lam[b] = sum(re(conj(slab(x, b)) * slab(hx, b)))

#: LOBPCG_X, LOBPCG_HX, LOBPCG_P, LOBPCG_HP  where each block sits in a LOBPCG
#:     state, a [4, nb, n] complex array: the current approximations x, H x,
#:     the previous step p, and H p
LOBPCG_X = 0
LOBPCG_HX = 1
LOBPCG_P = 2
LOBPCG_HP = 3

#: LOBPCG_BLOCKS  how many blocks a state holds, the four above
LOBPCG_BLOCKS = 4

#: LOBPCG_LOCK_FLOOR  a row whose residual is under max(tol, this) stops
#:     offering search directions, however small tol is
LOBPCG_LOCK_FLOOR = 1.0e-9

#| How many entries each row of a block has.
def row_length(a) {
    dims = shape(a)
    return dims[1]
}

#| The Euclidean norm of each row of `a`: one number per vector in a block.
def row_norms(a) {
    norms = zeros(len(a))
    for b to len(a) {
        v = slab(a, b)
        norms[b] = sqrt(sum(abs(v) * abs(v)))
    }
    return norms
}

#| Divide each row of `a` by its entry in `norms`, in place. A zero row stays
#| zero instead of becoming a row of NaN.
def divide_rows(a, norms) {
    for b to len(a) {
        if norms[b] > 1.0e-300 {
            setslab(a, b, slab(a, b) / norms[b])
        }
    }
}

#| The rows of `top` followed by the rows of `bottom`, as one block.
def stacked(top, bottom) {
    ntop = len(top)
    both = czeros(ntop + len(bottom), row_length(top))
    for b to ntop {
        setslab(both, b, slab(top, b))
    }
    for b to len(bottom) {
        setslab(both, ntop + b, slab(bottom, b))
    }
    return both
}

#| The Rayleigh quotient of each row, $\lambda_b = x_b^H H x_b$ for a unit row
#| $x_b$: the best estimate of an eigenvalue that the vector alone can give.
def rayleigh_quotients(x, hx) {
    lam = zeros(len(x))
    for b to len(x) {
        lam[b] = sum(re(conj(slab(x, b)) * slab(hx, b)))