Reparameterized Low-Rank Updates

• Reparameterized low-rank updates are techniques that refine solutions to matrix equations by applying targeted low-rank corrections instead of full recomputation.
• They leverage rational and extended Krylov subspace methods to efficiently compress correction terms, ensuring rapid convergence and rigorous singular value decay guarantees.
• This approach underpins scalable solvers for Sylvester, Lyapunov, and Riccati equations, yielding significant speedups and memory savings in large-scale control and PDE applications.

Reparameterized Updates (Low-Rank)

A reparameterized low-rank update is a computational strategy in which the solution to a linear or nonlinear matrix equation—in particular, Sylvester, Lyapunov, and Riccati equations—is updated when the coefficients undergo structured low-rank changes. Rather than recomputing the full solution, the method targets only the correction, exploiting the low effective dimension induced by the perturbation. The approach leverages rational or extended block Krylov subspace methods for efficient large-scale computation, enables rigorous singular value decay guarantees, and underpins scalable algorithms for divide-and-conquer over hierarchical low-rank structures such as HODLR and HSS matrices. Reparameterized low-rank updating is now fundamental to large-scale control, PDE discretizations, and matrix function applications, yielding nearly linear time and memory complexity when appropriately embedded (Kressner et al., 2017).

1. Sylvester/Lyapunov Equation and Low-Rank Modeling

Consider the Sylvester equation

$A\,X \;+\; X\,B \;=\; C,$

where $A\in\mathbb{C}^{n\times n}$, $B\in\mathbb{C}^{m\times m}$, and $C\in\mathbb{C}^{n\times m}$, with $$\operatorname{spec}(A)\cap[-\operatorname{spec}(B)]=\varnothing$$ to guarantee unique solvability. The Lyapunov case is recovered for $B=A^*$, $C=C^* \preceq 0$, and Hurwitz $A$.

Suppose a reference solution $X_0$ to $A_0 X_0 + X_0 B_0 = C_0$ is known. A low-rank perturbation modifies the problem: $$A_0 \mapsto A_0 + \delta A,\quad B_0 \mapsto B_0 + \delta B,\quad C_0 \mapsto C_0 + \delta C,$$ with each update possessing factored form, e.g., $\delta A = U_A V_A^*$, $\delta B = U_B V_B^*$, $\delta C = U_C V_C^*$, where the ranks are much smaller than $\min\{n,m\}$.

2. Correction Equation and Right-Hand Side Compression

Subtracting the old from the new Sylvester equation and collecting terms yields a correction equation for $\delta X$: $$(A_0+\delta A)\,\delta X + \delta X\,(B_0+\delta B) = D,$$ where

$D = \delta C - \delta A X_0 - X_0 \delta B,$

with rank at most $$s = \operatorname{rank}(\delta A) + \operatorname{rank}(\delta B) + \operatorname{rank}(\delta C)$$.

A skinny factorization $D = UV^*$ with

$$U=[\,U_C,\, -U_A,\, -X_0 U_B\,], \quad V=[\,V_C,\, X_0^* V_A,\, V_B\,]$$

enables efficient compressed low-rank representation. A thin SVD can further reduce this to an optimal (numerically) effective rank $\tilde s$.

3. Rational and Extended Block Krylov Subspaces

To solve

$\widetilde{A}\, X + X\,\widetilde{B} = U V^*,$

with $U, V$ of small column dimension, tensorized block (rational/extended) Krylov subspaces are constructed:

• Left: $$\mathcal{U}_t = \mathrm{span}\{ U,\, \widetilde{A}^{-1} U,\, \widetilde{A} U,\, \widetilde{A}^{-2} U,\, \ldots \} \subseteq \mathbb{C}^n$$
• Right: $$\mathcal{V}_t = \mathrm{span}\{ V,\, \widetilde{B}^{-*} V,\, \widetilde{B}^* V,\, \ldots \} \subseteq \mathbb{C}^m$$

After $t$ iterations, orthonormal bases $U_t \in \mathbb{C}^{n\times \ell}$, $V_t \in \mathbb{C}^{m\times \ell}$ are formed ($\ell \approx 2s t$), yielding the approximation $X_t \approx U_t\,Y_t\,V_t^*$. The "compressed" Sylvester equation is small ($\ell \times \ell$) and solved directly, with convergence monitored via the residual norm.

4. Algorithmic Workflow, Computational Complexity, and Storage

The procedure divides into two stages:

1. Right-Hand Side Compression: Compute and compress $U, V$ factors on the updated right-hand side.
2. EKSM Solution: Apply the extended Krylov subspace method to produce the low-rank correction $\delta X$.

High-level pseudocode:

• Input: $(A_0,\delta A, B_0, \delta B, C_0, \delta C, X_0)$
• Compute $U, V$ so $D \approx UV^*$, recompress as needed.
• Solve $$(A_0+\delta A) \delta X + \delta X (B_0 + \delta B) = UV^*$$ via block EKSM.
• Return the updated solution $X_0 + \delta X$.

Complexity:

• For HODLR coefficients, a single EKSM invocation costs $\mathcal{O}(k^2 n \log^2 n)$ flops (or $\mathcal{O}(k^2 n \log n)$ for HSS).
• Memory requirements: storing $U_t, V_t$ in $\mathcal{O}(\ell n + \ell m)$; HODLR structure needs $\mathcal{O}(kn\log n)$.

5. Singular Value Decay and Stability Guarantees

Under mild spectral separation, the singular values of the update $\delta X$ decay rapidly. For the correction equation

$$A \delta X + \delta X B = D, \quad \operatorname{rank}(D) = s,$$

with numerical ranges $E = \mathcal{W}(A)$, $F = \mathcal{W}(-B)$ disjoint and compact, for any rational $r \in \text{Rat}_{h,h}$,

$$\sigma_{ks+1}(\delta X) \leq (1+\sqrt{2})^2 Z_h(E,F) \|\delta X\|,$$

where $Z_h(E, F)$ is the $h$th Zolotarev number and decays exponentially with $h$. For Hermitian positive definite $A,B$ with spectra in $[a,b]$, this yields exponential decay of $\sigma_j$ in $j$, enabling $\epsilon$-rank approximations with $h = \mathcal{O}(\log(1/\epsilon))$ (Kressner et al., 2017).

6. Embedding in Newton-Kleinman Iteration for Riccati Equations

In continuous-time CARE: $$A X + X A^* - X B X = C, \quad B \succeq 0, \quad A - X B\ \text{Hurwitz},$$ the Newton-Kleinman iteration solves a sequence of Lyapunov equations: $$(A - X_k B) X_{k+1} + X_{k+1} (A^* - B X_k) = C - X_k B X_k.$$ If $B$ is low-rank, each right-hand side difference remains low-rank. The low-rank update machinery specifically computes $\delta X_k$ as the solution with low-rank right-hand side, so subsequent Lyapunov solves become computationally cheap. Newton's quadratic convergence ensures only a few expensive solves are needed, the rest benefiting fully from the low-rank update formalism.

7. Empirical Benchmarks and Role in Divide-and-Conquer Solvers

Numerical experiments demonstrate:

• 5–10× speedup for large (up to $1.3\times10^5$) 2D Poisson/Lyapunov equations, with residuals ≤ $10^{-10}$.
• Memory reductions: HODLR from ≈0.5 GB to HSS ≈0.27 GB.
• For convection–diffusion and heat equations, 10× memory savings when compared to sparse-CG.
• For CAREs with banded + low-rank $B$, 3–5× speedup relative to plain Newton+EKSM.

Ranks remain small, often $\mathcal{O}(\log n)$, and residuals are driven to $10^{-8}$–$10^{-10}$. Performance gains are pronounced when many related Sylvester/Lyapunov (or CARE) solves share similarities apart from low-rank perturbations.

The methodology serves as the foundation for divide-and-conquer algorithms over hierarchical low-rank formats (HODLR, HSS), achieving nearly linear complexity in $n$ by recursively splitting and updating only the low-rank couplings (Kressner et al., 2017).

---

In summary, the reparameterized low-rank update approach for large matrix equations enables efficient, stable incremental solution updating under low-rank data changes. It achieves near-optimal storage and computational complexity, underpins fast solvers for hierarchical matrix structures, and integrates naturally into Newton-type methods for more complex nonlinear matrix equations. This framework is a core algorithmic building block for large-scale linear algebra in scientific computing and control (Kressner et al., 2017).