Optimality properties of Galerkin and Petrov-Galerkin methods for linear matrix equations

TL;DR

This paper investigates the optimality properties of Galerkin and Petrov-Galerkin methods when applied to linear matrix equations, extending their known benefits from numerical analysis to this broader context.

Contribution

It demonstrates that Galerkin and Petrov-Galerkin methods retain their optimality features for linear matrix equations and introduces constrained minimization techniques within this framework.

Findings

Galerkin and Petrov-Galerkin methods preserve optimality in linear matrix equations

Novel considerations for applying these methods to general matrix equations

Proposal of constrained minimization techniques in Petrov-Galerkin methods

Abstract

Galerkin and Petrov-Galerkin methods are some of the most successful solution procedures in numerical analysis. Their popularity is mainly due to the optimality properties of their approximate solution. We show that these features carry over to the (Petrov-)Galerkin methods applied for the solution of linear matrix equations. Some novel considerations about the use of Galerkin and Petrov-Galerkin schemes in the numerical treatment of general linear matrix equations are expounded and the use of constrained minimization techniques in the Petrov-Galerkin framework is proposed.

Equations

$$\mathcal{M}x=f,$$

$$A_{1}X{B}_1+A_{2}X{B}_2+\dots +A_{\ell }X{B}_{\ell }=F,$$

$$AX+X{A}^T=F,$$

$$Y=AXB\iff y=\mathcal{T}x,x=\mathrm{vec}(X),y=\mathrm{vec}(Y),$$

$$V_m^T{r}_m=0\iff V_m^T\mathcal{M}{V}_m{y}_m=V_m^{T}f.$$

$$\Vert e_m{\Vert }_{\mathcal{M}}^2=\Vert {\mathcal{M}}^{1/2}(x-x_m){\Vert }^2=\Vert {\mathcal{M}}^{1/2}x-{\mathcal{M}}^{1/2}{V}_m{y}_m{\Vert }^2.$$

$$({\mathcal{M}}^{1/2}{V}_m{)}^T{\mathcal{M}}^{1/2}{V}_m{y}_m=({\mathcal{M}}^{1/2}{V}_m{)}^T{\mathcal{M}}^{1/2}x,$$

$$V_k^T{R}_k{V}_k=0\iff (V_k\otimes V_k{)}^T{r}_k=0,$$

$$V_k^T{R}_k{W}_k=0\iff (W_k\otimes V_k{)}^T{r}_k=0,$$

$$\begin{array}[]{lrll}{\cal S}:&\mathbb{R}^{n\times p}&\rightarrow&\mathbb{R}^{n\times p}\\ &X&\mapsto&\displaystyle\sum_{j=1}^{\ell}A_{j}XB_{j},\\ \end{array}$$

$$x^T{\mathcal{S}}_{\ell }x=\mathrm{trace}\left(\sum _{j=1}^{\ell }{X}^T{A}_{j}X{B}_j\right).$$

$$\Vert X-X_k{\Vert }_{\mathcal{S}}=\underset{\genfrac{}{}{0ex}{}{Z=V_{k}Y{W}_k^T}{Y\in {\mathbb{R}}^{k\times k}}}{\min }\Vert X-Z{\Vert }_{\mathcal{S}}.$$

$$\Vert e_k{\Vert }_{\mathcal{S}}^2=e_k^T{\mathcal{S}}_{\ell }{e}_k=\mathrm{trace}((X-X_k{)}^T\mathcal{S}(X-X_k))=\Vert X-X_k{\Vert }_{\mathcal{S}}^2,$$

$$\Vert X{\Vert }_{\mathcal{S}}^2=2\mathrm{trace}(X^{T}AX).$$

$$\Vert E_k{\Vert }_{\mathcal{S}}^2=\underset{\genfrac{}{}{0ex}{}{Z=V_{k}Y{W}_k^T}{Y\in {\mathbb{R}}^{k\times k}}}{\min }\Vert X-Z{\Vert }_{\mathcal{S}}^2=2\mathrm{trace}(E_k^{T}A{E}_k).$$

$$\Vert X-X_k{\Vert }_2⩽2{\int }_0^{\infty }{e}^{-t{\alpha }_{\min }(A)}\Vert x-x_m{\Vert }_{2}dt,{\alpha }_{\min }(A)={\lambda }_{\min }((A+A^T)/2),$$

$$\Vert X-X_k{\Vert }_2$$

$$AX+XB=b_1{b}_2^T,$$

$$X={\int }_0^{\infty }{e}^{-tA}{b}_1{b}_2^T{e}^{-tB}dt.$$

$$V_k^T{R}_k{W}_k=0\iff (V_k^{T}A{V}_k)Y_k+Y_k(W_k^{T}B{W}_k)-(V_k^T{b}_1)(b_2^T{W}_k)=0.$$

$$X_k=V_k{\int }_0^{\infty }{e}^{-t{A}_k}(V_k^T{b}_1)(b_2^T{W}_k)e^{-t{B}_k}dt{W}_k^T.$$

$$\Vert X-X_k{\Vert }_{\star }$$

$${\hat{\kappa }}_A=\frac{{\lambda }_{\max }(A)+{\lambda }_{\min }(B)}{{\lambda }_{\min }(A)+{\lambda }_{\min }(B)},{\hat{\kappa }}_B=\frac{{\lambda }_{\max }(B)+{\lambda }_{\min }(A)}{{\lambda }_{\min }(B)+{\lambda }_{\min }(A)},$$

$$\Vert X-X_k{\Vert }_{\star }⩽\frac{2}{{\lambda }_{\min }(A)+{\lambda }_{\min }(B)}\left(\frac{\sqrt{{\hat{\kappa }}_A}+1}{\sqrt{{\hat{\kappa }}_A}}{\left(\frac{\sqrt{{\hat{\kappa }}_A}-1}{\sqrt{{\hat{\kappa }}_A}+1}\right)}^k+\frac{\sqrt{{\hat{\kappa }}_B}+1}{\sqrt{{\hat{\kappa }}_B}}{\left(\frac{\sqrt{{\hat{\kappa }}_B}-1}{\sqrt{{\hat{\kappa }}_B}+1}\right)}^k\right),$$

$$\sum _{j,k=0}^{\ell }{\alpha }_{j,k}{A}^{j}X{B}^k=F,$$

$$\left(\sum _{j=1}^{\ell }{B}_j^T\otimes A_j\right)\text{vec}(X)=\text{vec}(F).$$

$$\mathcal{A}\text{vec}(X)=\text{vec}(b{b}^T),\mathcal{A}=A\otimes I+I\otimes A\in {\mathbb{R}}^{n^2\times n^2}.$$

$$\frac{\Vert x-x_k^{cg}{\Vert }_{\mathcal{A}}}{\Vert x{\Vert }_{\mathcal{A}}}⩽2{\left(\frac{\sqrt{\kappa }-1}{\sqrt{\kappa }+1}\right)}^k,$$

$$\sqrt{2{\lambda }_{\min }(A)}\Vert v{\Vert }_2⩽\Vert v{\Vert }_{\mathcal{A}}⩽\sqrt{2{\lambda }_{\max }(A)}\Vert v{\Vert }_2.$$

$$\frac{\Vert X-X_k^{cg}{\Vert }_F}{\Vert X{\Vert }_F}=\frac{\Vert x-x_k^{cg}{\Vert }_2}{\Vert x{\Vert }_2}⩽\sqrt{\kappa }\frac{\Vert x-x_k^{cg}{\Vert }_{\mathcal{A}}}{\Vert x{\Vert }_{\mathcal{A}}}⩽2\sqrt{\kappa }{\left(\frac{\sqrt{\kappa }-1}{\sqrt{\kappa }+1}\right)}^k.$$

Full text

Some novel considerations about the use of Galerkin and Petrov-Galerkin schemes in the numerical treatment of general linear matrix equations are expounded and the use of constrained minimization techniques in the Petrov-Galerkin framework is proposed.

Dedicated to Volker Mehrmann on the occasion of his 65th birthday

{AMS}

65F10, 65F30, 15A06

1 Introduction

Many state-of-the-art solution procedures for algebraic linear systems of the form

[TABLE]

where ${\cal M}\in{\mathbb{R}}^{N\times N}$ and $f\in{\mathbb{R}}^{N}$, are based on projection. Given a subspace ${\cal K}_{m}$ of dimension $m$, and a matrix ${\cal V}_{m}$ whose orthonormal columns span ${\cal K}_{m}$, these methods seek an approximate solution $x_{m}={\cal V}_{m}y_{m}$ for some $y_{m}\in{\mathbb{R}}^{m}$ by imposing certain conditions. The most successful projection procedures impose either a Galerkin or a Petrov-Galerkin condition on the residual $r_{m}=f-{\cal M}x_{m}$. See, e.g., [40]. These conditions are very general, and they are at the basis of many approximation methods, beyond the algebraic context of interest here; any approximation strategy associated with an inner product can determine the projected solution by one of such a condition. Finite element methods, both at the continuous and discrete levels, strongly rely on this methodology; see, e.g., [49], but also eigenvalue problems [39].

It is very important to realize that this is a methodology, not a single method: the approximation space can be generated independently of the condition, and in a way to make the computation of $y_{m}$ more effective, while obtaining a sufficiently accurate approximation with the smallest possible space dimension.

A fundamental property of the Galerkin methodology is obtained whenever the coefficient matrix ${\cal M}$ is symmetric and positive definite (spd): the Galerkin condition on the residual corresponds to minimizing the error vector in the norm associated with ${\cal M}$ over the approximation space. This property is at the basis of the convergence analysis of methods such as the Conjugate Gradient (CG) [21], and it ensures monotonic convergence, in addition to finite termination, in exact precision arithmetic.

When $\mathcal{M}$ is not spd, the application of the Galerkin method does not automatically imply a minimization of the error norm. Nevertheless, a certain family of Petrov-Galerkin procedures still fulfills an optimality property. Indeed, these methods minimize the residual norm over the space $\mathcal{MK}_{m}$. See, e.g., [40]. Some of the most popular solvers for linear systems such as MINRES [34] and GMRES [41] belong to this collection of methods.

In the past decades, projection techniques have been successfully used to solve linear matrix equations of the form

[TABLE]

that have arisen as a natural algebraic model for discretized partial differential equations (PDEs), possibly including stochastic terms or parameter dependent coefficient matrices [5, 12, 38, 35], for PDE-constrained optimization problems [48], data assimilation [19], and many other applied contexts, including building blocks of other numerical procedures [30]; see also [44, 14] for further references.

The general matrix equation (2) covers two well known cases, the (generalized) Sylvester equation (for $\ell=2$), and the Lyapunov equation

[TABLE]

which plays a crucial role in many applications such as control and system theory [13, 2], and in the solution of Riccati equations by the Newton method, in which a Lyapunov equation needs to be solved at each Newton step. See, e.g., [33].

The aim of this paper is to generalize the optimality properties of the Galerkin and Petrov-Galerkin methods to matrix equations, and to extend other convergence properties of CG and some related schemes to the matrix setting. Some of the proposed results are new, some others can be found in articles scattered in the literature in different contexts. We thus provide a more uniform presentation of these results.

To introduce a matrix version of the error and residual minimization, we first recall the relation between matrix-matrix operations and Kronecker products. Indeed, if $\otimes$ denotes the Kronecker product and ${\cal T}:=B^{T}\otimes A$, then

[TABLE]

where the usual “vec($\cdot$)” operator stacks the columns of the argument matrix one after the other into a long vector.

2 The Galerkin condition

In this section we first recall the result connecting the Galerkin condition on the residual with the minimization of the error norm when this is applied to the solution of linear systems, and then we show that similar results can be obtained also in the matrix equation setting. For the rest of the section we assume that ${\cal M}$ in (1) is symmetric and positive definite.

2.1 The linear system setting

Let $x_{m}=V_{m}y_{m}$ be an approximation to the true solution of (1), and let $e_{m}=x-x_{m}$, $r_{m}=f-{\cal M}x_{m}$ be the associated error and residual, respectively. We recall that imposing the Galerkin condition yields

[TABLE]

Note that the coefficient matrix $V_{m}^{T}{\cal M}V_{m}$ is symmetric and positive definite. Solving this system yields the “projected” vector $y_{m}$, so as to completely define $x_{m}$.

Let $\|e_{m}\|_{\cal M}^{2}:=e_{m}^{T}{\cal M}e_{m}$ be the ${\cal M}$-norm associated with the spd matrix ${\cal M}$. For the error we thus have

[TABLE]

The minimization of the error ${\cal M}$-norm thus corresponds to solving the least squares problem on the right, which gives

[TABLE]

which, upon simplifications of the transpositions yields $V_{m}^{T}{\cal M}V_{m}y_{m}=V_{m}^{T}f$, that is, using (4), $V_{m}^{T}r_{m}=0$.

2.2 Galerkin method and error minimization for matrix equations

To simplify the presentation, we first discuss Galerkin projection with the Lyapunov equation. Given the equation (3) with $A$ spd and $F=F^{T}$, then it can be shown that $X$ is symmetric. Letting ${\rm range}(V_{k})$ be an approximation space, we can determine an approximation to $X$ as $X_{k}=V_{k}Y_{k}V_{k}^{T}$, which in vector notation is written as $\text{vec}(X_{k})=(V_{k}\otimes V_{k}){\rm vec}(Y_{k})$. The matrix $Y_{k}$ is obtained by imposing the Galerkin condition in a matrix sense to the residual matrix $R_{k}=F-(AX_{k}+X_{k}A)$, that is

[TABLE]

where $r_{k}={\rm vec}(R_{k})$. Therefore, if one writes the Lyapunov equation by means of the Kronecker formulation, the obtained approximation space is ${\cal K}_{m}={\rm range}(V_{k}\otimes V_{k})$.

We explicitly notice that $X_{k}$ belongs to ${\rm range}(V_{k})$, which is much smaller than range$(V_{k}\otimes V_{k})$. Therefore, by sticking to the matrix equation formulation, we expect to build a much smaller approximation space than if a blind use of the Kronecker form were used. In other words, by exploiting the original matrix structure, no redundant information is sought after. In section section 4 we provide a rigorous analysis of this argument. See also [27]. To be able to exploit the derivation in (5) we will define an error matrix and the associated inner product.

The generalization to the multiterm linear equation (2) requires the definition of two approximation spaces, since the right and left coefficient matrices are not necessarily the same. Therefore, let range$(V_{k})$ and range$(W_{k})$ be two approximation spaces of dimension $k$ each111In principle, we can have $\text{dim}(\text{range}(V_{k}))\neq\text{dim}(\text{range}(W_{k}))$. Here we consider $\text{dim}(\text{range}(V_{k}))=\text{dim}(\text{range}(W_{k}))=k$ for the sake of simplicity in the presentation., and let us write the approximation to $X$ as $X_{k}=V_{k}Y_{k}W_{k}^{T}$. With the residual matrix $R_{k}=F-\sum_{j=1}^{\ell}A_{j}X_{k}B_{j}$, the Galerkin condition now takes the form

[TABLE]

where $r_{k}={\rm vec}(R_{k})$, so that ${\cal K}_{m}={\rm range}(W_{k}\otimes V_{k})$ with $m=k^{2}$ in the Kronecker formulation.

To adapt the error minimization procedure to the matrix equation setting we first introduce a matrix norm, that allows us to make a connection with the ${\cal M}$-norm of the error vector. A corresponding derivation for $\ell=2$ can be found, for instance, in [52, p. 2557] and [11, p. 149].

Definition 2.1**.**

Let

[TABLE]

and ${\cal S}_{\ell}=\sum_{j=1}^{\ell}B_{j}^{T}\otimes A_{j}$. We say that the operator ${\cal S}$ is symmetric and positive definite if for any $0\neq x\in{\mathbb{R}}^{np}$, $x={\rm vec}(X)$, with $X\in{\mathbb{R}}^{n\times p}$, it holds that ${\cal S}_{\ell}={\cal S}_{\ell}^{T}$ and $x^{T}{\cal S}_{\ell}x>0$, where

[TABLE]

*The norm induced by this operator will be denoted by $\|X\|_{\cal S}$. *

Note that any linear operator $\mathcal{L}:\mathbb{R}^{n\times p}\rightarrow\mathbb{R}^{n\times p}$ can be written in the form (6) with a uniquely defined minimum number of terms $\ell$ called the Sylvester index. See [25].

Assuming $\mathcal{S}$ to be spd, in the following proposition we show that the error matrix is minimized in the ${\cal S}$-norm.

Proposition 2.2**.**

Let ${\cal S}(X)=F$ with ${\cal S}:X\mapsto\sum_{j}A_{j}XB_{j}$ spd, and let $\text{range}(V_{k})$, $\text{range}(W_{k})$ be the constructed approximation spaces, so that $X_{k}=V_{k}Y_{k}W_{k}^{T}$ is the Galerkin approximate solution. Then

[TABLE]

Proof 2.3**.**

Let $e_{k}={\rm vec}(X-X_{k})$ be the error vector, $r_{k}={\rm vec}(F-\sum_{j}A_{j}X_{k}B_{j})$ the residual vector, ${\cal K}_{m}=\text{range}(W_{k}\otimes V_{k})$ the approximation space and ${\cal S}_{\ell}=\sum_{j}B_{j}^{T}\otimes A_{j}$ the coefficient matrix. Then, since $\mathcal{S}$ is spd by assumption, also ${\cal S}_{\ell}$ is spd, and the Galerkin condition ${\cal V}_{m}^{T}r_{k}=0$, $\mathcal{V}_{m}=W_{k}\otimes V_{k}$, corresponds to the minimization of the error. More precisely, it holds

[TABLE]

*and the proof is completed. *

Proposition 2.2 states that as long as the approximation spaces are expanded, the error will decrease monotonically in the considered norm. A Galerkin approach for a multiterm linear matrix equation was for instance employed in [38]; the proposition above thus ensures that under the stated hypotheses on the data the method will minimize the error as the approximation spaces grow. See also Example 2.5.

A result similar to the one stated in Proposition 2.2 can be found in [27] where the authors consider specific approximation spaces and assume $\mathcal{S}$ to be a so-called Laplace-like operator. Proposition 2.2 shows the strength of the Galerkin method, also in the general matrix equation setting. Indeed, the optimality condition of the Galerkin method does neither depend on the adopted approximation spaces nor on the definition of $\mathcal{S}$, as long as this is spd.

Given a general linear matrix equation (2) written in the form $\mathcal{S}(X)=F$, one would like to characterize the symmetry and positive definiteness of $\mathcal{S}$ by looking only at the properties of the matrices $A_{j}$ and $B_{j}$ and avoid the construction of the large matrix $\mathcal{S}_{\ell}$.

Assuming $\ell$ to be the Sylvester index of $\mathcal{S}$, it is easy to show that $\mathcal{S}$ is a symmetric operator if and only if the matrices $A_{j}$ and $B_{j}$ are symmetric for all $j=1,\ldots,\ell$, whereas, in general, it is not possible to identify the positive definiteness of $\mathcal{S}$ by examining the spectral distributions of $A_{j}$ and $B_{j}$, even when these are completely known. See, e.g., [28]. Note that for $\mathcal{S}$ to be spd it is not necessary for all the $A_{j}$’s and $B_{j}$’s to be positive definite. Nevertheless, if $A_{j}$, $B_{j}$ are positive definite for all $j=1,\ldots,\ell$, then $\mathcal{S}$ is positive definite; see, e.g., [52, Proposition 3.1] for $\ell=2$. Therefore, in the case of the Lyapunov equation with $A$ spd, also the operator $\mathcal{S}$ is spd and it holds that

[TABLE]

Another case where the properties of $\mathcal{S}$ can be determined in terms of the (symmetric) coefficient matrices $A_{j}$ and $B_{j}$ is the Sylvester operator ${\cal S}:X\mapsto AX+XB$. By exploiting the property of the Kronecker product, it holds that $\mathcal{S}$ is positive definite if and only if $\lambda_{i}(A)+\lambda_{j}(B)>0$ for all $i$s and $j$s. Moreover, the norm $\|\cdot\|_{\cal S}$ can be written as $\|X\|_{\cal S}^{2}={\rm trace}(X^{T}AX)+{\rm trace}(XBX^{T})$.

Remark 2.4**.**

Consider the Lyapunov equation (3) with the spd coefficient matrix $A$, and let $E_{k}:=X-X_{k}$ be the corresponding error matrix. Then, the previous discussion shows that

[TABLE]

In the remark above we have not specified whether the known term $F$ in (3) needs to be symmetric. If $F$ is symmetric, then indeed the two spaces can coincide, and $E_{k}$ is also symmetric. On the other hand, if $F$ has the form $F=F_{1}F_{2}^{T}$, possibly low rank, natural choices as approximation spaces are such that ${\rm range}(F_{1})\subseteq{\rm range}(V_{k})$ and ${\rm range}(F_{2})\subseteq{\rm range}(W_{k})$, so that the (vector) residual is orthogonal to $F_{2}\otimes F_{1}$. A possible alternative could use $V_{k}=W_{k}$ such that ${\rm range}(F_{1}),{\rm range}(F_{2})\subseteq{\rm range}(V_{k})$, where however in general we expect ${\rm range}(V_{k})$ to have larger dimension than in the previous case.

Example 2.5**.**

By applying the stochastic Galerkin methodology for the discretization of elliptic stochastic PDEs [3], the resulting algebraic formulation can be written as the linear matrix equation (2) with typically $\ell>2$. When dealing with the stochastic steady-state diffusion problem with homogeneous Dirichlet boundary conditions, the symmetric matrices $A_{j}$ and $B_{j}$ may not all be positive definite; nonetheless, the associated operator $\mathcal{S}$ is symmetric and indeed positive definite (see, e.g., [37]), so that the previous theory applies. In the following we consider the Galerkin approach developed in [38] – based on the rational Krylov subspace – to illustrate the monotonic decrease of the error $\mathcal{S}$-norm as the approximation space increases222The Matlab code is available at http://www.dm.unibo.it/~simoncin/software.html..We generate $A_{j}$ and $B_{j}$ as the second test case in the S-IFISS package [43] with the default setting for all the requested parameters. This yields* a linear matrix equation of the form (2) with $\ell=6$, $A_{j}\in\mathbb{R}^{n\times n}$, $n=225$, and $B_{j}\in\mathbb{R}^{p\times p}$, $p=56$. The right-hand side $F$ has rank 1. Thanks to the small problem dimension, we could compute the vectorized solution $x\in\mathbb{R}^{np}$ as $x={\rm vec}(X)={\cal S}_{\ell}^{-1}f$ (Matlab function “$\setminus$”), to be used as a reference “exact” solution. In particular, if $X_{k}$ denotes the approximate solution obtained after $k$ iterations of the Galerkin method, we compute $\|X-X_{k}\|_{\mathcal{S}}/\|X\|_{\mathcal{S}}$ until this falls below $10^{-6}$. Figure 1 displays the history of this relative error $\mathcal{S}$-norm, illustrating the expected monotonically non-increasing curve.*

3 Convergence properties

In the previous section we have shown that the Galerkin condition leads to a minimization of the error $\mathcal{S}$-norm and this property does not depend on the selected space $\mathcal{K}_{k}=\text{range}(V_{k})$. In actual computations, a measurable estimate of the error is needed and in [45] an upper bound on the Euclidean norm of the error is provided in the case of the Lyapunov equation (3) with rank-one right-hand side $F=bb^{T}$ with $\|b\|=1$, and a positive definite but not necessarily symmetric $A$. By exploiting the closed-form of the solution $X$, the authors showed that

[TABLE]

where $x=e^{-tA}b$, $x_{k}=V_{k}e^{-tA_{k}}e_{1}$, $A_{k}:=V_{k}^{T}AV_{k}$, and $\|\cdot\|_{2}$ denotes the Euclidean norm.

This led to the following proposition when the selected approximation space is the Krylov subspace $\text{range}(V_{k})=K_{k}(A,b)={\rm span}\{b,Ab,\ldots,A^{k-1}b\}$ and $A$ is symmetric.

Proposition 3.1** ([45]).**

Let $A$ be spd, and let $\lambda_{\max}$ and $\lambda_{\min}$ be the largest and the smallest eigenvalue of $A$, respectively. Denoting by $\hat{\kappa}=(\lambda_{\max}+\lambda_{\min})/2\lambda_{\min}$ the condition number of the spd matrix $A+\lambda_{\min}I$, then the Galerkin approximate solution $X_{k}=V_{k}Y_{k}V_{k}^{T}$ satisfies

[TABLE]

This bound, in terms of slope as $k$ increases, was shown to be sharp in [45]. Notice that the bound (7) holds also for the Frobenius norm of the error, namely $\|X-X_{k}\|_{F}$. Indeed, we can still write $\|X-X_{k}\|_{F}\leqslant 2\int_{0}^{\infty}e^{-t\alpha_{\min}(A)}\|x-x_{m}\|_{F}dt$ and the rest of the proof of Proposition 3.1 makes use of bounds for norms of vectors only, for which the Euclidean and the Frobenius norms coincide. See [45, Proposition 3.1] for more details. The bound can be generalized to the use of other spaces, such as rational Krylov subspaces, see, e.g., [6, 7, 24, 18].

We generalize the bound presented in Proposition 3.1 to the case of the Sylvester equation,

[TABLE]

with $A$ and $B$ symmetric and positive definite; without loss of generality, we can assume that $\|b_{1}\|_{\star}=\|b_{2}\|_{\star}=1$ where $\|\cdot\|_{\star}$ denotes either the Euclidean or the Frobenius norms.

We first recall the Cauchy representation of the solution matrix $X$ to (8). Let us for now only assume that $A$ and $B$ are positive definite, and not necessarily symmetric. We can write (see, e.g., [28])

[TABLE]

Consider the approximation $X_{k}=V_{k}Y_{k}W_{k}^{T}$ where $V_{k}$ and $W_{k}$ span suitable subspaces and both have orthonormal columns. The matrix $Y_{k}$ is obtained by imposing the Galerkin condition on $R_{k}=AX_{k}+X_{k}B-b_{1}b_{2}^{T}$, that is

[TABLE]

Let $A_{k}:=V_{k}^{T}AV_{k},\;B_{k}:=W_{k}^{T}BW_{k}$. Thus $Y_{k}$ is obtained by solving a reduced Sylvester equation, whose size depends on the approximation space dimensions. Since the spectrum of $A_{k}$ ($B_{k}$) is contained in the spectral region of $A$ ($B$) , we have that $\Lambda(A_{k})+\Lambda(B_{k})\subset\mathbb{C}_{+}$ and the matrix $Y_{k}$ can be written in integral form as $Y_{k}=\int_{0}^{\infty}e^{-tA_{k}}(V_{k}^{T}b_{1})(b_{2}^{T}W_{k})e^{-tB_{k}}dt$ so that

[TABLE]

Let $x:=e^{-tA}b_{1},x_{k}:=V_{k}e^{-tA_{k}}(V_{k}^{T}b_{1})$, $y:=e^{-tB}b_{2},y_{k}:=W_{k}e^{-tB_{k}}(W_{k}^{T}b_{2})$. Then, using $\|x\|_{\star}\leqslant e^{-t\alpha_{\min}(A)}$ (see, e.g., [16, Lemma 3.2.1]), and since $\alpha_{\min}(A_{k})\geqslant\alpha_{\min}(A)$, it holds that $\|x_{k}\|_{\star}\leqslant e^{-t\alpha_{\min}(A)}$. Similarly, $\|y\|_{\star},\|y_{k}\|_{\star}\leqslant e^{-t\alpha_{\min}(B)}$. Therefore, (see also [27, Lemma 4.7])

[TABLE]

where $\hat{y}=e^{-t(B+\lambda_{\min}(A)I)}b_{2}$, $\hat{x}=e^{-t(A+\lambda_{\min}(B)I)}b_{1}$, and analogously for $\hat{y}_{k},\hat{x}_{k}$. The inequality in (9) states that the $\star$-norm of the error associated with the Galerkin solution can be bounded by integrating over $[0,\infty)$ the errors obtained in the approximation of the exponential of the shifted matrices $B+\lambda_{\min}(A)I$ and $A+\lambda_{\min}(B)I$.

In the next proposition we specialize the bound above when the Krylov subspaces $\text{range}(V_{k})=K_{k}(A,b_{1})$ and $\text{range}(W_{k})=K_{k}(B,b_{2})$ are adopted as approximation spaces and $A$, $B$ are both symmetric and positive definite. To this end, let us define $\lambda_{\min}(A)$, $\lambda_{\max}(A)$, $\lambda_{\min}(B)$, and $\lambda_{\max}(B)$ to be the extreme eigenvalues of $A$ and $B$, respectively, and

[TABLE]

the condition numbers of $A+\lambda_{\min}(B)I$ and $B+\lambda_{\min}(A)I$, respectively.

Proposition 3.2**.**

Let $A$ and $B$ be spd and $\text{range}(V_{k})=K_{k}(A,b_{1})$, $\text{range}(W_{k})=K_{k}(B,b_{2})$. Then the Galerkin approximate solution $X_{k}=V_{k}Y_{k}W_{k}^{T}$ to (8) is such that

[TABLE]

*where $\|\cdot\|_{\star}$ denotes either the Euclidean or the Frobenius norm. *

Proof 3.3**.**

*The proof can be obtained by applying the same arguments of the proof of [45, Proposition 3.1] to the single integrals $\int_{0}^{\infty}\|\hat{y}-\hat{y}_{k}\|_{\star}dt$, $\int_{0}^{\infty}\|\hat{x}-\hat{x}_{k}\|_{\star}dt$ in (9). *

Convergence results for generic matrix equations of the form (2) are difficult to derive as no easy-to-handle closed-form solution is known in general. The main difficulty is given by the fact that the exponential of a Kronecker sum $\sum_{j=1}^{\ell}B_{j}^{T}\otimes A_{j}$ cannot be separated in the product of the exponentials of the single terms if no further assumptions on $A_{j}$ and $B_{j}$ are considered.

By adapting the reasonings proposed in this section, one may be able to deduct error estimates for some special equations of the form

[TABLE]

where the coefficient matrices are given as powers of two seed matrices $A$ and $B$, and $\alpha_{j,k}\in\mathbb{R}$ for all $j,k$. Indeed, in this case, the exact solution $X$ can be written in integral form as illustrated in [28, Theorem 4]. However, such derivations deserve a separate analysis.

4 Comparison with the Kronecker formulation

Given a linear matrix equation of the form (2), the simplest-minded numerical procedure for its solution consists in applying well-established iterative schemes to the vector linear system obtained from (2) by Kronecker transformations, namely

[TABLE]

Sometimes this is the only option as effective algorithms to solve (2) in its natural matrix equation form are still lacking in the literature in the most general case. The methods developed so far require some additional assumptions on the coefficient matrices $A_{j}$, $B_{j}$; see, e.g., [10, 42, 23, 26, 38].

In this section we show that exploiting the matrix structure of equation (2) not only leads to numerical algorithms with lower computational costs per iteration and modest storage demands, but they also avoid some spectral redundancy encoded in the problem formulation (10). Such a redundancy often leads to a delay in the convergence of the adopted solution scheme when iterative procedures are applied to (10). A similar discussion can be found in [27, Remark 4.5] for more general tensor structured problems.

To illustrate this phenomenon we consider a Lyapunov equation of the form (3) with $A\in\mathbb{R}^{n\times n}$ spd and $F=bb^{T}$, $b\in\mathbb{R}^{n}$, $\|b\|=1$. We compare the Galerkin method applied to the matrix equation (3) with the CG method applied to the linear system

[TABLE]

Notice that since $A$ is spd, $\mathcal{A}$ is also spd. Let $x={\rm vec}(X)$ be the exact solution to (11). Let the CG initial guess be equal to the zero vector, and let $x_{k}^{cg}$ be the approximate solution to $x$ obtained after $k$ CG iterations. Then the following classical bound for the energy-norm of the error $x-x_{k}^{cg}$ holds

[TABLE]

where $\kappa=\lambda_{\max}(\mathcal{A})/\lambda_{\min}(\mathcal{A})=\lambda_{\max}(A)/\lambda_{\min}(A)$. See, e.g., [20, Theorem 10.2.6]. This bound may be rather pessimistic since it takes into account neither the role of the right-hand side nor the actual spectral distribution of $\mathcal{A}$. See, e.g., [51, 9, 8, 29].

We want to compare the bound in (12) with the estimate proposed in Proposition 3.1, using the same norms and relative quantities. To this end, we recall that for any vector $v$ it holds that

[TABLE]

In particular, letting $X_{k}^{cg}\in\mathbb{R}^{n\times n}$ be such that $\text{vec}(X_{k}^{cg})=x_{k}^{cg}$, we have

[TABLE]

Therefore, to obtain a relative error (in Frobenius norm) of less than $\varepsilon$, a sufficient number $k_{*}^{(cg)}$ of CG iterations is given by

[TABLE]

If $X_{k}$ denotes the approximate solution computed after $k$ iterations of the Galerkin-based method with $K_{k}(A,b)$ as approximation space, the error norm bound in (7) can be written in relative terms as

[TABLE]

where we used $\|x\|_{2}\geqslant\lambda_{\min}(\mathcal{A}^{-1})\,\|\text{vec}(bb^{T})\|_{F}=1/(\lambda_{\max}+\lambda_{\min})$. Once again, to obtain a relative error (in Frobenius norm) of less than $\varepsilon$, a sufficient number $k_{*}^{(G)}$ of iterations is given by

[TABLE]

The bounds (13)–(14) show that the asymptotic behavior of the relative error norms of CG and the Galerkin method are guided by $\kappa$ and $\hat{\kappa}$, respectively, where $\hat{\kappa}$ is always smaller than $\kappa$, for $\kappa>1$. Indeed,

[TABLE]

The worse conditioning of the linear system formulation (11) may lead to a delay in the convergence of CG so that, for a fixed threshold, CG may require more iterations to converge than the Galerkin method applied to the matrix equation (3). This is numerically illustrated in the examples below.

We once again stress that the similarities of the two formulations (matrix equation and Kronecker form) highlight the fact that what makes the matrix equation context more efficient than CG on ${\cal A}x=b$ is the special choice of the approximation space, that is ${\cal K}_{m}=\text{range}(V_{k}\otimes V_{k})$, which heavily takes into account the Kronecker sum structure of ${\cal A}$. On the other hand, CG applied blindly on ${\cal A}$ generates a redundant approximation space.

Example 4.1**.**

We consider the spd matrix $A=QDQ^{T}\in\mathbb{R}^{n\times n}$, where $D$ is a diagonal matrix whose diagonal entries are uniformly distributed (in logarithmic scale) values between 1 and 100, and $Q$ is orthogonal. This means that $\kappa=100$ and $\hat{\kappa}=50.5$ for any $n$. The vector $b\in\mathbb{R}^{n}$ is a random vector with unit norm.

For $\varepsilon=10^{-6}$, a direct computation shows that $k_{*}^{(G)}=68$ iterations of the Galerkin method are sufficient to get $\|X-X_{k_{*}^{(G)}}\|_{F}/\|X\|_{F}\leqslant\varepsilon$, whereas according to the bound (13), $k_{*}^{(cg)}=84$ iterations are required for CG to reach the same accuracy when solving ${\cal A}x=b$. In practice, the number of actual iterations can be lower, since this estimate is obtained from a bound.

Figure 2 reports the error convergence history of the two iterations, using logarithmic scale for $n=1000$. The two methods are stopped as soon as the relative error norm becomes smaller than $\varepsilon$. The “exact” solution $X$ was computed with the Bartels-Stewart method [4], which was feasible due to the small problem size.

Both methods require slightly fewer iterations than predicted by the bounds. Nonetheless, we can still appreciate that CG applied to the linear system (11) requires more iterations than the Galerkin method applied to the matrix equation (3) to achieve the same prescribed accuracy. * *

Example 4.2**.**

We modify the data of Example 4.1 by replacing $\lambda_{\min}(A)=1$ with $\lambda_{\min}(A)=\lambda_{1}=0.001$, while the other eigenvalues are such that $\lambda_{2}=2,\ldots,\lambda_{n}=n$. Here $b\in\mathbb{R}^{n}$ is the vector of all ones normalized. The relative error energy norm obtained by CG and the Galerkin method is reported in the left plot of Figure 3 for $n=100$. Notice that with such a $n$ we have $\kappa=10^{5},\hat{\kappa}\approx 5\cdot 10^{4}$.

Both methods stagnate in the initial phase of the solution process, followed by a rapid convergence afterwards. The stagnation phase is significantly longer for CG, contributing to the overall CG delay. A closer look at the convergence history of the Ritz values towards the eigenvalues of the corresponding coefficient matrices provides a better understanding. In this setting, the Ritz values for the Galerkin and CG methods are the eigenvalues of $V_{k}^{T}AV_{k}$ and of $Q_{k}^{T}{\cal A}Q_{k}$, respectively, where the columns of $Q_{k}$ are the orthonormal basis of the space generated by CG (here we used the Arnoldi procedure to compute $Q_{k}$). Recalling that $\lambda_{\min}(A)=0.001$, $\lambda_{\max}(A)=100$ so that $\lambda_{\min}(\mathcal{A})=0.002$, $\lambda_{\max}(\mathcal{A})=200$, the right plot of Figure 3 reports the convergence history of the extreme Ritz values computed at the $k$-th iteration of both CG and the Galerkin method for $k=1,\ldots,82$ (the dashed lines indicate the target eigenvalues). The Ritz value tending to the largest eigenvalue converges in very few iterations. For each approach, the Ritz value approximating the smallest eigenvalue takes many more iterations to converge, and these iterations seem to match the stagnation phase observed in the left plot. It appears that the matrix Galerkin approximation space is able to implicitly capture the Kronecker structure of the eigenvector associated with $\lambda_{\min}({\cal A})$ much earlier than what CG can do by using the unstructured basis $Q_{k}$. Once again, this emphasizes the importance of the Kronecker basis determined by the matrix Galerkin method.

5 Petrov-Galerkin method and residual minimization

Whenever the linear operator $\mathcal{S}$ is not spd, the Galerkin method does not necessarily lead to a minimization of the error norm. As in the linear system setting, a numerical procedure fulfilling an optimality condition can be obtained by imposing a Petrov-Galerkin condition on the residual also when solving linear matrix equations. For the case of Lyapunov and Sylvester equations, this strategy has been already explored in, e.g., [30, 22], and in this section we are going to present some considerations about the application of Petrov-Galerkin methods to the solution of generic linear matrix equations of the form (2).

We first recall the Petrov-Galerkin framework applied to the solution of the linear system (1). If the columns of $\mathcal{V}_{m}\in\mathbb{R}^{N\times m}$ constitute an orthonormal basis for the selected trial space $\mathcal{K}_{m}$, we want to compute a solution $x_{m}=\mathcal{V}_{m}y_{m}$, where $y_{m}\in\mathbb{R}^{m}$ is calculated by imposing a Petrov-Galerkin condition on the residual vector $r_{m}=f-\mathcal{M}x_{m}$. In its full generality, such a condition reads

[TABLE]

where $\mathcal{L}_{m}$ is the chosen test space. See, e.g., [40, Chapter 5].

For the particular choice $\mathcal{L}_{m}=\mathcal{M}\mathcal{K}_{m}$, the condition in (15) is equivalent to computing $x_{m}$ as the minimizer of the residual norm over $\mathcal{K}_{m}$, namely

[TABLE]

See, e.g., [40, Proposition 5.3]. With the selection $\mathcal{K}_{m}=K_{m}(\mathcal{M},f)$, the minimization problem above can be significantly simplified by exploiting the Arnoldi relation; this is the foundation of some of the most popular minimal residual methods for linear systems such as, e.g., MINRES [34] and GMRES [41].

A similar approach can be pursued for the solution of linear matrix equations. Indeed, let $N=np$ and consider $V_{k}\in\mathbb{R}^{n\times k}$, $W_{k}\in\mathbb{R}^{p\times k}$ with full column rank 333Once again, the two matrices may have different column dimensions, that is $V_{k_{1}}\in\mathbb{R}^{n\times k_{1}}$, $W_{k_{2}}\in\mathbb{R}^{p\times k_{2}}$. For the sake of clarity in the exposition, we limit our presentation to the case $k_{1}=k=k_{2}$., and let $\text{range}(V_{k})$, $\text{range}(W_{k})$, be the corresponding left and right approximation spaces. With $\mathcal{S}_{\ell}$ as in Definition 2.1, we can formally set $\mathcal{K}_{m}=\text{range}(W_{k}\otimes V_{k})$ and $\mathcal{L}_{m}=\mathcal{S}_{\ell}\mathcal{K}_{m}$. An approximate solution in the form $X_{k}=V_{k}Y_{k}W_{k}^{T}$, with $Y_{k}\in\mathbb{R}^{k\times k}$, can be determined by imposing the condition (15) to the vector form of the residual matrix $R_{k}=\mathcal{S}(V_{k}Y_{k}W_{k}^{T})-F$.

Petrov-Galerkin methods for (2) thus seek a solution $X_{k}=V_{k}Y_{k}W_{k}^{T}$ by solving

[TABLE]

that is

[TABLE]

In spite of their appealing minimization property, minimal residual methods are not very popular in the matrix equation literature. This is mainly due to the difficulty in dealing with the numerical solution of the minimization problem (16). In general, one can apply an operator-oriented (preconditioned) CG method to the normal equations as

[TABLE]

where $\mathcal{S}^{*}$ is the adjoint of $\mathcal{S}$, namely

[TABLE]

If $\text{range}(V_{k})$ and $\text{range}(W_{k})$ are general spaces, the solution of (17) can be very expensive in terms of both computational time and memory requirements.

In [30], the authors consider (17) in the case of the Lyapunov equation (3) with $F$ low-rank and negative semidefinite. In particular, if $F=-bb^{T},$ $b\in\mathbb{R}^{n\times q}$, $q\ll n$, they employ the approximation spaces $\text{range}(V_{k})=\text{range}(W_{k})$ such that $b=V_{1}L_{b}$ for some $L_{b}\in\mathbb{R}^{q\times q}$, $q=\text{rank}(C)$, and satisfying an Arnoldi-like relation of the form

[TABLE]

for $[V_{k},\breve{V}_{k+1}]\in\mathbb{R}^{n\times(k+1)q}$ having orthonormal columns and $\underline{H}_{k}\in\mathbb{R}^{(k+1)q\times kq}$. In this case, the minimization problem (16) can be written as

[TABLE]

and three different methods for its solution are illustrated.

If the coefficient matrix $A$ in (3) is stable (antistable) and $F$ is symmetric negative semidefinite, the exact solution $X$ is symmetric positive (negative) semidefinite. See, e.g., [47]. However, as reported in [30], the numerical solution $X_{k}=V_{k}Y_{k}V_{k}^{T}$ is not guaranteed to be semidefinite if $Y_{k}$ is computed as in (18).

In [30, Section 3.4] it is shown that (18) is equivalent to computing $Y_{k}$ as the solution of the generalized Sylvester equation

[TABLE]

where

[TABLE]

so that $D$ is symmetric but indefinite. This is one of the main obstacles in proving the semidefiniteness of $Y_{k}$ through the matrix formulation (19). Without further hypotheses, the symmetric matrix $Y_{k}$ solving (19) is indefinite in general, thus preventing $Y_{k}$ from preserving the semidefiniteness property of the solution to be approximated.

From a computational viewpoint, if resorting to a Kronecker form is excluded, the generalized Sylvester equation (19) can be solved by means of the methods described in [30] and its references. In addition, setting $\mathfrak{L}(Z)=\underline{H}_{k}^{T}\underline{H}_{k}Z+Z\underline{H}_{k}^{T}\underline{H}_{k}$ and $\mathfrak{N}(Z)=H_{k}ZH_{k}+H_{k}^{T}ZH_{k}^{T}$, fixed point iterations can be used whenever the spectral radius of the operator $\mathfrak{L}^{-1}(\mathfrak{N}(\cdot))$ is less than one; see, e.g., [17, 42, 23] for various implementations.

5.1 A constrained residual minimization approach for Lyapunov equations

To cope with the lack of semidefiniteness in the least squares problem approach, we propose to explicitly impose the semidefiniteness as a constraint. For instance, if a negative semidefinite solution is sought, the problem becomes

[TABLE]

To numerically solve this inequality constrained least squares problem, we consider a linear matrix inequalities (LMI) approach, which suits very well the matrix equation framework [15, 46]; other general purpose methods could also be considered [1, 32].

In the LMI context, (20) can be stated as the following semidefiniteness matrix inequalities

[TABLE]

for the unknown matrix $Y$ and scalar $\gamma>0$; here $J=[I_{kq};0]$ and $M=[L_{b}L_{b}^{T},0;0,0]$.

Example 5.1**.**

We consider the Lyapunov equation (3) with $A=QDQ^{-1}\in\mathbb{R}^{n\times n}$, $D$ as in Example 4.1, $Q$ a random matrix, and $F=-bb^{T}$, where $b\in\mathbb{R}^{n}$ is a random vector with unit norm.

Since $A$ is antistable and the right-hand side is symmetric negative semidefinite, the solution $X$ is symmetric negative semidefinite and we thus expect the approximate solution $X_{k}=V_{k}Y_{k}V_{k}^{T}$ to be so as well.

We apply the Petrov-Galerkin method discussed in this section in the solution process and we adopt the Krylov subspace as approximation space, i.e., $\text{range}(V_{k})=K_{k}(A,b)$. The matrix $Y_{k}$ is computed in two different ways. We first solve the unconstrained minimization problem (18) getting the matrix $Y_{k}^{\text{uncon}}$. In particular, $Y_{k}^{\text{uncon}}$ is computed by applying a (preconditioned) CG method to the matrix equation (19). See, e.g., [30]. Then, we compute $Y_{k}^{\text{const}}$ by solving the constrained minimization problem (20). The Petrov-Galerkin method is stopped as soon as the relative residual norm becomes smaller than $10^{-6}$.

In Figure 4 we plot the intervals $[\min_{j}\{\lambda_{j}(Y_{k}^{\text{uncon}})\geqslant 0\},\max_{j}\{\lambda_{j}(Y_{k}^{\text{uncon}})\geqslant 0\}]$ of the undesired positive eigenvalues of $Y_{k}^{\text{uncon}}$ for all $k$ for the case $n=1000$. For $k=1,2$, $Y_{k}^{\text{uncon}}$ has all negative eigenvalues, while it starts being indefinite for $k\geqslant 3$ so that $X_{k}^{\text{uncon}}=V_{k}Y_{k}^{\text{uncon}}V_{k}^{T}$ is an indefinite approximation to the negative semidefinite $X$. Nonetheless, for $k=68$, the (undesired) positive eigenvalues of $X_{k}^{\text{uncon}}$ are small enough so as to still allow a sufficiently accurate approximation, in terms of relative residual norm. On the other hand, this problem is not encountered with $Y_{k}^{\text{constr}}$, thanks to the explicit negative semidefiniteness constraint in the formulation (20).

From the legend of Figure 4, we can see that the number of positive eigenvalues of $Y_{k}^{\text{uncon}}$ increases as the iteration proceed, even though they diminish in magnitude. The latter trend is not surprising. Indeed, even if $Y_{k}^{\text{uncon}}$ is computed by (18), the Petrov-Galerkin method is converging towards the negative semidefinite solution $X$ and, for an approximation space spanning the whole ${\mathbb{R}}^{n}$, the method would retrieve the exact solution, regardless of the minimization problem (18).

We would like to point out that both tested variants of the Petrov-Galerkin method needed 68 iterations to converge and the actual values of the residual norm provided by (18) and (20) were always very similar to each other, during the whole convergence history. This phenomenon surely deserves further studies as, in principle, (20) leads to a residual norm that is greater or equal than the one provided by (18), while the two solutions (constrained and unconstrained) do not necessarily have to be close to each other.

In our computational experiments, we have used the Yalmip software [31] running the algorithm Sedumi in Matlab [50]. This algorithm is rather expensive and computing the solution $Y_{k}$ to (20) at each Krylov iteration $k$ often leads to a very time consuming solution procedure. We think this issue can be fixed in different ways. For instance, one may compute $Y_{k}$, and thus check the residual norm, only periodically, say every $d\geqslant 1$ iterations. Moreover, the explicit solution $Y_{k}$ is required only at convergence while we just need the value of the residual norm during the Krylov routine. It may be possible to compute such a residual norm without calculating the whole $Y_{k}$ as it is done in [36] for the Galerkin method and in [30] for the Petrov-Galerkin technique equipped with the unconstrained minimization problem (18).

The study of the aforementioned enhancements and, more in general, the employment of constrained minimization procedures in the solution of linear matrix equations will be the topic of future research.

6 Conclusions

We have shown that the optimality properties of Galerkin and Petrov-Galerkin methods naturally extend to the general linear matrix equation setting. Such features do not depend on the adopted approximation spaces even though, in actual computations, fast convergence depends on the suitable subspace selection. Identifying effective subspaces for general (multiterm) linear matrix equations depends on the problem at hand, and it may seem easier to recast the solution in terms of a large vector linear system. On the other hand, the vector form can be extremely memory consuming, while the vector linear system encodes some spectral redundancy which may cause a delay in the converge of the adopted iterative solution scheme.

Petrov-Galerkin schemes require to solve a matrix minimization problem at each iteration and we have suggested to explicitly incorporate a semidefiniteness constraint in its formulation. To the best of our knowledge, such approach has never been proposed in the literature and the employment of constrained optimization techniques in the context of Petrov-Galerkin methods for linear matrix equations opens many new research directions.

Acknowledgements

Both authors are members of the Italian INdAM Research group GNCS.

We thank the two anonymous reviewers for their insightful remarks.