contrib/Eigen/StructuredMatrices)Fast operators and direct solvers for matrices determined by far fewer parameters than they have entries: circulant, Toeplitz, Hankel, block-circulant, Kronecker products and sums, diagonal-plus-low-rank, Vandermonde and Cauchy matrices, and the diagonal-plus-rank-one symmetric eigenproblem. Each operator stores only its generators or factors, applies itself without forming the matrix, and solves, factors or inverts itself through closed forms or algorithms that exploit the structure, for example in $O(n \log n)$ for circulant systems where a dense LU would take $O(n^3)$.
#include <contrib/Eigen/StructuredMatrices> // or the legacy <unsupported/Eigen/StructuredMatrices>
The module needs the FFT module (included by the umbrella header) and nothing outside Eigen. The class reference is the Doxygen block of each class; this file gives the overview, the mathematics, and the costs.
| Class | Matrix | Stored as | Product | Direct solve | Also |
|---|---|---|---|---|---|
Circulant | $C_{ij} = c_{(i-j) \bmod n}$ | first column | $O(n \log n)$ | solve() (pseudo-inverse), $O(n \log n)$ | eigen, SVD, rank, inverse, det |
Toeplitz | $T_{ij} = t_{i-j}$ | first column, first row | $O(p \log p)$ | LookAheadLevinson, $O(n^2)$ | solver conditionEstimate() |
Hankel | $H_{ij} = h_{i+j}$ | $h$, length $m+n-1$ | $O(p \log p)$ | solve(), $O(n^2)$ | toToeplitz() |
Bccb | 2-D circular convolution | $n_2 \times n_1$ array | $O(N \log N)$ | solve() (pseudo-inverse), $O(N \log N)$ | eigen, SVD, rank, inverse, det |
KroneckerOperator | $A \otimes B$ | $A$, $B$ | $O(N (n_1 + n_2))$ | solve(), leastSquaresSolve() | eigen, SVD, rank, inverse, det |
KroneckerSum | $A \oplus B$ | $A$, $B$ | $O(N (n_1 + n_2))$ | BartelsStewart, $O(N \sum_k n_k)$ | eigen |
DiagonalPlusLowRank | $D + U V^H$ | $d$, $U$, $V$ | $O(nk)$ | solve() (Woodbury), $O(nk^2 + k^3)$ | inverse, det |
Vandermonde | $V_{ij} = x_i^{\,j}$ | nodes | $O(mn)$ | BjorckPereyra, $O(n^2)$ | det |
Cauchy | $C_{ij} = 1/(x_i - y_j)$ | nodes | $O(mn)$ | CauchyLU, $O(n^2)$ | det |
DPR1EigenSolver | $D + \rho z z^T$ | — | — | — | eigen, $O(n^2)$ |
Sizes: $n$ for a square operator, $m \times n$ for a rectangular one, $p$ the 5-smooth embedding length $\ge m+n-1$, $N = n_1 n_2$ for the two-level operators, $k$ the rank of a low-rank correction. Product costs are per right-hand side; the Kronecker entries are for dense square factors.
Every class in the table but DPR1EigenSolver is an operator: it derives from EigenBase, stores its own copy of its generators, and provides rows(), cols() and coeff(i, j). It is built by its constructor or by a make* factory that deduces the compile-time sizes from the arguments (makeCirculant, makeToeplitz, makeHankel, makeBccb, makeKroneckerOperator, makeKroneckerSum, makeDiagonalPlusLowRank, makeVandermonde, makeCauchy).
Products. A * x returns an Eigen Product expression for a vector or a multi-column x. Assigning it behaves like any dense product: a temporary resolves aliasing between the destination and x, and .noalias() skips it. A real operator applied to a complex right-hand side yields a complex result.
Materialization. Assigning an operator to a dense matrix forms it explicitly (MatrixXd M = A;, and also M += A;, M -= A;). KroneckerOperator and KroneckerSum also materialize into a SparseMatrix (SparseMatrix<double> S; S = K;).
Transposition. transpose(), conjugate() and adjoint() return an operator of the same kind for every class except Vandermonde, whose transpose has no fast product. The FFT-based operators transform their cached symbol (index reversal, conjugation, or a phase factor) instead of recomputing an FFT.
Iterative solvers. Since A * x is an ordinary product expression, the operators run matrix-free in ConjugateGradient, BiCGSTAB, GMRES, MINRES, and, through adjoint(), the least-squares solvers LSMR and LeastSquaresConjugateGradient. They expose no coefficient storage, so the solver must be instantiated with IdentityPreconditioner:
Toeplitz<double> T(col, row); // rectangular, m > n LSMR<Toeplitz<double>, IdentityPreconditioner> lsmr(T); VectorXd x = lsmr.solve(b); // O(p log p) per iteration
LookAheadLevinson, BjorckPereyra, CauchyLU and BartelsStewart derive from SolverBase and follow the usual decomposition style: compute(op) or the constructor, a lazy solve(b) with any number of right-hand sides, info(), and transpose().solve(b) and adjoint().solve(b), which reuse the same factorization at the same cost. Circulant, Bccb, Hankel, KroneckerOperator, KroneckerSum and DiagonalPlusLowRank also have a solve(b) member that returns a plain matrix. Circulant and Bccb solve with their cached symbol; the others redo their setup on every call (a Levinson factorization, factor LUs, Schur forms, or the capacitance LU), so for repeated solves use LookAheadLevinson on toToeplitz() or BartelsStewart directly.
\mu \cdot 2^e$, renormalizing every factor and the running product to $\max(|\Re\mu|, |\Im\mu|) \in [\tfrac12, 1)$: the mantissa/exponent split of LINPACK's xGEDI [4], in base 2 rather than 10 so that the rescaling is exact [12]. A partial product therefore cannot overflow or underflow when the determinant itself is representable, whatever the ordering of the factors; zeros and non-finite factors propagate exactly.SVDBase: a singular value $\sigma$ counts as zero exactly when $\sigma < \max(c\,\varepsilon\,\sigma_{\max},\ \text{smallest normal})$ with $c = \min(\text{rows}, \text{cols})$. The pseudo-inverse solves of Circulant and Bccb invert exactly the symbol entries that rank() counts.FFT module's default backend, kissfft unless EIGEN_FFTW_DEFAULT, EIGEN_MKL_DEFAULT, EIGEN_POCKETFFT_DEFAULT or EIGEN_DUCCFFT_DEFAULT selects another. kissfft is fast only for 2-, 3- and 5-smooth lengths, so products pad the transform to the next length $2^a 3^b 5^c$; spectral quantities (eigenvalues, solves, rank, determinant) use the exact length. Products use the FFT only when a dimension exceeds 32 (for Bccb, when $N$ does), and never for a single-row or single-column Hankel. Each thread reuses one FFT engine and its plan tables (a fresh engine per call under EIGEN_AVOID_THREAD_LOCAL).Circulant<Scalar, Size = Dynamic> is the $n \times n$ matrix generated by its first column $c$:
C_{ij} = c_{(i-j) \bmod n}, \qquad C = F^H \operatorname{diag}(\hat c)\, F, \qquad \hat c_k = \sum_{j=0}^{n-1} c_j\, \omega^{jk}, \quad \omega = e^{-2\pi i/n},
with $F_{jk} = \omega^{jk}/\sqrt{n}$ the unitary DFT matrix. The eigenvalues are the DFT $\hat c$ of the generator, the operator's symbol, computed once at construction; the eigenvectors are the columns $f_k$ of $F^H$, with entries $e^{2\pi i jk/n}/\sqrt{n}$, the same for every circulant matrix. Everything follows from this factorization:
| Member | Result | Cost |
|---|---|---|
C * x | $\mathcal{F}^{-1}(\hat c \odot \mathcal{F} x)$ | $O(n \log n)$ |
solve(b) | $C^+ b = \mathcal{F}^{-1}(\hat c^{+} \odot \mathcal{F} b)$, $\hat c_k^{+} = 1/\hat c_k$ for the entries rank() counts, else 0 | $O(n \log n)$ |
eigenvalues(), symbol() | $\hat c$, in DFT order | precomputed for $n > 32$ |
eigenvectors() | $[f_0, \dots, f_{n-1}]$, dense | $O(n^2)$ |
singularValues() | $\lvert \hat c_k \rvert$, decreasing | $O(n \log n)$ |
matrixU(), matrixV() | columns $e^{i \arg \hat c_k} f_k$ and $f_k$, in the order of singularValues(), dense | $O(n^2)$ |
rank(), determinant() | count of retained $\hat c_k$; $\prod_k \hat c_k$ | $O(n \log n)$ |
inverse() | the Circulant with column $\mathcal{F}^{-1}(1/\hat c)$ | $O(n \log n)$ |
transpose(), conjugate(), adjoint() | Circulant, symbol $\hat c_{-k}$, $\overline{\hat c_{-k}}$, $\overline{\hat c_k}$ (indices mod $n$) | $O(n)$ |
For a real generator the symbol is conjugate-symmetric, $\hat c_{n-k} = \overline{\hat c_k}$, the two modes of each pair carry identical singular values, and determinant() returns the real part of the product. inverse() requires a nonsingular operator; solve() handles singular ones.
VectorXd c = VectorXd::Random(n); Circulant<double> C(c); // or makeCirculant(c) VectorXd y = C * x; VectorXd z = C.solve(y); // minimum-norm least-squares solution VectorXcd lambda = C.eigenvalues();
Toeplitz<Scalar, Rows = Dynamic, Cols = Dynamic> is the $m \times n$ matrix with constant diagonals, generated by its first column $c$ and first row $r$ ($r_0$ is ignored; the diagonal comes from $c$):
T_{ij} = \begin{cases} c_{i-j}, & i \ge j, \\ r_{j-i}, & i < j. \end{cases}
The product embeds $T$ in the circulant matrix of 5-smooth size $p \ge m+n-1$ with first column $[c_0, \dots, c_{m-1}, 0, \dots, 0, r_{n-1}, \dots, r_1]^T$; $Tx$ is the leading $m$ entries of that circulant applied to $[x; 0]$, at $O(p \log p)$ with the embedding symbol computed once at construction. transpose() swaps the generators; it and conjugate() and adjoint() reuse the embedding symbol. column() and row() return the generators.
LookAheadLevinson<Scalar> solves square systems $Tx = b$ in $O(n^2)$ by the look-ahead Levinson algorithm of Chan and Hansen [1, 2]. The classical Levinson recursion breaks down at a singular or ill-conditioned leading principal submatrix, which general (indefinite, nonsymmetric) Toeplitz matrices often have; the look-ahead variant steps over up to $p_{\max} - 1$ consecutive ill-conditioned leading submatrices with a block step, and remains weakly stable when there are no longer runs.
| Member | Meaning |
|---|---|
setMaxBlockSize(p) | $p_{\max}$, default 4; call before compute() |
info() | Success, or NumericalIssue when conditionEstimate() reaches $1/\varepsilon$; a run of ill-conditioned leading submatrices longer than the look-ahead is reported only through that estimate |
conditionEstimate() | estimate of $\kappa_2(T)$, the “algorithm condition number” of [1] |
transpose().solve(b), adjoint().solve(b) | reuse the factorization through the persymmetry $T^T = E T E$ ($E$ the exchange matrix) |
Toeplitz<double> T(col, row); // square LookAheadLevinson<double> levinson(T); VectorXd x = levinson.solve(b); // O(n^2) if (levinson.info() != Success) { /* condition estimate reached 1/eps */ }
Hankel<Scalar, Rows = Dynamic, Cols = Dynamic> is the $m \times n$ matrix with constant anti-diagonals, $H_{ij} = h_{i+j}$. It is built from its first column and last row (the shared corner entry comes from the column) and stores the single sequence $h$ of length $m+n-1$ (generator(), column(), lastRow()). A real square Hankel matrix is symmetric.
A Hankel matrix is a column-reversed Toeplitz matrix, $H = TE$ with $T_{ij} = h_{i+n-1-j}$ (toToeplitz()), so its product is a convolution evaluated in $O(p \log p)$ like the Toeplitz one. Single-row and single-column operators skip the FFT. transpose() keeps $h$ and swaps the dimensions; its symbol is the cached one multiplied by the phases $e^{2\pi i f (m-n)/p}$, with no new FFT.
solve(b) solves a square system in $O(n^2)$ as $T (Ex) = b$ with LookAheadLevinson. It factorizes on every call and eigen_asserts that the recursion did not break down; to reuse the factorization or check info() instead, run LookAheadLevinson on toToeplitz() and reverse the rows of its solution.
Bccb<Scalar, BlockSize = Dynamic, NumBlocks = Dynamic> is a block circulant matrix with circulant blocks, the matrix of a two-dimensional circular convolution. It is $N \times N$, $N = n_1 n_2$: an $n_1 \times n_1$ block circulant whose blocks are $n_2 \times n_2$ circulants. It is generated by the $n_2 \times n_1$ array $G$ whose column $k$ is the first column of block $k$; with $i = i_1 n_2 + i_2$ and $j = j_1 n_2 + j_2$,
C_{ij} = G_{(i_2 - j_2) \bmod n_2,\ (i_1 - j_1) \bmod n_1}, \qquad C \operatorname{vec}(X) = \operatorname{vec}(G \circledast X),
where $\operatorname{vec}$ stacks columns and $\circledast$ is 2-D circular convolution. BCCB matrices are diagonalized by the 2-D DFT $F_{n_1} \otimes F_{n_2}$ [5, 6]: the symbol $\hat G$, the 2-D DFT of $G$, holds the eigenvalues, and eigenvalue $f_1 n_2 + f_2$ is $\hat G_{f_2 f_1}$. The API mirrors Circulant (solve, eigenvalues, eigenvectors, singularValues, matrixU, matrixV, rank, inverse, determinant, and the transposition family) at $O(N \log N)$, with $O(N^2)$ for the dense vector matrices. blockSize() and numBlocks() return $n_2$ and $n_1$.
BCCB matrices model blurring with periodic boundary conditions in image deblurring, and are the standard preconditioners for two-level Toeplitz (BTTB) systems [6].
KroneckerOperator<Lhs, Rhs> is the Kronecker product $A \otimes B$, $A$ of size $m_1 \times n_1$ and $B$ of size $m_2 \times n_2$, as an operator that is never formed. With $\operatorname{vec}$ stacking the columns of an $n_2 \times n_1$ matrix $X$ (whatever the storage order) and $\operatorname{mat}$ its inverse,
(A \otimes B) \operatorname{vec}(X) = \operatorname{vec}(B X A^T),
at $O(m_2 n_1 (n_2 + m_1))$ per right-hand side instead of $O(m_1 m_2 n_1 n_2)$; right-hand sides are processed in cache-sized batches. Everything else factors through $A$ and $B$:
| Member | Identity | Notes |
|---|---|---|
solve(b) | $X = B^{-1} \operatorname{mat}(b) A^{-T}$ | one LU per dense factor, SparseLU per sparse factor; factors square and invertible |
leastSquaresSolve(b) | $(A \otimes B)^+ = A^+ \otimes B^+$ | one CompleteOrthogonalDecomposition per factor, rank decided per factor; rectangular or rank-deficient factors |
rank() | number of $\sigma_i(A)\,\sigma_j(B) \ge \min(\text{rows}, \text{cols})\,\varepsilon\,\sigma_{\max}(A)\,\sigma_{\max}(B)$ | two SVDs; can be below $\operatorname{rank}(A)\operatorname{rank}(B)$ |
determinant() | $\det(A)^{n_2} \det(B)^{n_1}$ | balanced accumulation of the factor LU pivots |
inverse() | $A^{-1} \otimes B^{-1}$ | a KroneckerOperator |
eigenvalues(), eigenvectors() | $\lambda_i(A)\,\mu_j(B)$ at index $i n_2 + j$; $V_A \otimes V_B$ | unsorted; the vectors are a KroneckerOperator |
singularValues(), matrixU(), matrixV() | $\sigma_i(A)\,\sigma_j(B)$; $U_A \otimes U_B$; $V_A \otimes V_B$ | thin SVD, unsorted; the vectors are KroneckerOperators |
transpose(), conjugate(), adjoint() | $A^T \otimes B^T$, $\bar A \otimes \bar B$, $A^H \otimes B^H$ |
The eigenvalues and singular values stay in Kronecker order because sorting them would destroy the Kronecker structure of the vector matrices.
Each factor can be any of:
| Factor passed | Stored as | Its side of a product | Its solve |
|---|---|---|---|
| dense expression | Matrix | GEMM | PartialPivLU |
DiagonalMatrix or .asDiagonal() expression | its diagonal, $O(n)$ | scaling | entrywise division |
SparseMatrix | compressed SparseMatrix | sparse-dense product, $O(\operatorname{nnz})$ per column | SparseLU, factorized once |
MatrixXd::Identity(m, n) | dimensions only | skipped when square; a rectangular one keeps the leading $\min(m, n)$ rows or columns | identity |
KroneckerOperator | as is (nested) | its own vec identity | its own factor solves |
KroneckerSum | as is (nested) | $B X + X A^T$ | BartelsStewart |
makeKroneckerOperator(a, b, c, ...) nests to the right, so products of three or more factors stay implicit. The decomposition family (eigenvalues through matrixV, rank, leastSquaresSolve) works on dense copies of diagonal, identity and sparse factors; rank and leastSquaresSolve also materialize a nested factor. A sparse factor whose SparseLU fails solves and inverts to NaN.
SparseMatrix<double> A = ...; // n x n auto K = makeKroneckerOperator(MatrixXd::Identity(p, p), A, MatrixXd::Identity(q, q)); // I_p (x) A (x) I_q VectorXd y = K * x; // O(p q nnz(A)), nothing of size p n q formed
Unlike kroneckerProduct() from the KroneckerProduct module, an expression meant to be evaluated into a matrix, KroneckerOperator is meant to be applied and solved with as it is.
KroneckerSum<Lhs, Rhs> is the Kronecker sum of square factors, $A$ of size $n_1$ and $B$ of size $n_2$,
A \oplus B = A \otimes I_{n_2} + I_{n_1} \otimes B, \qquad (A \oplus B) \operatorname{vec}(X) = \operatorname{vec}(B X + X A^T),
the operator of separable discretizations on tensor-product grids: with $D_x, D_y$ the 1-D second-difference matrices and $x$ the fast index of the unknown vector, the 2-D Laplacian is $D_y \oplus D_x$ and the 3-D one $D_z \oplus D_y \oplus D_x$. The factors may be of any kind KroneckerOperator accepts, including a KroneckerSum (makeKroneckerSum(a, b, c, ...) nests to the right). The product costs one product with each factor, $O(N (n_1 + n_2))$ for dense factors and $O(n_1 \operatorname{nnz}(B) + n_2 \operatorname{nnz}(A))$ for sparse ones, with no identity ever formed. eigenvalues() returns the sums $\lambda_i(A) + \mu_j(B)$ at index $i n_2 + j$ and eigenvectors() the KroneckerOperator $V_A \otimes V_B$.
BartelsStewart<KroneckerSumType> solves $(A_1 \oplus \cdots \oplus A_d)\, x = b$ for any number of factors. With Schur forms $A_k = Q_k T_k Q_k^H$,
A_1 \oplus \cdots \oplus A_d = Q\, (T_1 \oplus \cdots \oplus T_d)\, Q^H, \qquad Q = Q_1 \otimes \cdots \otimes Q_d,
so a solve applies $Q^H$, back-substitutes the (quasi-)upper-triangular middle factor, and applies $Q$; for two factors this is the Bartels-Stewart algorithm [3] for the Sylvester equation $B X + X A^T = \operatorname{mat}(b)$. The Schur forms depend on the factors:
isHermitian()): eigendecompositions, so the middle factor is diagonal and the solve divides by the eigenvalue sums, the fast diagonalization method [7], in real arithmetic for real factors;Q_k$ orthogonal and $T_k$ quasi-upper-triangular with a $2 \times 2$ block per complex conjugate eigenvalue pair; the coupled block rows are solved jointly, ending in dense solves at most $2^d$ wide, as LAPACK's xTRSYL does for two factors;T_k$ upper triangular.Setup costs one $O(n_k^3)$ decomposition per factor (sparse factors are densified for it) and each solve $O(N \sum_k n_k)$ per right-hand side, plus up to $O(4^d N)$ for the coupled solves of the real path. KroneckerSum::solve(b) builds a BartelsStewart for that one call.
The system is singular exactly when some sum $\lambda_{i_1}(A_1) + \cdots + \lambda_{i_d}(A_d)$ vanishes, which, as with PartialPivLU, nothing detects. The residual satisfies $\lVert b - Mx \rVert = O\big((\textstyle\sum_k n_k)\,\varepsilon\,(\sum_k \lVert A_k \rVert)\,\lVert x \rVert\big)$, relative to $\sum_k \lVert A_k \rVert$ rather than $\lVert M \rVert$, which matters when a shift is split across the factors as in $(A + cI) \oplus (B - cI)$. info() reports InvalidInput for a non-finite factor and NoConvergence when a Schur or eigenvalue iteration fails; the solve then returns NaN.
// Implicit Euler for u' = (Dy (+) Dx) u: (I - tau L) u_new = u, with I (x) I = I absorbed into one factor. SparseMatrix<double> Iy(ny, ny); Iy.setIdentity(); auto M = makeKroneckerSum(Iy - tau * Dy, -tau * Dx); BartelsStewart<decltype(M)> step(M); // Schur forms, once for (int k = 0; k < steps; ++k) u = step.solve(u);
DiagonalPlusLowRank<Scalar, Size = Dynamic, Rank = Dynamic> is $A = D + U V^H$ with $D = \operatorname{diag}(d)$ and $U, V$ of size $n \times k$ (diagonal(), factorU(), factorV(), correctionRank(); $k = 0$ is allowed). Products cost $O(nk)$. The Woodbury identity [8] and the matrix determinant lemma [9] reduce everything else to the $k \times k$ capacitance matrix $S = I_k + V^H D^{-1} U$ (capacitance()):
A^{-1} = D^{-1} - D^{-1} U S^{-1} V^H D^{-1}, \qquad \det A = \det(D) \det(S).
solve(b) costs $O(nk^2 + k^3)$ to set up and $O(nk)$ per right-hand side; determinant() is $O(nk^2 + k^3)$; inverse() returns the DiagonalPlusLowRank $D^{-1} + \tilde U \tilde V^H$ with $\tilde U = -D^{-1} U S^{-1}$ and $\tilde V = D^{-H} V$; transpose() is $D + \bar V \bar U^H$ and adjoint() is $\bar D + V U^H$.
solve and inverse need $D$ and $S$ nonsingular and determinant needs $D$ nonsingular, all with finite reciprocals $1/d_i$; this is not checked beyond NaN/Inf propagation. The splitting routes the solution through $D^{-1} b$, so accuracy degrades when $\max_i |1/d_i|$ greatly exceeds $\lVert A^{-1} \rVert$, even for a well-conditioned $A$ [10].
Vandermonde<Scalar, Rows = Dynamic, Cols = Dynamic> is the $m \times n$ matrix $V_{ij} = x_i^{\,j}$ of the nodes $x$ (nodes()); Vandermonde(x) is square and Vandermonde(x, n) has $n$ columns. V * a evaluates the polynomial with ascending coefficients $a$ at every node by Horner's rule,
(Va)_i = \sum_{j=0}^{n-1} a_j x_i^{\,j}: \qquad p_{n-1} = a_{n-1}, \quad p_j = a_j + x_i\, p_{j+1},
at $O(mn)$, the cost of a dense product but with $O(m)$ storage. determinant() uses $\det V = \prod_{i<j} (x_j - x_i)$ in $O(n^2)$. The transpose has no fast product, so Vandermonde is not closed under transposition; least-squares problems are best solved by a dense QR of the materialized matrix. Integer scalar types are rejected.
BjorckPereyra<Scalar> solves square systems in $O(n^2)$ operations and $O(n)$ storage [11, §4.6]: $Va = f$ is polynomial interpolation, solved by divided differences in the Newton basis followed by the change to the monomial basis, and the dual (moment) system $V^T w = b$ runs the transposed recurrences through transpose().solve(b). There is no factorization: compute() stores the nodes and info() reports NumericalIssue for an exactly repeated node or InvalidInput for a non-finite one. Genuinely complex nodes are put in Leja order [13] to limit growth in the Newton representation; real nodes keep their input order.
Real-node Vandermonde matrices are exponentially ill-conditioned: $\kappa_2(V)$ grows at least like $2^n$ for every real node set [14]. Björck-Pereyra solves are nevertheless often far more accurate than that suggests: for monotone nodes and a sign-alternating right-hand side the forward error obeys a small relative bound independent of $\kappa(V)$ [15, ch. 22]. Nodes on the unit circle are the well-conditioned case; for the $n$-th roots of unity $V/\sqrt{n}$ is unitary.
Cauchy<Scalar, Rows = Dynamic, Cols = Dynamic> is $C_{ij} = 1/(x_i - y_j)$ for row nodes $x$ and column nodes $y$ (rowNodes(), colNodes()), stored as the $m + n$ nodes. Every $x_i$ must differ from every $y_j$; this is not checked. Products are evaluated directly at $O(mn)$ with $O(m)$ extra storage. The class is closed under transposition, $C(x, y)^T = C(-y, -x)$ and $C(x, y)^H = C(-\bar y, -\bar x)$, and the determinant of a square Cauchy matrix has the closed form
\det C = \frac{\prod_{i<j} (x_j - x_i)(y_i - y_j)}{\prod_{i,j} (x_i - y_j)},
evaluated in $O(n^2)$. The Hilbert matrix is the Cauchy matrix with $x_i = i + 1$, $y_j = -j$.
CauchyLU<Scalar> factors $PC = LU$ with partial pivoting in $O(n^2)$ by the Gohberg-Kailath-Olshevsky algorithm [16]. A Cauchy matrix satisfies the rank-one displacement equation $D_x C - C D_y = \mathbf{1}\mathbf{1}^T$, which survives row permutations, so every Schur complement is again Cauchy-like, $S_{ij} = a_i b_j / (x_i - y_j)$, and elimination step $k$ only updates the generators:
a_i \leftarrow a_i \frac{x_i - x_k}{x_i - y_k}, \qquad b_j \leftarrow b_j \frac{y_j - y_k}{y_j - x_k}.
Each column of the Schur complement is generated from $O(n)$ data, so the pivot search costs nothing extra: this is genuine partial pivoting at fast-algorithm cost, which the Levinson family for Toeplitz matrices lacks. The factors are stored densely ($O(n^2)$ memory); info() reports NumericalIssue when a zero or non-finite entry prevents the factorization.
DPR1EigenSolver<RealScalar> (float, double or long double) computes all eigenvalues and eigenvectors of the real symmetric diagonal-plus-rank-one matrix $A = D + \rho z z^T$ in $O(n^2)$. It is the kernel at the heart of divide-and-conquer symmetric eigensolvers (LAPACK's xLAED2/xLAED3/xLAED4), as a standalone class. After deflation (entries with negligible $|z_i|$ are already eigenpairs, and nearly equal $d_i$ are merged by Givens rotations whose dropped coupling is below a backward-stability threshold), the remaining eigenvalues are the roots of the secular equation
f(\lambda) = 1 + \rho \sum_i \frac{z_i^2}{d_i - \lambda} = 0,
one between each pair of consecutive poles and, for $\rho > 0$, one beyond the largest pole. Each root is bracketed and bisected in coordinates shifted to its nearest pole, so every $\lambda - d_i$ is an exact data difference plus a small offset. The eigenvectors are built from the Gu-Eisenstat vector [17]
\hat z_i^2 = \frac{\prod_j (\lambda_j - d_i)}{\rho \prod_{j \ne i} (d_j - d_i)}, \qquad v_j \propto (D - \lambda_j I)^{-1} \hat z,
for which the computed $\lambda_j$ are the exact eigenvalues of $D + \rho \hat z \hat z^T$; this makes the computed $V$ numerically orthogonal without reorthogonalization. Either sign of $\rho$, $\rho = 0$, zero $z$, repeated $d_i$ and any ordering of $d$ are supported; the problem is scaled internally by an exact power of two.
DPR1EigenSolver<double> es(d, rho, z); // or es.compute(d, rho, z, EigenvaluesOnly) VectorXd lambda = es.eigenvalues(); // ascending MatrixXd V = es.eigenvectors(); // orthogonal
info() reports InvalidInput for non-finite input, an overflowing $\rho \lVert z \rVert^2$ or a non-finite eigenvalue, and NoConvergence when a secular root could not be resolved.
The tests are contrib/test/structured_*.cpp: structured_matrices (Circulant, Toeplitz, LookAheadLevinson, Hankel), structured_bccb, structured_kronecker, structured_kronecker_sum, structured_dplr, structured_vandermonde (and structured_vandermonde_int_index), structured_cauchy and structured_dpr1.
cmake -G Ninja -S . -B build cmake --build build --target structured_matrices structured_bccb structured_kronecker structured_kronecker_sum \ structured_dplr structured_vandermonde structured_vandermonde_int_index structured_cauchy structured_dpr1 ctest --test-dir build -R '^structured_' --output-on-failure --no-tests=error
Benchmarks are in contrib/benchmarks/StructuredMatrices (bench_structured_*, one per operator and for the solves, transposition, materialization and the Kronecker spectral paths) and benchmarks/StructuredMatrices (FFT operators, Vandermonde and Cauchy, diagonal-plus-low-rank). Each tree is a standalone CMake project:
cmake -G Ninja -S contrib/benchmarks -B build-contrib-bench -DCMAKE_BUILD_TYPE=Release cmake --build build-contrib-bench --target bench_structured_circulant