Cauchy-like Displacement Structure

• Cauchy-like displacement structure is a representation principle for dense matrices that exploits a Sylvester-type equation to yield a low-rank residual from diagonal operator interactions.
• It enables fast arithmetic operations such as inversion, factorization, and matrix-vector multiplication by updating compact low-rank generators instead of processing full matrices.
• The methodology supports efficient linear system solvers and robust structured pivoting strategies, with applications in stability analysis and approximation through companion-based algorithms.

Cauchy-like displacement structure is a representation principle for dense matrices whose algebraic interaction with simple diagonal operators produces only a low-rank residual. In its classical form, a matrix $C$ is Cauchy-like when it satisfies a Sylvester-type equation

$D_t C - C D_s = GH^*,$

with diagonal displacement matrices $D_t=\operatorname{diag}(t_1,\dots,t_n)$, $D_s=\operatorname{diag}(s_1,\dots,s_n)$, and generators $G,H\in\mathbb{C}^{n\times r}$ of small column dimension $r$. The quantity $r$ is the displacement rank. This viewpoint treats the matrix not as an arbitrary $n\times n$ array, but as the solution of a low-rank operator equation, which is why Cauchy-like matrices are central in fast structured linear algebra (Chandrasekaran et al., 2018, Arico' et al., 2010).

1. Displacement-operator formulation

The general displacement paradigm starts from a linear operator applied to a matrix. One common form is

$$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$

and a matrix $T$ is said to have displacement structure with respect to $D_t C - C D_s = GH^*,$0 when $D_t C - C D_s = GH^*,$1 is small compared with $D_t C - C D_s = GH^*,$2. If

$D_t C - C D_s = GH^*,$3

then $D_t C - C D_s = GH^*,$4 are generators, and the matrix can be stored using $D_t C - C D_s = GH^*,$5 rather than $D_t C - C D_s = GH^*,$6 numbers. An analogous Sylvester formulation uses

$D_t C - C D_s = GH^*,$7

or, in the notation of fast solvers,

$D_t C - C D_s = GH^*,$8

In both settings, the essential object is the low-rank residual, not the density pattern of the matrix itself (Chandrasekaran et al., 2018, Bostan et al., 2017).

For the Cauchy-like case, the displacement matrices are diagonal: $D_t C - C D_s = GH^*,$9 The displacement rank is exactly the rank of the right-hand side, and the compact data structure is the generator pair $D_t=\operatorname{diag}(t_1,\dots,t_n)$0, which requires only $D_t=\operatorname{diag}(t_1,\dots,t_n)$1 data when $D_t=\operatorname{diag}(t_1,\dots,t_n)$2. This representation is particularly convenient because diagonal operators are trivial to apply, and because individual matrix entries can be reconstructed directly from the generators when the knot sets are separated (Arico' et al., 2010).

A structural identity that governs much of the algorithmics is

$D_t=\operatorname{diag}(t_1,\dots,t_n)$3

It explains why inversion often preserves low displacement rank, and hence why structured solvers can recurse on Schur complements or inverses without immediately leaving the class of efficiently representable matrices (Chandrasekaran et al., 2018).

2. Canonical Cauchy-like representations and neighboring classes

The pure Cauchy matrix is

$D_t=\operatorname{diag}(t_1,\dots,t_n)$4

and satisfies a rank-one displacement relation with diagonal operators. One form given in the literature is

$D_t=\operatorname{diag}(t_1,\dots,t_n)$5

Equivalently, the canonical Cauchy form is often written as

$D_t=\operatorname{diag}(t_1,\dots,t_n)$6

so that entrywise

$D_t=\operatorname{diag}(t_1,\dots,t_n)$7

Matrices of the latter form are called Cauchy-like: they generalize the pure Cauchy kernel by allowing multiplicative row and column generators while retaining very small displacement rank, typically $D_t=\operatorname{diag}(t_1,\dots,t_n)$8 in the pure case (Chandrasekaran et al., 2018).

In solver-oriented notation, a Cauchy-like matrix has entries

$D_t=\operatorname{diag}(t_1,\dots,t_n)$9

with $D_s=\operatorname{diag}(s_1,\dots,s_n)$0 and $D_s=\operatorname{diag}(s_1,\dots,s_n)$1. This is a generalized Cauchy matrix modulated by low-rank generators. If $D_s=\operatorname{diag}(s_1,\dots,s_n)$2 for all $D_s=\operatorname{diag}(s_1,\dots,s_n)$3, every entry is reconstructable from the displacement data; if some coincidence $D_s=\operatorname{diag}(s_1,\dots,s_n)$4 occurs, the corresponding entry cannot be reconstructed through the formula because the denominator vanishes, and the matrix is only partially reconstructable (Arico' et al., 2010).

Cauchy-like structure belongs to a larger family of displacement-structured classes. Toeplitz-like matrices arise from the down-shift operator through

$D_s=\operatorname{diag}(s_1,\dots,s_n)$5

with displacement rank at most $D_s=\operatorname{diag}(s_1,\dots,s_n)$6 for a Toeplitz matrix. Vandermonde matrices satisfy

$D_s=\operatorname{diag}(s_1,\dots,s_n)$7

Generalized Pick matrices are defined by the condition that

$D_s=\operatorname{diag}(s_1,\dots,s_n)$8

is small. In a unified operator language, Toeplitz-like, Vandermonde-like, Cauchy-like, and Pick matrices are all instances in which interaction with a simple operator pair leaves a low-rank residue (Chandrasekaran et al., 2018).

A broader algebraic generalization replaces diagonal or shift operators by block-diagonal companion matrices associated with pairwise-coprime polynomial families. In that framework, Cauchy-like matrices are simply the diagonal-diagonal specialization obtained when the relevant polynomials are linear factors. This places classical Cauchy-like structure inside a continuum of companion-based displacement operators that also covers Toeplitz, Hankel, and Vandermonde types in a unified manner (Bostan et al., 2017).

3. Fast arithmetic from generators

The main computational consequence of Cauchy-like displacement structure is that many matrix operations can be expressed in terms of generator updates rather than dense matrix manipulation. Sums behave particularly simply: $D_s=\operatorname{diag}(s_1,\dots,s_n)$9 Products are more delicate, but when one can apply the displacement matrices efficiently and handle $G,H\in\mathbb{C}^{n\times r}$0, product structure can still be updated at low cost. For inversion and linear solving, generalized Schur elimination proceeds recursively on Schur complements while preserving low displacement rank, which is why factorization can often be carried out in $G,H\in\mathbb{C}^{n\times r}$1 rather than dense $G,H\in\mathbb{C}^{n\times r}$2 time when the operator itself is cheap to manipulate (Chandrasekaran et al., 2018).

For matrix reconstruction and matrix-vector multiplication, the displacement equation can be inverted explicitly in favorable cases. One series representation is

$G,H\in\mathbb{C}^{n\times r}$3

valid when the spectral radii of $G,H\in\mathbb{C}^{n\times r}$4 and $G,H\in\mathbb{C}^{n\times r}$5 are small. For nilpotent shift matrices this reduces to Gohberg–Semencul-type formulas. In the Cauchy-like setting, the operators are diagonal, and this makes fast matrix-vector multiplication for the generalized Cauchy or Pick class especially natural because the action of the displacement matrices is entrywise and the entire matrix is encoded by node vectors and generators rather than by $G,H\in\mathbb{C}^{n\times r}$6 stored entries (Chandrasekaran et al., 2018).

The same principle underlies asymptotically faster algebraic algorithms. A central subproblem is computing $G,H\in\mathbb{C}^{n\times r}$7, where $G,H\in\mathbb{C}^{n\times r}$8 is structured with displacement rank $G,H\in\mathbb{C}^{n\times r}$9 and $r$0 is an arbitrary $r$1 matrix. Classical bounds range from $r$2 to $r$3, depending on the structure, whereas companion-operator methods reduce this to bounds of order $r$4 or $r$5. The resulting inversion and linear-system algorithms are Las Vegas algorithms obtained by combining fast structured multiplication with randomized regularization (Bostan et al., 2017).

4. Linear-system solvers: Schur-based and modular-polynomial methods

A practical solver for reconstructable Cauchy-like matrices is built around the generalized Schur algorithm, a structured analogue of Gaussian elimination that updates generators instead of full matrix entries. For a matrix $r$6 satisfying

$r$7

the first Schur complement is again Cauchy-like, with explicitly updated generators. This recursive generator update is the core mechanism by which elimination remains quadratic in $r$8 rather than cubic. The same framework extends to multiple right-hand sides by embedding the system $r$9 into the augmented matrix

$r$0

whose $r$1th Schur complement equals the solution: $r$2 The resulting solver requires $r$3 floating-point operations and $r$4 memory, and for $r$5 right-hand sides the cost remains proportional to $r$6 (Arico' et al., 2010).

Stability considerations lead to structured pivoting strategies. Partial pivoting preserves the Cauchy-like displacement structure well. Sweet–Brent pivoting approximates complete pivoting by comparing the largest candidate in the current row and column. Gu’s pivoting uses a compact QR factorization of the left generator to control growth. Complete pivoting is included for comparison, but it generally destroys the intended complexity by requiring reconstruction of large submatrices and raising the cost to $r$7. A crucial structural fact is that Cauchy-like matrices are invariant under row and column permutations, so these pivoting schemes can be incorporated without abandoning the structured framework (Arico' et al., 2010).

Repeated knots require additional care. If the knot vector has no repetitions, the augmented matrix is reconstructable except for known diagonal entries in the $r$8-block; if repetitions occur, some off-diagonal entries become nonreconstructable. Two remedies are described: redundant-injective splitting, which separates repeated and nonrepeated knots, and gathering repeated knots into contiguous blocks while storing the nonreconstructable blocks explicitly. The latter yields a general solver for multiple knots without changing the overall complexity (Arico' et al., 2010).

A different, fully deterministic line of work reformulates structured linear systems as modular equations on univariate polynomials. For a Cauchy-like matrix represented by

$r$9

with repetition-free and disjoint node sets $n\times n$0 and $n\times n$1, one uses

$n\times n$2

and the decomposition

$n\times n$3

The solve $n\times n$4 is rewritten as a polynomial congruence, then handled in three steps: vector M-Padé approximation over $n\times n$5, modular inversion over $n\times n$6, and simultaneous M-Padé approximation over $n\times n$7. For Toeplitz-like, Vandermonde-like, and Cauchy-like structures, this gives deterministic algorithms for linear systems and nullspace computation within $n\times n$8 operations, or more explicitly

$n\times n$9

for the Cauchy-like case (Khichane et al., 2 Mar 2026).

The asymptotically fastest companion-based algorithms are algebraic rather than numerical in emphasis. In particular, the 2017 framework is explicit that its algorithms assume exact arithmetic over a field and do not address numerical stability; its gains concern arithmetic complexity and the exploitation of low displacement rank, not floating-point backward stability (Bostan et al., 2017).

5. Singular values, separation, and numerical rank

Small displacement rank does not by itself imply small numerical rank. A standard caution is that the identity matrix is Toeplitz and the exchange matrix is Hankel, yet both are full rank. For matrices satisfying a Sylvester displacement equation

$$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$0

the singular values are controlled not merely by the number of columns of $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$1 and $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$2, but also by the relative location of the spectra of $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$3 and $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$4. When those spectral sets are well separated, singular values decay geometrically, and low-rank approximation becomes effective (Beckermann et al., 2016).

For normal $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$5 and $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$6, if $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$7 and $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$8, and the displacement rank is $$L[A,B](T)=T-ATB^T,\qquad A,B\in\mathbb{R}^{n\times n},$$9, then

$T$0

where $T$1 is a Zolotarev number defined by an extremal rational approximation problem. Since Cauchy matrices have displacement rank $T$2, this specializes to

$T$3

For separated intervals, explicit bounds show geometric decay with a rate governed by the absolute cross-ratio

$T$4

rather than by Euclidean distance alone (Beckermann et al., 2016).

This geometric decay yields explicit numerical-rank estimates. If

$T$5

then for Cauchy matrices

$T$6

A notable feature is that under fixed separation assumptions this bound depends on $T$7 and $T$8, but not directly on $T$9. This distinguishes Cauchy-type matrices from Krylov, real Vandermonde, and positive semidefinite Hankel matrices, for which the corresponding low-rank bounds scale like $D_t C - C D_s = GH^*,$00 (Beckermann et al., 2016).

6. Recognition, parameter recovery, and approximation

For the pure rank-one Cauchy case, the displacement equation can be inverted to recover the defining node sets directly from the entries. If

$D_t C - C D_s = GH^*,$01

then

$D_t C - C D_s = GH^*,$02

so the $D_t C - C D_s = GH^*,$03-displacement rank is $D_t C - C D_s = GH^*,$04. The equalities $D_t C - C D_s = GH^*,$05 are equivalent, when all entries are nonzero, to the linear system

$D_t C - C D_s = GH^*,$06

In matrix form,

$D_t C - C D_s = GH^*,$07

The kernel is one-dimensional, $D_t C - C D_s = GH^*,$08, so the Cauchy points are unique up to a global translation. Normalized Cauchy points are obtained by subtracting the mean shift, and exact recovery uses only the first row and first column. If the matrix is Cauchy, this yields normalized points in $D_t C - C D_s = GH^*,$09 operations, while deciding whether an arbitrary matrix is Cauchy requires $D_t C - C D_s = GH^*,$10 operations (Liesen et al., 2014).

When the matrix is only approximately Cauchy, one may solve the least-squares problem

$D_t C - C D_s = GH^*,$11

The associated structured matrix $D_t C - C D_s = GH^*,$12 has explicitly known nonzero singular values: $D_t C - C D_s = GH^*,$13 with multiplicity $D_t C - C D_s = GH^*,$14, $D_t C - C D_s = GH^*,$15 with multiplicity $D_t C - C D_s = GH^*,$16, and $D_t C - C D_s = GH^*,$17 with multiplicity $D_t C - C D_s = GH^*,$18. This leads to a closed-form minimum-norm least-squares solution computable in $D_t C - C D_s = GH^*,$19 time from row means, column means, and the global mean of $D_t C - C D_s = GH^*,$20. If the residual $D_t C - C D_s = GH^*,$21 satisfies $D_t C - C D_s = GH^*,$22, then the recovered parameters are valid Cauchy points and

$D_t C - C D_s = GH^*,$23

The same paper stresses that approximation quality and point-recovery accuracy are complementary rather than equivalent notions (Liesen et al., 2014).

A more recent recovery framework rewrites the problem in terms of the affine subspace

$D_t C - C D_s = GH^*,$24

with the reciprocal characterization

$D_t C - C D_s = GH^*,$25

This leads to projector-based algorithms

$D_t C - C D_s = GH^*,$26

which recover generators $D_t C - C D_s = GH^*,$27 from

$D_t C - C D_s = GH^*,$28

Special choices reproduce earlier methods: $D_t C - C D_s = GH^*,$29 gives the first-row/first-column recovery, while $D_t C - C D_s = GH^*,$30 gives the Frobenius-orthogonal projector onto $D_t C - C D_s = GH^*,$31. The general algorithm has $D_t C - C D_s = GH^*,$32 serial cost, and its approximation error is controlled by quasi-optimality bounds of the form

$D_t C - C D_s = GH^*,$33

with projector-dependent factor $D_t C - C D_s = GH^*,$34 (Boito et al., 16 Feb 2026).

The same work introduces a direct displacement-based criterion,

$D_t C - C D_s = GH^*,$35

which measures closeness to the exact Cauchy Sylvester equation without forming $D_t C - C D_s = GH^*,$36. In Frobenius norm this becomes a least-squares problem of cost $D_t C - C D_s = GH^*,$37. Numerical experiments reported there show that this displacement-based method is generally the most accurate when the goal is to approximate the observed matrix $D_t C - C D_s = GH^*,$38, whereas projector-based methods are cheaper and can be advantageous when perturbations are concentrated in specific rows or columns. This suggests a useful distinction between recovering a hidden exact Cauchy matrix and approximating noisy data by a Cauchy model (Boito et al., 16 Feb 2026).