# ALGEBRAIC PROPERTIES OF SOLUTIONS TO CERTAIN NONLINEAR MATRIX EQUATIONS *

Valeria Simoncini

► To cite this version:

Valeria Simoncini. ALGEBRAIC PROPERTIES OF SOLUTIONS TO CERTAIN NONLINEAR MATRIX EQUATIONS *. 2025. hal-04975082

HAL Id: hal-04975082

https://hal.science/hal-04975082v1

Preprint submitted on 6 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
ALGEBRAIC PROPERTIES OF SOLUTIONS TO CERTAIN
NONLINEAR MATRIX EQUATIONS *

VALERIA SIMONCINI†

This paper is dedicated to Daniel Szyld, on the occasion of his 70th birthday.

**Abstract.** We are interested in algebraic properties of the solution $X$ to the linear or mildly nonlinear symmetric matrix equation $AX+XA+U\Phi(X)U^T+BB^T=0$ with symmetric $A$. We analyze monotonicity and low-rank properties of closed form solutions, whenever available, with respect to the solution of the equation $AX+XA+BB^T=0$. We extend this analysis to approximation recurrences, for which monotonicity and singular value decay properties are discussed.

**Key words.** Linear matrix equations, matrix iterations, singular value decay.

**AMS subject classifications.** 65F30

**1. Introduction.** We are interested in investigating certain structural properties of the solution $X$ to the following class of nonlinear equations,

$$ (1.1) \qquad AX + XA + U\Phi(X)U^T + BB^T = 0, $$

where the $n \times n$ matrix $A$ is assumed to be symmetric and negative or positive definite, and $U, B$ are tall matrices with full column rank. The function $\Phi$ stands for a linear or mildly nonlinear symmetric function depending on $X$, and is such that the dimensions of $\Phi(X)$ are conforming with those of $U$. In the following we shall denote with $\mathcal{M}: \mathbb{R}^{n \times n} \to \mathbb{R}^{n \times n}$ the coefficient operator, $\mathcal{M}(X) = AX + XA + U\Phi(X)U^T$, and with $\mathcal{L}: \mathbb{R}^{n \times n} \to \mathbb{R}^{n \times n}$ the portion of $\mathcal{M}$ corresponding to the Lyapunov operator, $\mathcal{L}(X) = AX + XA$.

The equation in (1.1) includes the problem

$$ (1.2) \qquad AX + XA + \sum_{k=1}^{\ell} N_k XN_k^T + BB^T = 0 $$

as a special case. This equation arises in the stability analysis of bilinear dynamical systems, see, e.g., [5],[9],[13],[16], where usually $N_k$ has rank much smaller than its dimension. The equation in (1.2) has recently emerged as a natural formulation in different contexts, beyond Control. This is the case, for instance, in problems where the unmixing of heterogeneous variables represents a computationally attractive solution strategy, such as in stochastic problems and space-time computations; see, e.g., [30]. The numerical solution of (1.2) is significantly more challenging than the Lyapunov equation $AX+XA+BB^T=0$, since in general obtaining a closed form solution matrix is computationally expensive. Because of this, (1.2) has recently attracted great interest in the numerical community, and different strategies have been proposed, see, e.g., [4],[10],[15],[28],[8],[23] and their references.

The generalization to (1.1) further broadens the class of problems that can be treated with similar procedures, including problems arising in different engineering

<footnote>*Version of March 3, 2025.</footnote>

<footnote>†Dipartimento di Matematica and (AM)², Alma Mater Studiorum - Università di Bologna, Piazza di Porta S. Donato, 5, I-40127 Bologna, Italy; and also IMATI-CNR, Pavia, Italy (valeria.simoncini@unibo.it).</footnote>
applications, see, e.g., [6],[12],[25]. While several successful computational advances have appeared in the recent literature, structural properties of the solution to (1.1) have been less analyzed. The fact that the solution is a matrix, and that it is symmetric, poses questions on structural properties that have been nicely addressed for the Lyapunov equation. In this paper we are thus interested in analyzing how properties such as definiteness and singular value decay of the solution matrix to the Lyapunov equation carry over to the setting in (1.1). Assessing these properties is crucial when targeting a low rank approximation to the exact solution. Our analysis exploits the closed form of the solution that can be derived for certain choices of $\Phi$. In particular, when available, a closed form highlights the composition of the solution in terms of two distinct components, from which the singular value decay can be readily deduced.

In case a closed form is not available or computationally unfeasible, we analyze sequences of approximate solutions, some of which are well established. For these approximate solutions, we study how definiteness and monotonicity are preserved throughout the iteration. In doing so, we also derive a new iteration for the nonlinear case $\Phi(X) = V^T(X \otimes X)V$. This iteration leverages recently introduced closed form solution techniques for quasi-linear equations of type (1.1), where $\Phi(X)$ is replaced by a real valued function in $X$.

An outline of the paper is as follows. In the subsections of section 2 we single out three distinct but highly connected classes of problems that fall in the framework of (1.1). In section 3.1 we summarize some known and less exercised structural properties of approximation sequences, together with a new recurrence for one of the discussed problems, which builds upon the construction methodology of other solutions. In section 3.2 we recall certain optimality properties of recently developed projection methods. In section 4 we derive monotonicity results for the discussed recurrences, their norms and their error norm. In section 5 we briefly discuss decay properties that can be inferred by the already derived structure. Section 6 contains our final considerations.

**1.1. Notation.** Given two matrices $A = (A_{i,j}) \in \mathbb{R}^{m \times n}$ and $B$, the Kronecker product is defined in block form as

$$A \otimes B = \begin{bmatrix} A_{1,1}B & \cdots & A_{1,n}B \\ \vdots & \ddots & \vdots \\ A_{m,1}B & \cdots & A_{m,n}B \end{bmatrix}.$$

This matrix operator satisfies

$$(1.3) \qquad \mathrm{vec}(AXB) = (B^T \otimes A)\mathrm{vec}(X),$$

where $\mathrm{vec}(X)$ stacks the columns of $X$ one below the other.

We use the notation $X \succeq 0$ to say that the square matrix $X$ is symmetric and positive semi-definite, while $X \succeq Y$ means that $X - Y \succeq 0$. We define the energy-norm associated with a symmetric and positive definite operator $\mathcal{A}$ as $\|X\|_{\mathcal{A}}^2 = \mathrm{trace}(X \mathcal{A}(X))$. Moreover, $\|X\|$ denotes the (spectral) matrix norm induced by the Euclidean vector norm, while $\|X\|_F$ denotes the Frobenius norm.

We specifically use bold face letters for matrices and vectors of leading dimension $n^2$, while calligraphic fonts are used for matrix operators.

**2. Problems with a quasi-linear matrix equation structure.**
**2.1. Multiterm linear matrix equations with low-rank structure.** We consider the problem (1.2), where for the sake of the presentation, we work with $\ell = 1$, so that the sum consists of a single term. Let $N \equiv N_1$ be written as $N = UV^T$ and assume that $U, V$ are full column rank matrices, with $r$ columns.

Using the Kronecker products, $\mathbf{A} = \mathbf{A} \otimes I + I \otimes \mathbf{A}$, $\mathbf{b} = \text{vec}(\mathbf{B}\mathbf{B}^T)$, and we can define $\mathbf{U} = \mathbf{U} \otimes \mathbf{U}$, $\mathbf{V} = \mathbf{V} \otimes \mathbf{V}$, so that $\text{rank}(\mathbf{U}) = r^2 = \text{rank}(\mathbf{V})$. Note that $\mathbf{U} = [\mathbf{u}_1, \dots, \mathbf{u}_{r^2}]$ with $\mathbf{u}_j = \mathbf{u}_k \otimes \mathbf{u}_i$ where $j = (k-1)r + i$. The matrix equation can be written in vectorized form as the linear system $(\mathbf{A} + \mathbf{U}\mathbf{V}^T)\mathbf{x} = \mathbf{f}$, with $\mathbf{f} = -\mathbf{b}$ and $\mathbf{x} = \text{vec}(\mathbf{X})$. If $r^2$ is still much smaller than $n$, then $\mathbf{x}$ can be determined by means of the Sherman-Morrison-Woodbury (SMW) formula ([11],[17],[29]) as follows

$$ \mathbf{x} = (\mathbf{A} + \mathbf{U}\mathbf{V}^T)^{-1}\mathbf{f} = \mathbf{A}^{-1}\mathbf{f} - \mathbf{A}^{-1}\mathbf{U}(\mathbf{I} + \mathbf{V}^T\mathbf{A}^{-1}\mathbf{U})^{-1}\mathbf{V}^T\mathbf{A}^{-1}\mathbf{f}. $$

The problem can be solved without using this vector formulation, but solely relying on matrix equation solves. This key property was first observed and algorithmically developed in [10], further explored later in [4], and more recently in [15]. Let

$$ (2.1) \qquad \mathbf{w} = \mathbf{A}^{-1}\mathbf{f} \Leftrightarrow \mathbf{A}\mathbf{W} + \mathbf{W}\mathbf{A} = -\mathbf{B}\mathbf{B}^T, $$

and note that $\mathbf{W}$ is symmetric. Let $\mathbf{U} = [\mathbf{u}_1, \dots, \mathbf{u}_r]$, $\mathbf{V} = [\mathbf{v}_1, \dots, \mathbf{v}_r]$ and again let $\mathbf{u}_j = \text{vec}(\mathbf{u}_i\mathbf{u}_k^T)$ for $k, i = 1, \dots, r$ such that $j = (k-1)r + i$. Then

$$ (2.2) \qquad \mathbf{p}_j = \mathbf{A}^{-1}\mathbf{u}_j \Leftrightarrow \mathbf{A}\mathbf{P}_j + \mathbf{P}_j\mathbf{A} = \mathbf{u}_i\mathbf{u}_k^T, \quad \mathbf{p}_j = \text{vec}(\mathbf{P}_j). $$

Note that for each $j_1$ there exists a $j_2$ such that $\mathbf{P}_{j_1} = \mathbf{P}_{j_2}^T$, so that not all $\mathbf{P}_j$ need to be explicitly computed. The inner system has coefficient matrix

$$ (2.3) \qquad \mathbf{H} = \mathbf{I} + \mathbf{V}^T\mathbf{A}^{-1}\mathbf{U} = \mathbf{I} + \mathbf{V}^T\mathbf{P} = \mathbf{I} + \mathbf{V}^T[\mathbf{p}_1, \dots, \mathbf{p}_{s^2}], $$

which can equivalently be computed by noticing that $\mathbf{v}_j^T\mathbf{p}_t = v_i^T\mathbf{P}_t\mathbf{v}_k$, $j = (k-1)r + i$. The solution of the $r^2 \times r^2$ linear system

$$ (2.4) \qquad \mathbf{H}\mathbf{g} = \mathbf{V}^T\mathbf{w} $$

requires first computing the right-hand side. This can be obtained in terms of the original data and of (2.1) as $\mathbf{V}^T\mathbf{A}^{-1}\mathbf{f} = [v_1^T\mathbf{W}\mathbf{v}_1, v_2^T\mathbf{W}\mathbf{v}_1, \dots, v_r^T\mathbf{W}\mathbf{v}_r]^T$. The final solution is then given by $\mathbf{X} = \mathbf{W} - \sum_{j=1}^{r^2} \mathbf{P}_j(g)_j$, where $(g)_j$ is the $j$th component of the vector $\mathbf{g}$. The effective implementation of this procedure requires some care; we refer to [15] and its references for computational details.

The solution $\mathbf{X}$ is expressed as a linear combination of solutions to Lyapunov equations having low-rank known term. The following result emphasizes this fact while summarizing the above derivation. The proof consists of a rewriting of the closed form above.

**PROPOSITION 2.1.** Let $\mathcal{L}: X \to AX + XA$. Let $\mathcal{G}$ be such that $\mathcal{G} = \text{vec}(\mathcal{G})$ where $\mathcal{G}$ is the solution to (2.4), with $\mathcal{H}$ defined in (2.3). Then the solution $\mathbf{X}$ to (6.1) satisfies $\mathcal{L}(\mathbf{X}) = -\mathbf{B}\mathbf{B}^T - \mathbf{U}\mathbf{G}\mathbf{U}^T$.

We notice that since $\mathbf{X}$, $\mathcal{L}(\mathbf{X})$ and $\mathbf{B}\mathbf{B}^T$ are symmetric, then also $\mathcal{G}$ must be symmetric, so that $\mathbf{X} = \mathbf{W} - \sum_{j=1}^{r(r+1)/2}(\mathbf{P}_j + \mathbf{P}_j^T)(g)_j$, with the indexes properly ordered.
**2.2. Quasi-linear matrix equations.** A special setting is determined by the problem

$$AX + XA + f(X)C + BB^T = 0,$$

where the real valued function $f$ can be either a linear or a nonlinear function of $X$; a typical linear example is $f(X) = \text{trace}(GX)$ for some matrix $G$ [25]. This problem corresponds to (1.1) with $\Phi(X) := f(X)I$ and $C = UU^T$. For linear $f$, it was shown in [25] that the solution $X$ can be obtained in closed form as follows. Let $M$ be the solution to $AX+XA+BB^T = 0$ and $N$ be the solution to $AX+XA+C = 0$. Then, for $1-f(N) \neq 0$,

$$ (2.5) \qquad X = M + \sigma N, \qquad \sigma = \frac{f(M)}{1 - f(N)}. $$

The result can be generalized to more terms $f_i(X)C_i$, $i = 1, \dots, s$.

Similarly to the Sherman-Morrison-Woodbury formula of section 2.1, the solution is split into the sum of solutions to distinct Lyapunov equations. This is not surprising, since for $f(X) = v^T X v$ and $C = uu^T$, the two problems are equivalent [25].

**2.3. Quadratic-bilinear matrix equations.** A quite significant generalization of the linear equation in the previous section includes a special quadratic term in the unknown matrix. More precisely, the problem can be written as<fnref n="1" />

$$ (2.6) \qquad AX + XA + \mathbb{H}^T(X \otimes X)\mathbb{H} + \sum_{j=1}^{\ell} N_j X N_j^T + BB^T = 0. $$

with $\mathbb{H} \in \mathbb{R}^{n^2 \times n}$. Finding a definite solution $X$ provides information on the controllability of the associated quadratic-bilinear dynamical system; see, e.g., [6],[12] and references therein.

In the following, we show that if $\mathbb{H}$ has low rank, then the whole problem can be formalized as (1.1) with a simplified structure. Since the novelty of the equation resides in the quadratic term, in the following we assume that $N_j = 0$. We should say that the application studied in [6] does not lead to a matrix $\mathbb{H}$ of low rank, so that our discussion does not apply to that case. Nonetheless, we find the derivation interesting, and it may provide insight for the full-rank setting.

We first notice that for each column $h_j$ of $\mathbb{H}$, it holds that

$$ h_i^T (X \otimes X) h_j = \text{trace}(X H_j X^T H_i^T), \quad \text{with } H_i \text{ s.t. } h_i = \text{vec}(H_i). $$

A matrix $\mathbb{H}$ of low rank carries convenient features. Setting $\mathbb{H} = VU^T$ with $U \in \mathbb{R}^{n \times r}$, $V \in \mathbb{R}^{n^2 \times r}$, we have

$$ \mathbb{H}^T (X \otimes X) \mathbb{H} = U V^T (X \otimes X) V U^T, $$

which corresponds to the term $U\Phi(X)U^T$ in (1.1) with $\Phi(X) = V^T(X \otimes X)V$. Using these relations, and assuming first that $\mathbb{H}$ has rank equal to one (that is, $r = 1$), the following proposition provides a closed form for the sought after solutions.

<footnote n="1">¹We state the problem for symmetric $A$ for consistency with respect to our setting, however the whole derivation holds in the nonsymmetric case as well.</footnote>
PROPOSITION 2.2. If $\mathbb{H} = \mathbf{v}u^T$ has rank equal to one, with $\mathbf{v} \in \mathbb{R}^{n^2}$ and $u \in \mathbb{R}^n$, then (2.6) can be written as

$$AX + XA + \chi(X)uu^T + BB^T = 0, \quad (2.7)$$

with $\chi(X) = \text{trace}(XSX^TS^T)$, $\mathbf{v} = \text{vec}(S)$. Let $M, N$ be solutions to $AX + XA + BB^T = 0$ and $AX + XA + uu^T = 0$, respectively. If $\alpha := \text{trace}(NSNS^T) \neq 0$, then equation (2.7) admits the following two solutions,

$$X_1 = M + \chi_1 N, \quad X_2 = M + \chi_2 N,$$

where $\chi_{1,2} = \frac{1}{2\alpha}(-\beta \pm \sqrt{\beta^2 - 4\alpha\gamma})$, with $\beta = \text{trace}(MSNS^T + NSMS^T) - 1$ and $\gamma = \text{trace}(MS\bar{M}S^T)$.

Before we proceed with the proof, we notice that the way we have expressed the closed form solution, yields information on the existence of real solutions, which are obtained only for real $\chi_{1,2}$. Moreover, for real $\chi$, both solutions are symmetric, as expected. Definiteness can be easily detected once the sign of $\mathcal{L}$ is known, by associating the corresponding choice of $\chi_{1,2}$.

*Proof.* We have

$$\mathbb{H}^T(X \otimes X)\mathbb{H} = u\mathbf{v}^T(X \otimes X)\mathbf{v}u^T = u\text{trace}(XSX^TS^T)u^T = u\chi(X)u^T.$$

Hence, the original problem can be written as in (2.7). The problem is now in the form discussed in section 2.2, with $\chi$ nonlinear, and we can proceed analogously. More precisely, we first notice that

$$X = M + \chi(X)N,$$

with $X$ symmetric. We multiply by $S$ and by $S^T$, and then multiply the resulting matrices so as to obtain

$$\begin{aligned} (XS)(XS^T) &= (MS + \chi(X)NS)(MS^T + \chi(X)NS^T) \\ &= MSMS^T + \chi(X)MSNS^T + \chi(X)NSMS^T + \chi(X)^2NSNS^T. \end{aligned}$$

Taking the trace,

$$\chi(X) = \text{trace}(MSMS^T) + \chi(X)\text{trace}(MSNS^T + NSMS^T) + \chi(X)^2\text{trace}(NSNS^T).$$

Setting $\alpha, \beta$ and $\gamma$ like in the statement, we obtain the quadratic algebraic equation $\alpha\chi^2 + \beta\chi + \gamma = 0$, whose solutions are the given $\chi_{1,2}$. $\square$

The generalization to $\mathbb{H}$ of (small) rank greater than one gives the more general form in (1.1). We write $\mathbb{H} = \mathbf{V}U^T$ with $\mathbf{V} = [\mathbf{v}_1, \dots, \mathbf{v}_r]$ and $U = [u_1, \dots, u_r]$. Then

$$\mathbb{H}^T(X \otimes X)\mathbb{H} = U\mathbf{V}^T(X \otimes X)\mathbf{V}U^T = U\Phi(X)U^T,$$

where the entries of the $r \times r$ matrix $\Phi(X)$ are given by

$$(\Phi(X))_{i,j} = \text{trace}(XS_jX^T S_i^T), \quad \text{vec}(S_k) = \mathbf{v}_k, \quad k = i, j,$$

so that $\Phi$ is nonlinear in $X$. Thus, we can write

$$AX + XA^T + U\Phi(X)U^T + BB^T = 0,$$

and proceed like in section 2.1, leading to a feasible computational strategy as long as $U$ has few columns.

We conclude this section by saying that disregarding the rank of $\mathbb{H}$, a general closed form solution can be described by means of an integral series, though this expression does not seem to lead to a computationally feasible strategy [6].
**3. Sequences of approximations.** Except when a computable closed form can be obtained, an approximation to the solution of (1.1) is derived by generating a sequence ${X^{(k)}}_{k \ge 0}$ approaching $X$ as $k$ goes to infinity. In particular, whenever the number of columns $r$ is not very small, say a few tens, also the Sherman-Morrison-Woodbury formula for $\Phi(X) = V^T XV$ becomes infeasible, and an approximation recurrence is required.

**3.1. Fixed point iterations.** A classical approach is given by a fixed point iteration. For the generic case (1.1), this can be written as

$$AX^{(k+1)} + X^{(k+1)}A = -U\Phi(X^{(k)})U^T - BB^T,$$

leading to a sequence of solves with Lyapunov equations, that is,

$$X^{(k+1)} = \mathcal{L}^{-1}(-U\Phi(X^{(k)})U^T - BB^T);$$

see, e.g., [10][28]. Computationally, the problem then becomes that of efficiently solving a sequence of Lyapunov equations, which can be a challenge for large problems, especially considering that the Lyapunov equations cannot be solved at high accuracy. This procedure gives rise to inexact iterates, for which the convergence analysis needs to be revisited, and associated with the accuracy threshold of the Lyapunov solver; we refer to [28] for this discussion. The idea is very general, so that the same approach can be used in a number of linear and nonlinear matrix equations, such as the algebraic Riccati equation [7]; see [18] for other examples related to the approximation of matrix functions. Ad-hoc iterations can be derived for specific problems, see below.

Fixed point iterations have also been proposed in [6] for solving the Quadratic-bilinear equation in section 2.3 in the large scale setting. More precisely, in [6, (41)] the following recurrence is proposed,

$$\mathcal{L}(X^{(k+1)}) = -BB^T - \mathbb{H}^T(X^{(k)} \otimes X^{(k)})\mathbb{H} - \sum_{j=1}^{\ell} U_j \Phi_j(X^{(k)}) U_j^T.$$

This general iteration is very relevant for $\mathbb{H}$ having large or column full-rank. On the other hand, for $\mathbb{H}$ having low rank, a possibly more effective fixed point iteration can be derived. To simplify the presentation, we set $\Phi_j(X^{(k)}) = 0$ for all $j$ and focus on the quadratic term, which makes the problem special.

If $\mathbb{H}$ has rank equal to one, Proposition 2.2 ensures a closed form solution. If $\mathbb{H}$ is low-rank, that is $\mathbb{H} = VU^T$ with $U = [u_1, \dots, u_r] \in \mathbb{R}^{n \times r}$ and $V = [v_1, \dots, v_r] \in \mathbb{R}^{n \times r}$, then the original equation can be written as

$$(3.1) \quad AX + XA^T + \sum_{i,j=1}^{r} u_i \Phi(X)_{i,j} u_j^T + BB^T = 0, \quad \text{with} \quad \Phi(X)_{i,j} = v_i^T (X \otimes X) v_j.$$

We recall that $\Phi(X)_{i,j} = \text{trace}(XV_j X^T V_i^T)$, where $v_i = \text{vec}(V_i)$. We can thus define a linearized fixed point iteration as follows. We replace $\Phi(X)$ with $\Psi_{i,j}(X^{(k)}, X) = \text{trace}(XV_j(X^{(k)})^T V_i^T)$. For a computed $X^{(k)}$, we set $K_{i,j}^{(k)} := V_j(X^{(k)})^T V_i^T$ so that

$$\Psi_{i,j}(X^{(k)}, X) = \text{trace}(X K_{i,j}^{(k)}),$$

which is now linear in $X$. The iteration (3.1) is transformed into the quasi-linear problem

$$(3.2) \quad \mathcal{L}(X^{(k+1)}) = -BB^T - \sum_{i,j=1}^{r} \text{trace}(X^{(k+1)} K_{i,j}^{(k)}) u_i u_j^T,$$
in the unknown matrix $X^{(k+1)}$. We can obtain $X^{(k+1)}$ by a procedure similar to that recalled in section 2.2. Indeed, let $P_{\{i,j\}} = \mathcal{L}^{-1}(u_i u_j^T) \in \mathbb{R}^{n \times n}$, for $i,j = 1, \dots, r$. and let $M$ solve $AX + XA + BB^T = 0$. Then

$$X + \sum_{i,j=1}^{r} \omega_{i,j}^{(k)} P_{\{i,j\}} = M, \quad \omega_{i,j}^{(k)} = \text{trace}(X K_{i,j}^{(k)}).$$

For each $\ell, s = 1, \dots, r$, multiply the equation by $K_{\ell,s}^{(k)}$ and take the trace, so as to obtain

$$\text{trace}(X K_{\ell,s}^{(k)}) + \sum_{i,j=1}^{r} \omega_{i,j}^{(k)} \text{trace}(P_{\{i,j\}} K_{\ell,s}^{(k)}) = \text{trace}(M K_{\ell,s}^{(k)}),$$

that is,

$$\omega_{\ell,s}^{(k)} + \sum_{i,j=1,\dots,r} \omega_{i,j}^{(k)} \text{trace}(P_{\{i,j\}} K_{\ell,s}^{(k)}) = \text{trace}(M K_{\ell,s}^{(k)}).$$

Denote by $\tau_{i,j}^{(\ell,s)} = \text{trace}(P_{\{i,j\}} K_{\ell,s}^{(k)})$ the known coefficients, and let

$$\omega^{(k)} = \text{vec}(\text{trace}(X K_{\ell,s}^{(k)}))_{\ell,s=1,\dots,r}, \quad t^{(k)} = \text{vec}(\text{trace}(M K_{\ell,s}^{(k)}))_{\ell,s=1,\dots,r}$$

be the unknown traces and right-hand side vector of length $r^2$, respectively. This yields the linear system

$$(I_{r^2} + T)\omega^{(k)} = t^{(k)}, \quad T = \begin{bmatrix} 1 + \tau_{1,1}^{1,1} & \tau_{2,1}^{1,1} & \cdots & \tau_{1,2}^{1,1} & \cdots & \tau_{r,r}^{1,1} \\ \tau_{1,1}^{2,1} & 1 + \tau_{2,1}^{2,1} & \cdots & \tau_{1,2}^{2,1} & \cdots & \tau_{r,r}^{2,1} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ \tau_{1,1}^{1,2} & \tau_{2,1}^{1,2} & \cdots & \tau_{1,2}^{1,2} & \cdots & \tau_{r,r}^{1,2} \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ \tau_{1,1}^{r,r} & \tau_{2,1}^{r,r} & \cdots & \tau_{1,2}^{r,r} & \cdots & 1 + \tau_{r,r}^{r,r} \end{bmatrix}.$$

Let $\Omega^{(k)}$ be the $r \times r$ matrix reshaping of the solution $\omega^{(k)}$. Then the next iterate $X^{(k+1)}$ is obtained as

$$X^{(k+1)} = M - \sum_{i,j=1}^{r} (\Omega^{(k)})_{i,j} P_{\{i,j\}}. \quad (3.3)$$

As a computational side remark, we notice that only $\Omega^{(k)}$ changes at each iteration, and that all other matrices could be computed once for all before starting the recurrence.

**EXAMPLE 3.1.** In this example we report on our experience with the iteration (3.3), compared with the standard procedure in (3.1). To this end we consider a built-up, fully reproducible example, where we make the problem less definite by changing a problem parameter. The matrix $A$ is given as $A = I \otimes T + T \otimes I$ of size $n \times n$ with $n = 2500$, where $T = n_h^2 \text{tridiag}(1, -2, 1)$, with $n_h^2 = n$. The nonlinear term includes $\mathbb{H} = \alpha V U^T$ with both $U, V$ chosen randomly with $r = 4$ columns, and $\alpha$ controls the definiteness of the operator $\mathcal{M}$. The known term $E = BB^T$ is a random positive definite matrix. These definitions correspond to the following matlab ([22]) commands
<table><caption>TABLE 1</caption><thead><tr><th rowspan="2">α</th><th colspan="2">fixed point</th><th colspan="2">quasi-linear</th></tr><tr><th># its</th><th>CPU time</th><th># its</th><th>CPU time</th></tr></thead><tbody><tr><td>4</td><td>17</td><td>18.88</td><td>10</td><td>35.00</td></tr><tr><td>5</td><td>42</td><td>36.40</td><td>13</td><td>43.64</td></tr><tr><td>5.5</td><td>188</td><td>147.13</td><td>19</td><td>60.64</td></tr><tr><td>6</td><td>—</td><td>—</td><td>88</td><td>309.27</td></tr></tbody></table>

Example 3.1. Performance of standard fixed point iteration and quasi-linear variant. Time is in seconds. ‘-’ stands for divergence.

204
205 alpha=6;
206 nh=40; n=nh*nh; r=4;
207 T=toeplitz([-2,1,zeros(1,nh-2)]); I=speye(nh);
208 A=n*(kron(I,T)+kron(T,I));
209 rng(1);
210 E=rand(n,n);E=E*E';E=E/norm(E);
211 V=alpha*randn(n*n,r);
212 U=randn(n,r);

The stopping criterion was

$$ \frac{\|X^{(k+1)} - X^{(k)}\|_F}{\|BB^T\|_F} < 10^{-6}, $$

213 and a maximum of 200 iterations was considered. At completion the true residual 214 norm was also computed.

215 Table 1 reports the performance of both iterations for $\alpha \in [4, 6]$. For small $\alpha$, the 216 standard iteration (3.1) is faster in terms of CPU time, whereas the new approach 217 is able to get good information on the solution in fewer iterations. For the second 218 value of $\alpha$, the gap in the number of iterations increases, although the CPU time of 219 the new iteration remains larger. After that, the standard fixed point iteration starts 220 showing convergence difficulties, until a blow up for $\alpha = 6$, while good properties 221 of (3.1) can still be observed. Although not reported here, we notice that up to 222 $\alpha = 5.5$ convergence of the iteration in (3.3) is monotonic, in terms of $\|X^{(k+1)} - X^{(k)}\|$ 223 whereas for $\alpha = 6$ oscillations occur. Larger values of $\alpha$ than those considered here 224 will eventually lead to lack of convergence also for the new recurrence.

225 We consider these results very promising. However, assessing the overall quality 226 of the new method requires a deeper investigation, which is beyond the scope of this 227 work. ◇

The previous examples all generate matrix sequences {$X^{(k)}$}$_{k \ge 0}$ such that

$$ AX^{(k+1)} + X^{(k+1)}A = -BB^T - \sum_{i,j} u_i w_{i,j}^{(k)} u_j^T $$

228 for some $w_{i,j}^{(k)}$, so that

$$ (3.4) \qquad X^{(k+1)} = -\mathcal{L}^{-1}(BB^T) - \sum_{i,j} w_{i,j}^{(k)} \mathcal{L}^{-1}(u_i u_j^T). $$
We stress once again that pairs of solutions of $\mathcal{L}^{-1}(u_i u_j^T)$ are the transpose of each other, and that overall the sum is symmetric; see the similar discussion in section 2.1.

The solution matrix is a linear combination of solutions to distinct Lyapunov equations, where the first term is the solution one would obtain with $\Phi \equiv 0$. Hence, inspecting the magnitude of the coefficients $w_{i,j}^{(k)}$ enables one to appreciate the influence on the solution of the extra term in (1.1) compared with the plain Lyapunov problem.

**3.2. Projection methods.** An alternative to fixed point iterations is given by projection-type methods for the original equation (1.1), which have been mostly explored in the case of linear matrix equations, that is for $\Phi(X) = V^T X V$. Given an approximation space of dimension $s_k$ and a matrix $W_k$ with orthonormal columns with the same range, an approximate solution is obtained as $X^{(k)} = W_k Y^{(k)} W_k^T$, where $Y^{(k)} \in \mathbb{R}^{s_k \times s_k}$ is determined, e.g., by imposing an orthogonality (Galerkin) condition to the residual $R^{(k)} = \mathcal{M}(X^{(k)}) + BB^T$ [30]. Such condition can be formalized as requiring that $W_k^T R^{(k)} W_k = 0$. It can be shown, see [23] and its references, that for linear and symmetric positive definite $\mathcal{M}$ in (1.1), the approximate solution $X^{(k)}$ obtained with the Galerkin condition minimizes the error in the energy-norm, that is

$$ (3.5) \qquad X^{(k)} = \arg \min_{\substack{Z=W_k \tau W_k^T \\ \tau \in \mathbb{R}^{s_k \times s_k}}} \|X-Z\|_{\mathcal{M}}. $$

This result will be used in the next section. In [23] a two-term-type recurrence was also derived for this approach, so that the approximation can be written as

$$ X^{(k+1)} = X^{(k)} + P_k \alpha^{(k)} P_k^T. $$

This solution update generalizes the well known two-term recurrence of the Conjugate Gradient method for vector linear systems. In spite of the similarity with respect to the previous recurrences, here $P_k$ has low rank, which grows at each iteration, and the dimensions of the matrix $\alpha_k$ grow accordingly. As opposed to (stationary) fixed point iterations, this matrix update is nonstationary, since $\alpha_k$ and $P_k$ are chosen at each iteration so as to satisfy certain optimality properties; see [23] for more details.

**4. Monotonicity results.** In this section we are concerned with monotonicity properties of the solution matrix or of its approximation iterates. These may be expressed in terms of definiteness, or in terms of some norm.

If $A$ is symmetric and negative definite, we have that the Lyapunov solution $X_{\mathcal{L}} = \mathcal{L}^{-1}(-BB^T)$ is positive semi-definite, and that the solution $X$ to (1.1) satisfies

$$ X \succeq X_{\mathcal{L}} \quad \text{if and only if} \quad \Phi(X) \succeq 0. $$

Positive definiteness can be checked in all instances where we have derived a closed form of the solution and it readily depends on the involved scalars (see, e.g., (2.5)), whereas the property $X \succeq X_{\mathcal{L}}$ needs to be derived case by case.

For the sequences $\{X^{(k)}\}_{k \ge 0}$ we can show a similar result.

**PROPOSITION 4.1.** Let the problem $AX + XA + U\Phi(X)U^T + BB^T = 0$ be given. In the iteration (3.4), assume that $\Phi$ is linear and such that $\Phi(Y) \succeq 0$ for all $Y \succeq 0$. Assume that $X^{(0)} = 0$. Then,

(i) If $\mathcal{L}$ is negative definite, then for any $k \ge 0$, it holds that $X^{(k+1)} \succeq X^{(k)} \succeq 0$;
(ii) If $\mathcal{L}$ is positive definite, the semidefiniteness alternates between successive iterates.

Proof. We notice that $\mathcal{L}(X^{(k+1)}) = -BB^T - U\Phi(X^{(k)})U^T$, so that
$$X^{(k+1)} - X^{(k)} = -\mathcal{L}^{-1}(U\Phi(X^{(k)} - X^{(k-1)})U^T).$$

Hence, the definiteness sign of $X^{(k+1)} - X^{(k)}$ depends on that of $\mathcal{L}$ and on the sign of the approximate solutions difference at the previous step $k$. i) Taking $X^{(0)} = 0$ ensures $X^{(1)} \succeq 0$, from which all subsequent terms satisfy $X^{(k+1)} - X^{(k)} \succeq 0$. ii) Taking $X^{(0)} = 0$ ensures $0 \succeq X^{(1)}$ and from there on, the signs alternate. $\square$

For given symmetric matrices $Z_1, Z_2$ such that $Z_1 \succeq Z_2$ it holds that $\lambda_i(Z_1) \ge \lambda_i(Z_2)$ where $\lambda_i(\cdot)$ are the increasingly ordered eigenvalues of the argument matrix [19, Corollary 7.7.4]. Hence, the result of Proposition 4.1(i) implies the following norm inequalities, under the same hypotheses

$$\|X^{(k+1)}\| \ge \|X^{(k)}\|, \quad \|X^{(k+1)}\|_F \ge \|X^{(k)}\|_F.$$

In the linear setting, by just reformulating classical results, we can explicitly write a strictly decreasing bound for the error. We consider the linear case $AX + XA + UV^T XVU^T + BB^T = 0$, which in Kronecker form corresponds to

$$Ax + b = 0 \quad \text{with} \quad A = L + N,$$

where $L = A \otimes I + I \otimes A$ and $N = (U \otimes U)(V \otimes V)^T$. Let $x^*$ be the exact solution to the linear system above. Given the recurrence $Lx^{(k+1)} = -b - Nx^{(k)}$, the following result provides a bound for the error $e^{(k)} = x^* - x^{(k)}$ for a stationary iteration. We include the proof for completeness, though it may be derived similarly to other proofs in the literature$^2$, [26, Corollary 4.1].

**PROPOSITION 4.2.** Let $A = L + N$ with both $A$ and $L$ symmetric and positive definite. If $\rho(L^{-1}A) < 2$, then

$$\|e^{(k+1)}\|_{A} \le \rho(I - L^{-1}A)\|e^{(k)}\|_{A}.$$

Proof. Let $B = L^{-1}N = I - L^{-1}A$. We first observe that $\hat{B} := A^{\frac{1}{2}}BA^{-\frac{1}{2}} = I - A^{\frac{1}{2}}L^{-1}A^{\frac{1}{2}}$ is symmetric and that thanks to the similarity transformations, $\rho(\hat{B}) = \rho(B)$. Hence, using $\|e^{(k+1)}\|_A = \|Be^{(k)}\|_A$, we have

$$\begin{aligned} \|e^{(k+1)}\|_{A} &= \|A^{\frac{1}{2}}Be^{(k)}\| = \|A^{\frac{1}{2}}BA^{-\frac{1}{2}}A^{\frac{1}{2}}e^{(k+1)}\| \\ &\le \|\hat{B}\| \|e^{(k+1)}\|_{A} = \rho(B)\|e^{(k)}\|_{A}. \end{aligned} \quad \square$$

In our setting, $L$ corresponds to $\mathcal{L}$, and the result above reads

$$\|E^{(k+1)}\|_{\mathcal{M}} \le \rho(I - \mathcal{L}^{-1}(\mathcal{M}(\cdot)))\|E^{(k)}\|_{\mathcal{M}},$$

where $\mathcal{M}: X \mapsto AX + XA + U\Phi(X)U^T$ is symmetric and positive definite.

Projection methods satisfy similar inequalities, though stemming from quite different properties. More precisely, the minimization property in (3.5) ensures monotonicity of the approximate solution obtained via a projection method. The following result makes this statement more precise, under the only hypothesis that $\mathcal{M}$ be symmetric and positive definite.

<footnote>$^2$We thank Silvia Noschese for pointing us to [26] for this type of results.</footnote>
**PROPOSITION 4.3.** Let $\mathcal{M}$ be symmetric and positive definite. Then the Galerkin approximate solution $X^{(k)} = W_k Y^{(k)} W_k^T$ to (1.1) converges monotonically increasing towards $X$ in the $\mathcal{M}$-norm, that is

$$ \|X^{(k+1)}\|_{\mathcal{M}} \geq \|X^{(k)}\|_{\mathcal{M}}. $$

288 Proof. Let $E^{(k)} = X - X^{(k)}$ and $R^{(k)} = \mathcal{M}(X^{(k)}) + BB^T$. Observe that $R^{(k)} = -\mathcal{M}(E^{(k)})$ and that $\text{trace}(\mathcal{M}(Y)X) = \text{trace}(Y\mathcal{M}(X))$ for symmetric matrices $X, Y$ of conforming dimensions. From (3.5) we have

$$ (4.1) \qquad \|E^{(k+1)}\|_{\mathcal{M}} \leq \|E^{(k)}\|_{\mathcal{M}}. $$

292 We have $\|E^{(k)}\|_{\mathcal{M}}^2 = \text{trace}(E^{(k)}\mathcal{M}(E^{(k)})) = -\text{trace}(E^{(k)}R^{(k)}).$ The Galerkin condition on the residual ensures that $\text{trace}(X^{(k)}R^{(k)}) = 0$, since $\text{trace}(X^{(k)}R^{(k)}) = \text{trace}(Y^{(k)}W_k^T R^{(k)}W_k)$ and $W_k^T R^{(k)}W_k = 0$. Then,

$$ \begin{aligned} \text{trace}(E^{(k)} R^{(k)}) &= \text{trace}(X R^{(k)}) = -\text{trace}(\mathcal{M}(X) E^{(k)}) \\ &= -\text{trace}(\mathcal{M}(X)X) + \text{trace}(\mathcal{M}(X)X^{(k)}). \end{aligned} $$

294 We recall that $\text{trace}(\mathcal{M}(X)X^{(k)}) = \text{trace}(-BB^T X^{(k)}) = \text{trace}(\mathcal{M}(X^{(k)})X^{(k)})$, where once again we used that $\text{trace}(X^{(k)}R^{(k)}) = 0$. Hence,

$$ \|E^{(k)}\|_{\mathcal{M}}^2 = -\text{trace}(E^{(k)}R^{(k)}) = \text{trace}(\mathcal{M}(X)X) - \text{trace}(\mathcal{M}(X^{(k)})X^{(k)}). $$

297 The result follows from substituting into the inequality (4.1). □

298 This property has applications in different directions, for instance in approximating the $\mathcal{H}$-energy norm of a linear dynamical system, see, e.g., [1].

299 **5. Decay properties of the solution's singular values.** For $BB^T$ of very low rank, it is known that the singular values of the solution to the Lyapunov equation $AX + XA + BB^T = 0$ exhibit an exponentially decaying pattern, which depends on the spectral properties of $A$ and to a lesser extent on the right-hand side matrix $B$. In addition to being very helpful when questing for low-rank approximate solutions, the topic is fascinating and has attracted a lot of attention, especially in the case of non-normal $A$, see, e.g., [24],[3],[31],[2],[20],[27],[14].

300 Let $\sigma_1 \geq \sigma_2 \geq \dots \geq \sigma_n$ be the singular values of the solution matrix $X$, and recall that (Eckart-Young-Mirsky theorem)

$$ \frac{\sigma_k}{\sigma_1} = \min_{\text{rank}(Y)=k} \frac{\|X - Y\|}{\|X\|}. $$

301 As a consequence of this property, for any rank-$k$ approximation $X^{(k)}$ to $X$, any bound for the error $\|X - X^{(k)}\|$ can be used as an upper estimate for the decay of the $k$th singular value. This has been a common strategy in the literature cited above. 302 A very realistic estimate stemming from rational function approximations of $X$ was derived in [21],

$$ (5.1) \qquad \lambda_{pk+1}(X) \leq \frac{16\|B\|_2^2}{\lambda_{\min}(A)} \exp\left(-\frac{k\pi^2}{\log(8\kappa(A))}\right), \quad 1 \leq pk < n, $$

310 where $\lambda_{\min}(A)$, $\lambda_{\max}(A)$ are the extreme eigenvalues of the symmetric and positive definite matrix $A$, and $\kappa(A) = \frac{\lambda_{\max}(A)}{\lambda_{\min}(A)}$ is the condition number of $A$; a similar estimate 311 explicitly using elliptic functions was earlier established in [27, sec.2.1.2].
The question thus arises whether the decaying properties carry over to the general form (1.1). The answer is in the affirmative, despite the extra term $U\Phi(X)U^T$. In fact, the main role is played by the rank of $[B, U]$. This can be readily seen by writing (1.1) as $AX + XA = -BB^T - U\Phi(X)U^T$, with

$$BB^T + U\Phi(X)U^T = [B, U] \begin{bmatrix} I & \\ & \Phi(X) \end{bmatrix} [B, U]^T.$$

Hence, independently of the action of $\Phi$ on $X$, the solution $X$ has a worst possible decay completely described by $A$ and a full rank $[B, U]$. The decay can be faster in case $\Phi(X)$ is rank deficient. We refer to [32] for bounds that explicitly take into account the factor $[B, U]$, though their relevance can be better appreciated for $A$ non-symmetric. Interestingly, the formalization via $[B, U]$ was used in the proof of Theorem 4.1 in [4] to derive a closed form approximation to (1.1) as a sum of exponentials.

In the previous sections we gave a solution in closed form for some of the considered problems. This makes the decay analysis more precise, while avoiding explicit reference to $\Phi(X)$. Indeed, in all cases, we have

$$X = M + \sum_{i} \gamma_{i} P_{i},$$

where $M$ and all $P_i$ have similar spectral decay properties, so that something similar to (5.1) still holds, where only $[B, U]$ affects the worst case decay.

EXAMPLE 5.1. We consider the problem $AX + XA + uv^T Xvu^T + BB^T = 0$ with $A$ diagonal with diagonal entries logarithmically distributed in the interval $[10^{-10}, 1]$ and dimension $n = 1000$, $u$ a vector with normally distributed random components (Matlab seed `rng=1`), $v = e_1$, and $B = 1$. In Figure 1 we report the distribution of the first few tens singular values of the solutions $Y_1 = -\mathcal{L}^{-1}(BB^T)$, $Y_2 = -\mathcal{L}^{-1}(uu^T)$, together with that of $X = Y_1 + \sigma Y_2$, obtained with the procedure described in section 2.2. Moreover, the normalized estimate in (5.1) is also shown. The estimate is quite descriptive of the true decay for $p = 2$, since $[B, u]$ is full rank. The decay of the two solutions $Y_1, Y_2$ is very similar, and since $\sigma = O(1)$, $Y_1$ and $Y_2$ equally contribute to the spectral decay of $X$. ◇

<figure><img src="image_1.png" /><figcaption>FIG. 1. Example 5.1. Singular values of the solutions $Y_1, Y_2$ and $X = Y_1 + \sigma Y_2$, and the estimate in (5.1).</figcaption></figure>
**6. Generalizations and outlook.** Many of the results we have stated can be generalized to the case

$$ (6.1) \qquad A_1 X + X A_2 + \sum_{i=1}^{\ell} U_i \Phi_i(X) U_i^T + B_1 B_2^T = 0, $$

as long as the associated coefficient operator is symmetric and positive definite. This setting may have broader applications, such as in the discretization of time- and parameter-dependent partial differential equations and in the numerical treatment of problems dealing with uncertainty in the data; see, e.g., [30] and the numerical experiments in [23].

We have shown that writing the solution to (1.1) as a sum of terms, allows us to devise a simple indicator of the effect on the solution of the extra addend in the equation, that is not associated with the Lyapunov equation. Moreover, easily characterized positive definiteness and monotonicity properties of the approximate solutions can help calibrate the approximation phase, thus monitoring computational costs and dynamic rank truncations.

**Acknowledgments.** We thank Peter Benner for a discussion on the problem analyzed in [6].

The author is a member of the INdAM Research Group GNCS. Moreover, her work was partially supported by the European Union - NextGenerationEU under the National Recovery and Resilience Plan (PNRR) - Mission 4 Education and research - Component 2 From research to business - Investment 1.1 Notice Prin 2022 - DD N. 104 of 2/2/2022, entitled “Low-rank Structures and Numerical Methods in Matrix and Tensor Computations and their Application”, code 20227PCCKZ – CUP J53D23003620006.

## REFERENCES

[1] A. C. ANTOULAS, *Approximation of large-scale Dynamical Systems*, Advances in Design and Control, SIAM, Philadelphia, 2005.

[2] A. C. ANTOULAS, D. C. SORENSEN, AND Y. ZHOU, *On the decay rate of Hankel singular values and related issues*, Systems & Control Letters, 46 (2002), pp. 323–342.

[3] J. BAKER, M. EMBREE, AND J. SABINO, *Fast singular value decay for Lyapunov solutions with nonnormal coefficients*, SIAM Journal on Matrix Analysis and Applications, 36 (2015), p. 656–668.

[4] P. BENNER AND T. BREITEN, *Low Rank Methods for a Class of Generalized Lyapunov Equations and Related Issues*, Numerische Mathematik, 124 (2013), pp. 441–470.

[5] P. BENNER AND T. DAMM, *Lyapunov equations, energy functionals, and model order reduction of bilinear and stochastic systems*, SIAM J. Control Optim., 49 (2011), pp. 686–711.

[6] P. BENNER AND P. GOYAL, *Balanced truncation for quadratic-bilinear control systems*, Advances in Computational Mathematics, 50 (2024), p. 88.

[7] D. A. BINI, B. IANNAZZO, AND B. MEINI, *Numerical Solution of Algebraic Riccati Equations*, Society for Industrial and Applied Mathematics, 2011.

[8] I. BIOLI, D. KRESSNER, AND L. ROBOL, *Preconditioned Low-Rank Riemannian Optimization for Symmetric Positive Definite Linear Matrix Equations*, Tech. Rep. 2408.16416, ArXiv, 2024.

[9] R. CHOUDHARY AND K. AHUJA, *Stability analysis of bilinear iterative rational Krylov algorithm*, Lin. Alg. Appl., 538 (2018), pp. 56–88.

[10] T. DAMM, *Direct methods and ADI-preconditioned Krylov subspace methods for generalized Lyapunov equations*, Num. Lin. Alg. with Appl., 15 (2008), pp. 853–871. Special issue on Matrix equations.

[11] G. GOLUB AND C. F. VAN LOAN, *Matrix Computations*, The Johns Hopkins University Press, Baltimore, 4th ed., 2013.
383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

[12] I. V. GOSEA AND A. C. ANTOULAS, *Data-driven model order reduction of quadratic-bilinear systems*, Numerical Linear Algebra with Applications, 25 (2018), p. e2200.

[13] W. GRAY AND J. MESKO, *Energy functions and algebraic gramians for bilinear systems*, in IFAC Proceedings Volumes, vol. 31 of 17, Elsevier, 1998, pp. 101–106.

[14] L. GRUBISIĆ AND D. KRESSNER, *On the eigenvalue decay of solutions to operator Lyapunov equations*, Systems & Control Letters, 73 (2014), pp. 42–47.

[15] Y. HAO AND V. SIMONCINI, *The Sherman-Morrison-Woodbury formula for generalized linear matrix equations and applications*, Numerical Linear Algebra w/Appl., 28 (2021), p. e2384.

[16] C. HARTMANN, A. ZUEVA, AND B. SCHAEFER-BUNG, *Balanced model reduction of bilinear systems with applications to positive systems*, tech. rep., Institut fuer Mathematik, Freie Universitaet, Berlin, 2010.

[17] H. V. HENDERSON AND S. R. SEARLE, *On deriving the inverse of a sum of matrices*, SIAM Rev., 23 (1981), pp. 53–60.

[18] N. J. HIGHAM, *Functions of Matrices - Theory and Computation*, SIAM, Philadelphia, USA, 2008.

[19] R. A. HORN AND C. R. JOHNSON, *Matrix Analysis*, Cambridge University Press, Cambridge, II ed., 2013.

[20] N. KOMAROFF, *Simultaneous eigenvalue lower bounds for the Lyapunov matrix equation*, IEEE Trans. on Autom. Control, 33 (1988), pp. 126–128.

[21] D. KRESSNER AND C. TOBLER, *Krylov subspace methods for linear systems with tensor product structure*, SIAM J. Matrix Anal. Appl., 31 (2010), pp. 1688–1714.

[22] THE MATHWORKS, INC., *MATLAB 7, r2020b* ed., 2020.

[23] D. PALITTA, M. IANNACITO, AND V. SIMONCINI, *A subspace-conjugate gradient method for linear matrix equations*, Tech. Rep. 2501.02938, ArXiv, January 2025. Available from http://arxiv.org/abs/2501.02938.

[24] T. PENZL, *Eigenvalue decay bounds for solutions of Lyapunov equations: the symmetric case*, Systems and Control Letters, 40 (2000), pp. 139–144.

[25] M. PORCELLI AND V. SIMONCINI, *Numerical solution of a class of quasi-linear matrix equations*, Linear Algebra and its Applications, 664 (2023), pp. 349–368.

[26] A. QUARTERONI, R. SACCO, AND F. SALERI, *Matematica Numerica*, Springer, 2008.

[27] J. SABINO, *Solution of Large-Scale Lyapunov Equations via the Block Modified Smith Method*, PhD thesis, Rice University, 2006.

[28] S. D. SHANK, V. SIMONCINI, AND D. B. SZYLD, *Efficient low-rank solutions of generalized Lyapunov equations*, Numerische Mathematik, 134 (2016), pp. 327–342.

[29] J. SHERMAN AND W. J. MORRISON, *Adjustment of an inverse matrix corresponding to a change in one element of a given matrix*, Ann. Math. Statist., 21 (1950), pp. 124–127.

[30] V. SIMONCINI, *Computational methods for linear matrix equations*, SIAM Rev., 58 (2016), pp. 377–441.

[31] D. SORENSEN AND Y. ZHOU, *Bounds on eigenvalue decay rates and sensitivity of solutions to Lyapunov equations*, Tech. Rep. 02-07, Rice University, Houston, Texas, 2002.

[32] N. TRUHAR AND K. VESELIĆ, *Bounds on the trace of a solution to the Lyapunov equation with a general stable matrix*, Systems & Control Letters, 56 (2007), pp. 493–503.