🔢 mathematics

On the data-sparsity of the solution of Riccati equations with applications to feedback control

This paper demonstrates that solutions to large-scale continuous-time algebraic Riccati equations with quasiseparable coefficients inherit numerical quasiseparability, enabling the development of two efficient solvers for general and banded cases that are validated through applications in partial differential equation control and agent-based models.

Original authors: Stefano Massei, Luca Saluzzi

Published 2026-01-30

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Stefano Massei, Luca Saluzzi

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to solve a massive, complex puzzle to control a giant system, like a fleet of drones or the temperature of a whole city. In the world of mathematics, this puzzle is called a Riccati equation. Usually, solving this puzzle for huge systems is like trying to drink from a firehose; the amount of data is so large that standard computers choke on it, taking forever or running out of memory.

This paper introduces a clever trick: The solution to this puzzle is naturally "sparse" or "organized," even if the inputs look messy.

Here is the breakdown of the paper's ideas using everyday analogies:

1. The Hidden Order (Quasiseparability)

Think of the matrices (the grids of numbers) in these equations as a giant spreadsheet.

The Problem: Usually, we assume these spreadsheets are full of random numbers, making them impossible to compress.
The Discovery: The authors found that if the input spreadsheets have a specific structure (called quasiseparable), the solution spreadsheet also has a hidden structure.
The Analogy: Imagine a giant wall of lightbulbs. If the switches controlling them are arranged in a specific way, the pattern of lights that turns on isn't random. Instead, the lights in the corners of the wall are very dim (almost zero), and the brightness fades away as you move away from the center. This means you don't need to store the brightness of every single bulb; you only need to store the bright ones and a simple rule for how they fade. This "fading" property is what the authors call numerical quasiseparability.

2. The Two New Tools (Algorithms)

Because they discovered this hidden order, the authors built two new "machines" (algorithms) to solve the puzzle much faster.

Tool #1: The Divide-and-Conquer Chef (Algorithm 2)
- How it works: Imagine you have a giant, heavy pizza to eat. Instead of trying to eat it whole, you cut it in half, then cut those halves in half, and so on, until you have tiny, manageable slices. You solve the problem for the tiny slices and then glue the answers back together.
- The Magic: This tool uses a special format (HSS) that acts like a compression algorithm for the pizza slices. It allows the computer to handle massive problems by only focusing on the "important" parts of the data, ignoring the empty space.
Tool #2: The Pruning Gardener (Algorithm 3)
- How it works: This tool is designed for problems where the data is already somewhat organized (like a garden with rows of plants). It uses a method called "Newton-Kleinman" which is like taking a step toward the solution, checking how close you are, and taking another step.
- The Twist: As it takes steps, the data can get messy and grow too wide. The "Gardener" uses a trimming tool (thresholding). If a number is very small (like a tiny weed), it cuts it out and sets it to zero. This keeps the data "banded" (narrow and tidy) so the computer doesn't get overwhelmed.

3. Why This Matters (The Applications)

The paper tests these tools on two real-world scenarios where the "puzzle" is huge:

Controlling Fluids (The Allen-Cahn Equation): Imagine trying to control the temperature of a fluid flowing through a pipe to keep it from freezing or boiling. The math here is like a 1D or 2D grid of points. The new tools allowed the computer to calculate the perfect control signal in seconds, whereas old methods would take hours or fail completely.
Controlling a Swarm (The Cucker-Smale Model): Imagine a flock of birds or a swarm of robots trying to agree on a direction (consensus). Each agent has its own control. The math here involves a massive grid representing all the agents. The new tools successfully calculated how to steer the whole swarm to a stop or a specific formation efficiently.

4. The "Secret Sauce" (Mathematical Proof)

Before building the tools, the authors had to prove why the solution would be organized.

They used a concept called Zolotarev numbers (a fancy way of measuring how well you can approximate a curve with a simple fraction).
The Metaphor: They proved that the "fading" of the data (the dimming lightbulbs) happens so fast that you can approximate the solution with very few numbers. It's like saying, "Even though this wall of lights is huge, 99% of it is dark, so we only need to describe the 1% that is bright."

Summary

In short, this paper says: "Don't be scared of the size of the data. If the inputs have a certain structure, the answer will naturally be simple and organized. We built two fast tools to exploit this, allowing us to solve control problems for massive systems (like PDEs and agent swarms) that were previously too big to handle."

The paper does not claim these tools work for medical diagnoses or predicting stock markets; it strictly focuses on control theory (steering systems) and partial differential equations (modeling physical phenomena like heat or fluid flow).

Technical Summary: On the data-sparsity of the solution of Riccati equations with applications to feedback control

1. Problem Statement

The paper addresses the computational challenge of solving large-scale continuous-time algebraic Riccati equations (CAREs) of the form:
$R(X) := A^\top X + XA - XFX + Q = 0$
where $X \in \mathbb{R}^{n \times n}$ is the unknown stabilizing solution, and $F, Q$ are symmetric positive semidefinite matrices. While standard methods relying on Hamiltonian matrix decompositions are feasible for small $n$ , they become prohibitively expensive for large-scale systems ( $n \gg 1$ ).

The specific focus of this work is on scenarios where the coefficient matrices $A, F, Q$ possess quasiseparable structures (i.e., their off-diagonal blocks are low-rank) or are banded. The authors investigate whether the stabilizing solution $X$ inherits these structural properties, thereby enabling efficient approximation and computation. This is particularly relevant for infinite-horizon optimal control problems arising in the control of partial differential equations (PDEs) and agent-based models, where standard low-rank approximation methods (like RADI or rational Krylov subspace methods) often fail because the quadratic coefficient $F$ is full-rank rather than low-rank.

2. Methodology

The paper employs a dual approach combining theoretical analysis of solution structure with the development of two distinct numerical algorithms.

2.1 Theoretical Analysis: Decay of Singular Values

The authors establish that if $A, F, Q$ are quasiseparable, the off-diagonal blocks of the solution $X$ exhibit rapid decay in their singular values.

Connection to Sylvester Equations: The analysis relates generic off-diagonal blocks of $X$ to the solution of a Sylvester equation with a low-rank right-hand side.
Zolotarev Numbers: The decay rate is bounded by Zolotarev numbers, $Z_h(E, F)$ , which arise from rational approximation problems involving the numerical ranges or spectra of the closed-loop matrices.
Improved Bounds: The authors derive new decay bounds that improve upon existing results (specifically [41, Theorem 2.7]) for Sylvester equations with quasiseparable coefficients. These bounds are extended to the CARE setting, covering cases where $F$ is the identity, $F$ is full-rank quasiseparable, and $F$ is general symmetric positive definite.
Tensor Train (TT) Ranks: The decay results are used to derive tighter upper bounds for the Tensor Train ranks of the associated value function $V(y) = y^\top X y$ , improving upon previous estimates that were limited to low-rank $F$ .

2.2 Algorithmic Development

Two solvers are proposed based on the structural properties of the coefficients:

Algorithm 2: Divide-and-Conquer for General Quasiseparable Coefficients

Approach: This method utilizes a divide-and-conquer strategy based on the Hierarchically Semiseparable (HSS) matrix format.
Mechanism: The CARE is split into a block-diagonal part (solved recursively) and a low-rank correction part. The correction term $\delta X$ is approximated using an Extended Krylov Subspace Method (EKSM) projected onto a tensorized Krylov subspace.
Complexity: Under ideal assumptions, the complexity is $O(n \log n \cdot r^2)$ , where $r$ is the HSS rank.

Algorithm 3: Truncated Inexact Newton-Kleinman (tink) for Banded Coefficients

Approach: Tailored for the specific case where coefficients are banded. It employs an inexact Newton-Kleinman (NK) iteration.
Mechanism:
- The Lyapunov equations arising in the NK iteration are solved inexactly using the GMRES method.
- A thresholding mechanism is applied at each iteration to truncate the bandwidth of the iterate, preventing the bandwidth from growing uncontrollably.
- A line search strategy ensures the monotonic decrease of the Riccati residual and the stability of the closed-loop matrices.
Convergence: The authors provide convergence guarantees showing that with appropriate choices of inner iterations and truncation parameters, the sequence converges to the stabilizing solution.

3. Key Results

3.1 Theoretical Findings

Numerical Quasiseparability: The paper proves that the stabilizing solution of a CARE with quasiseparable coefficients is numerically quasiseparable, even when the quadratic coefficient $F$ is full-rank.
Improved Decay Bounds: The derived bounds for off-diagonal singular values are tighter than previous literature, particularly for the case of full-rank $F$ .
TT Rank Estimates: The analysis provides sharper estimates for the TT ranks of the value function, showing an $O(\log(\epsilon^{-1}))$ dependence on the approximation error $\epsilon$ , which is asymptotically tighter than previous $O(\log(\epsilon^{-1})^{7/2})$ bounds.

3.2 Numerical Experiments

Synthetic Tests:
- Algorithm 2: Demonstrated linear-logarithmic scaling with respect to matrix size $n$ and quadratic scaling with respect to the quasiseparable order. It achieved speed-ups of 65x to 700x compared to the dense MATLAB solver icare for large $n$ .
- Algorithm 3 (tink): Showed that for banded coefficients, the method is highly efficient. However, for ill-conditioned $F$ , the intermediate bandwidths in the NK iteration can grow, affecting performance. The "1-LS" variant (line search only at the first iteration) proved most efficient.
Application to PDEs (Allen-Cahn Equation):
- The methods were applied to the optimal control of the 1D and 2D Allen-Cahn equation.
- In 1D, tink achieved a 16x speed-up over dense evaluation of the exact solution formula while maintaining accuracy.
- In 2D, a nested dissection ordering combined with the divide-and-conquer solver (Algorithm 4) allowed for the solution of large-scale problems where dense solvers would require thousands of seconds. The computational cost per control step remained manageable as the grid resolution increased.
Application to Agent-Based Models (Cucker-Smale):
- The d&c_care algorithm was used to control interacting particle systems.
- Sorting the agent positions prior to solving significantly reduced the HSS rank of the solution, leading to lower CPU times. The method successfully stabilized the system to the consensus state.

4. Significance and Claims

The paper claims to provide the first theoretical results demonstrating the numerical quasiseparability of CARE solutions in the case of a full-rank quadratic coefficient $F$ . This theoretical insight justifies the use of hierarchical matrix formats for a broader class of control problems than previously understood.

The authors emphasize that their proposed algorithms enable the solution of large-scale optimal control problems (specifically those involving PDEs and agent-based models) that are currently intractable for state-of-the-art solvers like RADI or rational Krylov methods, which rely on low-rank $F$ . The work bridges the gap between theoretical decay properties and practical, scalable algorithms for feedback control synthesis.

5. Limitations and Future Directions

The authors modestly note several areas for future investigation:

The parameter $t$ in the decay bounds could potentially be improved by refining the splitting strategies (as suggested in Remarks 1 and 2).
The dependence on the condition number $\kappa(F)$ in the theoretical bounds appears to be an artifact of the proof, as numerical experiments suggest the decay behavior is robust even for large $\kappa(F)$ .
In the banded case, the tink algorithm sometimes exhibits a transient phase with large bandwidths before converging; alternative thresholding techniques to mitigate this are suggested as a potential improvement.