Nightmare Matrices 2: Symmetric Positive Definite Case

June 01, 2026

I present here a monstrous counter-example for sparse linear solvers which operate only on an input sparse matrix. This counter-example is closely related to the nightmare matrices of a previous post, but with a key strengthening: the matrices I present here will be symmetric and positive definite (“SPD” for short). In my prior post I relied on indefiniteness of the matrix to make iterative solutions a non-viable “shortcut” around direct methods. It remained an open question to me whether the expander symbolic structure of a nightmare matrix could even admit an SPD formulation that remains difficult to solve.

The reason it is not obvious is because the classic archetype SPD matrix, the graph Laplacian tends to have ideal spectral distributions for Krylov methods. Independent of conditioning I have observed that expander matrices that are SPD tend to have eigenvalue distributions that are heavy in the spectral midpoint and not clustered at the endpoints. This is close to the best-case scenario for Krylov methods as they quickly isolate the midpoint from the small and large modes, and then rapidly converge from there. So the open question to me was could I force the worst-case spectral distribution while still preserving SPD and underlying expander structure? The answer turns out to be yes. The construction below gives matrices with the following properties:

  1. Sparse: \( O(1) \) nonzeros per row.
  2. Symmetric and positive definite
  3. \( O(n^2) \) fill in sparse direct factorization
  4. Can take any prescribed condition number with eigenvalues in [min,max]
  5. Eigenvalues can take Chebyshev nodal distribution in [min,max] (worst case for Krylov solvers like Conjugate Gradients).

Such a matrix will resist sparse direct solvers because of the excessive fill-in, it will resist many structure-based preconditioning efforts as the expander property ruins the underlying structure, and it will also resist attempts to solve it iteratively because the spectrum has been hand-picked to be the worst possible for Krylov iterative methods.

I will describe the construction, give numerical results confirming the above properties, and then I’ll give some mitigations a sparse solver library could potentially adopt to overcome my construction although I suspect the approach here could be strengthened to make this prohibitively expensive.

The Construction

The construction works by taking an initial SPD matrix with underlying expander structure and embedding it into a larger matrix that preserves the expander structure but has an enforced spectral distribution to make solving it with Krylov methods difficult.

Suppose \(A_0\in\mathbb{R}^{n\times n}\) is a symmetric positive definite matrix. Then define

$$ D=\operatorname{diag}(\delta_1,\dots,\delta_n), \qquad \delta_i>0. $$

And then take the orthogonal block matrix

$$ Q = \frac{1}{\sqrt{2}} \begin{bmatrix} I & I\\ I & -I \end{bmatrix}. $$

and set

$$ M = Q \begin{bmatrix} A_0 & 0\\ 0 & D \end{bmatrix} Q^T. $$

Equivalently,

$$ \boxed{ M = \frac12 \begin{bmatrix} A_0+D & A_0-D\\ A_0-D & A_0+D \end{bmatrix}. } $$

If \( A_0 \) is sparse, symmetric, and positive definite then so is \( M \). \(M\) is twice as large as \(A_0\) and has twice as many nonzeros per row, but by appropriately selecting \( D \) we can explicitly specify half of the spectrum. In fact, the spectrum of \(M\) is the union of the spectrum of \(A_0\) and the diagonal entries of \(D\).

In words this is reversing a block orthogonal decomposition where our starting point contains the planted eigenvalues as well as an expander matrix.

Sample Code

The following code can compute \(M \) given an \( A_0 \)

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla


def chebyshev_nodes(n, lo, hi, *, sort=True):
    """
    First-kind Chebyshev nodes mapped to [lo, hi].

    These are dense near both endpoints but do not include the endpoints.
    """
    j = np.arange(n, dtype=float)
    x = np.cos((2 * j + 1) * np.pi / (2 * n))
    vals = 0.5 * (lo + hi) + 0.5 * (hi - lo) * x
    return np.sort(vals) if sort else vals


def estimate_extreme_eigs_spd(A, *, dense_cutoff=512, tol=1e-8):
    """
    Estimate lambda_min and lambda_max of a real symmetric sparse SPD matrix.

    For small matrices, this uses dense eigvalsh for reliability.
    For larger matrices, this uses scipy.sparse.linalg.eigsh.
    """
    A = sp.csr_matrix(A)
    n, m = A.shape
    if n != m:
        raise ValueError("A must be square.")

    if n <= dense_cutoff:
        vals = np.linalg.eigvalsh(A.toarray())
        return float(vals[0]), float(vals[-1])

    lmin = spla.eigsh(
        A,
        k=1,
        which="SA",
        return_eigenvectors=False,
        tol=tol,
    )[0]

    lmax = spla.eigsh(
        A,
        k=1,
        which="LA",
        return_eigenvectors=False,
        tol=tol,
    )[0]

    return float(lmin), float(lmax)


def planted_spd_lift(
    A0,
    planted_eigs,
    *,
    fmt="csr",
    check_symmetric=True,
    symmetry_tol=1e-10,
):
    """
    Given an n x n sparse SPD matrix A0 and n positive target eigenvalues,
    build a 2n x 2n sparse SPD matrix M whose spectrum is

        spectrum(M) = spectrum(A0) union planted_eigs.

    The construction is

        M = 1/2 [[A0 + D, A0 - D],
                 [A0 - D, A0 + D]]

    where D = diag(planted_eigs).

    The planted eigenvectors are exactly

        [e_i; -e_i] / sqrt(2),

    with eigenvalue planted_eigs[i].

    The original eigenvectors v of A0 lift to

        [v; v] / sqrt(2).
    """
    A0 = sp.csr_matrix(A0)
    n, m = A0.shape

    if n != m:
        raise ValueError("A0 must be square.")

    if check_symmetric:
        diff = (A0 - A0.T).tocoo()
        sym_err = 0.0 if diff.nnz == 0 else float(np.max(np.abs(diff.data)))
        if sym_err > symmetry_tol:
            raise ValueError(
                f"A0 does not look symmetric; max |A0 - A0.T| = {sym_err:g}."
            )

    planted_eigs = np.asarray(planted_eigs, dtype=float)

    if planted_eigs.shape != (n,):
        raise ValueError(f"planted_eigs must have shape ({n},).")

    if np.any(planted_eigs <= 0):
        raise ValueError("All planted eigenvalues must be positive for SPD output.")

    D = sp.diags(planted_eigs, format="csr")

    M = 0.5 * sp.bmat(
        [
            [A0 + D, A0 - D],
            [A0 - D, A0 + D],
        ],
        format=fmt,
    )

    return M

The initial \( A_0 \) can be computed using the networkx library to make an expander:

import networkx as nx
graph = nx.random_regular_graph(degree, n, seed=seed)
edges = {(min(i, j), max(i, j)) for i, j in graph.edges()}

and a matrix can be made SPD by giving random positive off-diagonal nonzeros (initially setting the diagonal to zero) and then overwriting the diagonal as A@np.ones(m)*(1+eps).

Results

The matrix \(M\) constructed like above looks like the following, structurally:

matrix

you can visually get an idea for the expander structure by attempting a reordering like Reverse Cuthill-McKee (RCM) which attempts to move nonzeros into a band:

matrix

The RCM ordering doesn’t prove anything about the hardness of this system, but it is one useful datapoint and it is easy to visualize.

What I did next was the following, for both \( A_0 \) and \( M \) for different matrix sizes

  1. I solved a system with CHOLMOD and reported the best achieved fill reduction after trying all available methods
  2. I solved a system with Conjugate Gradients and reported number of iterations
  3. I computed the resulting bandwidth after RCM reordering
  4. I computed the “equivalent bandwidth” which is the bandwidth of a hypothetical matrix of the same size that would take the same memory as CHOLMOD’s computed \(L\) factor

I present a table for these below. The examples were chosen so that the resulting system would have \( \kappa(M) = 1e9 \)

The main things to look for are that CG takes essentially \(0.5n\) iterations on the lifted matrices, the fill grows rapidly with problem size, and RCM still leaves a very large bandwidth.

Matrix n nnz / row CG iters / n COLAMD fill COLAMD equiv bw / n Best fill Best equiv bw / n RCM bw / n
A0_for_M10245129.00.0644539.4560.1819.456 colamd0.1810.632812
M1024102418.00.5000009.4290.1819.429 colamd0.1810.636719
A0_for_M204810249.00.03320317.8740.17117.874 colamd0.1710.631836
M2048204818.00.50000017.8470.17117.847 colamd0.1710.632812
A0_for_M409620489.00.01660234.1720.16334.172 colamd0.1630.624512
M4096409618.00.50000034.1440.16334.144 colamd0.1630.625732
A0_for_M819240969.00.00830167.1390.16067.139 colamd0.1600.624023
M8192819218.00.50000067.0190.16067.019 colamd0.1600.624390
A0_for_M1638481929.00.004150132.8870.158132.887 colamd0.1580.623169
M163841638418.00.500000132.8590.158132.859 colamd0.1580.623474

The trend above appears to validate all of the properties I set out to satisfy: sparsity, hardness for CG, and O(n^2) fill-in. Not surprisingly the number of CG iterations it takes to solve this system ends up becoming \( \frac{1}{2}n \) and this is because once it has resolved the planted Chebyshev modes it only has to do comparatively little to get the remaining modes.

Potential Mitigations

Most solver libraries which only have access to the underlying sparse matrix and perhaps a few hints from the user (e.g. SPD, stencil structure, etc.) will not have preprocessing that can detect the block structure from the above matrix. But it is feasible that this structure could be detected, and the block orthogonal decomposition could be undone and the Chebyshev modes explicitly deflated away. We could achieve this by vertex matching algorithms since \( A_0 \) is replicated in a fairly straightforward block form into \( M \).

To defeat this, the \(D \) matrix could be made tridiagonal SPD with the desired spectrum instead of purely diagonal. This introduces more nonzeros into \(M \) and makes vertex matching less viable. More generally, the block form here is just the simplest version of the adversarial idea. I suspect it is possible to make discovering this block structure arbitrarily hard at the cost of introducing more nonzeros into \(D \).

Discussion

For expander graphs the graph Laplacian tends to be handled by Krylov methods very well. This is because the system itself is well conditioned, but also the eigenvalues are not clustered at the endpoints in a pathological way (i.e. like with Chebyshev nodes). The counter example here is designed to show that this is not a fundamental property of sparse expander SPD matrices.