xSTEMR: are zero eigenvalues sufficiently accurate?

[christoph-conrads]

Continuing from #997, the SciPy maintainers noted that xSTEMR computes zero eigenvalues less accurately than expected for a full real symmetric 3x3 matrix with eigenvalues 3, 1, and 0. I thought that this is just a side effect of the small size but testing with random matrices shows that xSTEMR can compute an approximation as large as 20 ε max(|λ|) to the smallest eigenvalue for a 6x6 matrix with eigenvalues 3, 3, 3, 3, 1, 0 (code below). Personally, I am not sure I would consider such an eigenvalue zero anymore because it is larger than ε n max(|λ|), where n is the dimension of the matrix.

Is this approximation of a zero eigenvalue considered a bug?

Democode (requires Python3, NumPy, and Python):

#!/usr/bin/python3

# The code below computes with various methods ("drivers" in SciPy lingo) the
# eigenvalues of a random matrix with fixed eigenvalues and computes the
# accuracy of one of the eigenvalues. A five-number summary is shown for each
# method.

import sys

import numpy as np
import numpy.linalg as la
import scipy
from scipy.linalg import pinvh, pinv


def sample_smallest_ev(rng, num_samples: int, driver: str, dtype=np.float64):
    assert num_samples > 5

    eigs = np.array([3, 3, 3, 3, 1, 0], dtype=dtype)
    errors = []

    for i in range(num_samples):
        n = len(eigs)
        a = rng.uniform(-1, 1, size=(n, n)).astype(dtype)
        q, r = np.linalg.qr(a)
        d = np.diag(r)
        q = q @ np.diag(d / abs(d))

        eps = np.finfo(dtype).eps
        assert la.norm(q.T @ q - np.eye(n)) <= 4 * n * eps * la.norm(q)

        x = q @ np.diag(eigs) @ q.T
        x = np.triu(x) + np.triu(x, 1).T

        assert np.all(x == x.T)

        if driver == "srf":
            es = scipy.linalg.eigvalsh(x)
        else:
            es = scipy.linalg.eigh(x, driver=driver)[0]

        f = min
        errors.append(f(es) - f(eigs))

        if driver == "srf":
            continue

        y1 = pinvh(x)
        y2 = pinv(x)

        np.allclose(np.matmul(np.matmul(x, y1), x), x)  # true
        np.allclose(np.matmul(np.matmul(y1, x), y1), y1)  # false
        np.allclose(np.conjugate(np.matmul(x, y1)), np.matmul(x, y1))  # true
        np.allclose(np.conjugate(np.matmul(y1, x)), np.matmul(y1, x))  #true

        np.allclose(np.matmul(np.matmul(x, y2), x), x)  # true
        np.allclose(np.matmul(np.matmul(y2, x), y2), y2)  # true
        np.allclose(np.conjugate(np.matmul(x, y2)), np.matmul(x, y2))  # true
        np.allclose(np.conjugate(np.matmul(y2, x)), np.matmul(y2, x))  #true

    errors = sorted(map(lambda x: abs(x) / (max(eigs) * eps), errors))
    percentiles = range(0, 101, 25)
    five_number_summary = list(map(lambda p: np.percentile(errors, p),
                                   percentiles))

    return five_number_summary


def main():
    num_samples = 1000
    # srf = square-root free / Pal-Walker-Kahan QR algorithm (eigenvalues only)
    for driver in ["evd", "evr", "srf"]:
        rng = np.random.RandomState(seed=1)
        summary = sample_smallest_ev(rng, num_samples, driver=driver)

        print("n=%d driver=%-3s %+8.2e %+8.2e +%8.2e %+8.2e %+8.2e" %
              (num_samples, driver, *summary))


if __name__ == "__main__":
    sys.exit(main())

Output on my x86-64 machine (SciPy 1.10.1, OpenBLAS 0.3.21):

n=1000 driver=evd +0.00e+00 +3.12e-01 +6.46e-01 +1.15e+00 +4.12e+00
n=1000 driver=evr +0.00e+00 +3.33e+00 +7.33e+00 +1.07e+01 +2.00e+01
n=1000 driver=srf +0.00e+00 +2.29e-01 +4.79e-01 +8.55e-01 +2.69e+00