Notebooks · lib/hf/eigen

Three symmetric eigensolvers that check each other

Every SCF iteration ends in a symmetric eigenproblem: the Fock matrix, written in an orthonormal basis, is diagonalized, and its lowest eigenvectors become the occupied orbitals. This module holds three solvers that share no numerical kernel (only the final sort and sign fix):

eig_ql (member 0, production): Householder reduction to a tridiagonal matrix, then the implicit QL iteration with Wilkinson shifts. This is the tred2/tql2 pair of EISPACK in the form JAMA gives it.
eig_jacobi: cyclic Jacobi rotations, each one zeroing one off-diagonal pair, until every off-diagonal entry is too small to change the diagonal it sits next to. It is slower, but it works on the full matrix and never forms a tridiagonal, so a bug in the reduction cannot hide in it. The method is the cyclic Jacobi of G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed. (Johns Hopkins, 2013), Algorithms 8.5.1 and 8.5.3, with the relative stopping test of J. Demmel and K. Veselic, SIAM J. Matrix Anal. Appl. 13, 1204 (1992).
eig_builtin: the language's own eig_sym, a Jacobi code written outside this library (in Rust for the interpreter, in C for the compiled runtime), with its own rotation formula and stopping rule.

Each member takes (A, B, opts) with a symmetric A[n, n]; B and opts are not read (the stage signature is shared with the orthogonalizers). It returns [n + 1, n]: row 0 holds the eigenvalues in ascending order, and row 1 + i column k holds component i of eigenvector k.

Eigenvectors are only defined up to sign, and inside a degenerate level up to a rotation. Every member fixes the sign the same way (the component of largest magnitude is positive; the first one wins a tie), so two members agree vector by vector wherever the eigenvalues are separated. Inside a degenerate level only the projector onto the whole level is comparable, and that is what the SCF uses: the density is the projector onto the occupied levels.

eig_of(A, k) runs member k on A; eig_values and eig_vectors take the eigenvalues and the eigenvector matrix out of its result.

Used with import hf.eigen.

Look inside this program

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

This program declares a record type (EigAgreement), and records run only in compiled C for now, so the Run buttons here are off. Open it in the workbench and use its C tab to compile it.

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 44–662

# the members of eig_methods (ScfOptions.eigensolver)
EIG_QL = 0
EIG_JACOBI = 1
EIG_BUILTIN = 2
EIG_EPS = 2.220446049250313e-16
# eig_tql2: the most QL iterations it spends on one eigenvalue
EIG_TQL2_MAXIT = 60
# cyclic Jacobi (eig_jacobi_raw): the most sweeps it runs
EIG_JACOBI_MAXSWEEP = 100
# Park-Miller's minimal standard generator, for the test matrices
EIG_MINSTD_A = 16807.0
EIG_MINSTD_M = 2147483647.0

#| How far the eigensolvers are apart on one matrix (eig_agreement): the
#| largest eigenvalue gap and occupied-projector gap to member 0, and the
#| worst residual `|A V - V diag(w)| / |A|` and orthogonality defect
#| `|V^T V - I|` of any member.
record EigAgreement {
    values: Num
    projector: Num
    residual: Num
    orth: Num
}

# ---- shared helpers --------------------------------------------------------------

#| Sort the eigenvalues in `d` ascending, moving the columns of `V` with them.
#| Selection sort: n is at most a few dozen, and a stable, simple order is
#| what makes the members comparable.
def eig_sort(d: Arr, V: Arr, n) {
    for i to n - 1 {
        k = i
        p = d[i]
        for j to n {
            if j > i and d[j] < p {
                k = j
                p = d[j]
            }
        }
        if k != i {
            d[k] = d[i]
            d[i] = p
            t = V[:, i]
            V[:, i] = V[:, k]
            V[:, k] = t
        }
    }
    return 0
}

#| Flip each column of `V` so its largest-magnitude component is positive.
def eig_fix_signs(V: Arr, n) {
    for k to n {
        big = 0.0
        at = 0
        for i to n {
            if abs(V[i, k]) > big {
                big = abs(V[i, k])
                at = i
            }

# the members of eig_methods (ScfOptions.eigensolver)
EIG_QL = 0
EIG_JACOBI = 1
EIG_BUILTIN = 2
EIG_EPS = 2.220446049250313e-16
# eig_tql2: the most QL iterations it spends on one eigenvalue
EIG_TQL2_MAXIT = 60
# cyclic Jacobi (eig_jacobi_raw): the most sweeps it runs
EIG_JACOBI_MAXSWEEP = 100
# Park-Miller's minimal standard generator, for the test matrices
EIG_MINSTD_A = 16807.0
EIG_MINSTD_M = 2147483647.0

#| How far the eigensolvers are apart on one matrix (eig_agreement): the
#| largest eigenvalue gap and occupied-projector gap to member 0, and the
#| worst residual `|A V - V diag(w)| / |A|` and orthogonality defect
#| `|V^T V - I|` of any member.
record EigAgreement {
    values: Num
    projector: Num
    residual: Num
    orth: Num
}

# ---- shared helpers --------------------------------------------------------------

#| Sort the eigenvalues in `d` ascending, moving the columns of `V` with them.
#| Selection sort: n is at most a few dozen, and a stable, simple order is
#| what makes the members comparable.
def eig_sort(d: Arr, V: Arr, n) {
    for i to n - 1 {
        k = i
        p = d[i]
        for j to n {
            if j > i and d[j] < p {
                k = j
                p = d[j]
            }
        }
        if k != i {
            d[k] = d[i]
            d[i] = p
            t = V[:, i]
            V[:, i] = V[:, k]
            V[:, k] = t
        }
    }
    return 0
}

#| Flip each column of `V` so its largest-magnitude component is positive.
def eig_fix_signs(V: Arr, n) {
    for k to n {
        big = 0.0
        at = 0
        for i to n {
            if abs(V[i, k]) > big {
                big = abs(V[i, k])
                at = i
            }