LME: Linear Matrix Equation#
The Linear Matrix Equation (LME) solver object encapsulates functionality intended for the solution of linear matrix equations such as Lyapunov, Sylvester, and their various generalizations. A matrix equation is an equation where the unknown and the coefficients are all matrices, and it is linear if there are no nonlinear terms involving the unknown. This excludes some types of equations such as the Algebraic Riccati Equation or other quadratic matrix equations, which are not considered here.
The general form of a linear matrix equation is
where \(X\) is the solution.
Since SLEPc is mainly concerned with iterative methods, so is the LME module. This implies that LME can address only the case where the solution \(X\) has low rank. In many situations, \(X\) does not have low rank, which means that the LME solver, if it converges, will truncate the solution to an artificially small rank, with a large approximation error.
Warning
The LME module should be considered experimental. As explained below, the currently implemented functionality is quite limited, and may be extended in the future.
Current Functionality#
The user interface of the LME module is prepared for different types of equations, see a list in table Problem types considered in LME. However, currently there is only basic support for continuous-time Lyapunov equations, and the rest are just a wish list.
Problem Type |
Equation |
\(A\) |
\(B\) |
\(D\) |
\(E\) |
|
|---|---|---|---|---|---|---|
Continuous-Time Lyapunov |
\(AX+XA^*=-C\) |
yes |
\(A^*\) |
- |
- |
|
Sylvester |
\(AX+XB=C\) |
yes |
yes |
- |
- |
|
Generalized Lyapunov |
\(AXD^*+DXA^*=-C\) |
yes |
\(A^*\) |
yes |
\(D^*\) |
|
Generalized Sylvester |
\(AXE+DXB=C\) |
yes |
yes |
yes |
yes |
|
Discrete-Time Lyapunov |
\(AXA^*-X=-C\) |
yes |
- |
- |
\(A^*\) |
|
Stein |
\(AXE-X=-C\) |
yes |
- |
- |
yes |
Continuous-Time Lyapunov Equation#
Given two matrices \(A,C\in\mathbb{C}^{n\times n}\), with \(C\) Hermitian positive-definite, the continuous-time Lyapunov equation, assuming that the right-hand side matrix \(C\) has low rank, can be expressed as
where we have written the right-hand side in factored form, \(C=C_1C_1^*\), \(C_1\in\mathbb{C}^{n\times r}\) with \(r = \operatorname{rank}(C)\). In general, \(X\) will not have low rank, but our intention is to compute a low-rank approximation \(\tilde X=X_1X_1^*\), \(X_1\in\mathbb{C}^{n\times p}\) for a given value \(p\). Note that the equation has a unique solution if and only if \(A\) is stable, i.e., all of its eigenvalues have strictly negative real part.
The reason why SLEPc provides a solver for equation (2) is that it is required in the eigensolver EPSLYAPII, which implements the Lyapunov inverse iteration [Elman and Wu, 2013, Meerbergen and Spence, 2010].
Krylov Solver#
Currently, the solver for the continuous-time Lyapunov equation available in SLEPc is based on the following procedure for each column \(c\) of \(C_1\):
Build an Arnoldi factorization for the Krylov subspace \(\mathcal{K}_m(A,c)\), \(AV_m=V_mH_m+h_{m+1,m}v_{m+1}e_m^*\).
Solve the compressed Lyapunov equation, \(H_mY+YH_m^*=-\tilde c\,\tilde c^* \quad\mathrm{with}\quad\tilde c=V_m^*c\).
Set the approximate solution to \(X_m=V_mYV_m^*\).
Furthermore, our Krylov solver incorporates an Eiermann-Ernst-type restart as proposed by Kressner [2008]. The restart will discard the Arnoldi basis \(V_m\) but keep \(H_m\), then continue the Arnoldi recurrence from \(v_{m+1}\). The upper Hessenberg matrices generated at each restart are glued together. This implies that at each restart the size of the compressed Lyapunov equation grows. After convergence, once the full \(Y\) is available, we perform a second pass to reconstruct the successive \(V_m\) bases.
The restarted algorithm is as follows:
for \(k=1,2,\ldots\)
Run Arnoldi \(AV_m^{(k)}=V_m^{(k)}H_m^{(k)}+h_{m+1,m}^{(k)}v_{m+1}^{(k)}e_m^*\).
Set \(H_{km}=\left[\begin{array}{cc}H_{(k-1)m}&0\\h_{m+1,m}^{(k-1)}e_1e_{(k-1)m}^*&H_m^{(k)}\end{array}\right]\).
Solve compressed equation \(H_{km}Y+YH_{km}^*=-\|c\|_2^2\,e_1e_1^*\).
Compute residual norm \(\|R_k\|_2=h_{m+1,m}^{(k)}\|e_{km}^*Y\|_2\).
end
The residual norm is used for the stopping criterion.
Projected Problem#
As mentioned in the previous section, the projected problem is a smaller Lyapunov equation, which must be solved in each outer iteration of the method (restart). The solution of this small dense matrix equation is also implemented within the LME module (there is no auxiliary class for this).
If \(A\) is stable, then \(X\) is symmetric positive semidefinite. Hence, looking at step 3 of the algorithm in the previous section, it is better to compute the Cholesky factor \(Y=LL^*\) to obtain the low rank factor \(X_1=V_mL\). The method of Hammarling will compute \(L\) directly, without computing \(Y\) first.
If SLEPc has been configured with the option --download-slicot (see Wrappers to External Libraries), it will be possible to use SLICOT subroutines to apply Hammarling’s method, otherwise a more basic algorithm will be used.
Note
Most of the SLICOT subroutines are implemented for real arithmetic only, so it is not possible to use them in a complex build of PETSc/SLEPc.
In the restarted algorithm described above, in the second pass to recalculate the successive \(V_m\)’s, the update is based on the truncated SVD of the Cholesky factor
Basic Usage#
Once the solver object has been created with LMECreate(), the user has to define the matrix equation to be solved. First, the type of equation must be set with LMESetProblemType(), choosing one from table Problem types considered in LME.
The coefficient matrices \(A\), \(B\), \(D\), \(E\) must be provided via LMESetCoefficients(), but some of them are optional depending on the matrix equation. For Lyapunov equations, only \(A\) must be set, which is normally a large, sparse matrix stored in MATAIJ format. Note that in table Problem types considered in LME, the notation \(A^*\) means that this matrix need not be passed, but the user may choose to pass an explicit transpose of matrix \(A\) (for improved efficiency). Also note that some of the equation types impose restrictions on the properties of the coefficient matrices (e.g., \(A\) stable in Lyapunov equations) and possibly on the right-hand side \(C\).
To conclude the definition of the equation, we must pass the right-hand side \(C\) with LMESetRHS(). Note that \(C\) is expressed as the outer product of a tall-skinny matrix, \(C=C_1C_1^*\). This can be represented in PETSc using a special type of matrix, MATLRC, that can be created with MatCreateLRC passing a dense matrix \(C_1\).
The call to LMESolve() will run the solver to compute the solution \(X\). The rank of \(X\) may be prescribed by the user or selected dynamically by the solver. Next, we discuss these two options:
To prescribe a fixed rank for \(X\), we must preallocate it, creating a MATLRC matrix, and then pass it via
LMESetSolution()before callingLMESolve(). In this way, the solver will be restricted to a rank equal to the number of columns provided in \(X_1\).If we do not call
LMESetSolution()beforehand, the solver will select the rank dynamically. Then afterLMESolve()we must callLMEGetSolution()to retrieve a MATLRC matrix representing \(X\).
References
H. Elman and M. Wu. Lyapunov inverse iteration for computing a few rightmost eigenvalues of large generalized eigenvalue problems. SIAM J. Matrix Anal. Appl., 34(4):1685–1707, 2013. doi:10.1137/120897468.
D. Kressner. Memory-efficient Krylov subspace techniques for solving large-scale Lyapunov equations. In 2008 IEEE International Conference on Computer-Aided Control Systems, 613–618. IEEE, 2008. doi:10.1109/cacsd.2008.4627370.
K. Meerbergen and A. Spence. Inverse iteration for purely imaginary eigenvalues with application to the detection of Hopf bifurcations in large-scale problems. SIAM J. Matrix Anal. Appl., 31(4):1982–1999, 2010. doi:10.1137/080742890.