Residual-based iterations for the generalized Lyapunov equation
Tobias Breiten Emil Ringh







Tobias Breiten Institute for Mathematics and Scientific Computing, Karl-Franzens-Universität
 
              Heinrichstrasse 36, 8010 Graz, Austria

              tobias.breiten@uni-graz.at
           Emil Ringh Department of Mathematics, KTH Royal Institute of Technology

              Lindstedtsvägen 25, 100 44 Stockholm, Sweden

              eringh@kth.se




Residual-based iterations for the generalized Lyapunov equation
    Tobias Breiten Emil Ringh

    December 30, 2023
===============================================================




This paper treats iterative solution methods to the generalized Lyapunov equation. 
Specifically it expands the existing theoretical justification for the alternating linear scheme (ALS) from the stable Lyapunov equation to the stable generalized Lyapunov equation. 
Moreover, connections between the energy-norm minimization in ALS and the theory to H2-optimality of an associated bilinear control system are established.
It is also shown that a certain ALS-based iteration can be seen as iteratively constructing rank-1 model reduction subspaces for bilinear control systems associated with the residual. 
Similar to the ALS-based iteration, the fixed-point iteration can also be seen as a residual-based method minimizing an upper bound of the associated energy norm. 
Lastly a residual-based generalized rational-Krylov-type subspace is proposed for the generalized Lyapunov equation.


65F10 58E25 65F30 65F35







§ INTRODUCTION










This paper concerns iterative ways to compute approximate solutions to what has become known as  the generalized Lyapunov equation

    (X) + Π(X) + BB^T = 0,

where X∈^n× n is unknown, and the operators , Π: ^n× n→^n× n are defined as

    (X) := AX + XA^T 
    Π(X) := ∑_i=1^m N_iXN_i^T,

with A, N_i∈^n× n for i=1,…,m given.
The operator  is commonly known as the Lyapunov operator, and Π is sometimes called a correction.
We further assume that A is stable, i.e., A has all its eigenvalues in the left-half plane, which implies that  is invertible <cit.>. Moreover, we assumed that ρ(^-1Π)<1, where ρ denotes the (operator) spectral radius. The assumption on the spectral radius implies that (<ref>) has a unique solution, see, e.g., <cit.>. Furthermore, the definition of Π in (<ref>) implies that it is non-negative, in the sense that Π(X) is positive semidefinite when X is positive semidefinite. Thus one can assert that the unique solution X is indeed positive definite, see, e.g., <cit.> or <cit.>.



 §.§ Related work

The standard Lyapunov equation, AX+XA^T+BB^T = 0, has been well studied for a long time and considerable research effort has been, and is still, put into finding efficient algorithms for computing the solution.
Some examples of methods are the classical  Bartels-Stewart algorithm <cit.> for small and dense problems based on factorizing the matrix A and backward substitution, the Smith method presented in <cit.>, and the Riemannian optimization method <cit.> which computes a low-rank approximation by minimizing an associated cost function over the manifold of rank-k matrices, where k≪ n.
The Lyapunov equation has a lose connection to control theoretic methods,

such as the iterative rational Krylov algorithm (IRKA) <cit.> which computes locally 2̋-optimal reduced order systems.
Related research is presented in a series of paper <cit.>, where Druskin and co-authors develop a strategy to choose shifts for the rational Krylov subspace for efficient subspace reduction when solving PDEs <cit.>, as well as for model reduction of linear single-input-single-output (SISO) systems and solutions to Lyapunov equations <cit.>. Instead of computing full spaces iteratively with a method such as IRKA, the idea is to construct an infinite sequence with asymptotically optimal convergence speed <cit.>. Then the subspace can be dynamically extended as needed, until required precision is achieved, rather than starting the process by deciding the dimension of the space, an a priori parameter choice with (complex) implications of the final approximation accuracy. The idea is also further developed by using tangential directions, proving especially useful for situations where the right-hand side is not of particularly low rank <cit.>, e.g., multiple-input-multiple-output (MIMO) systems.
For a more complete overview of results and techniques for Lyapunov equations, see, e.g., the review article <cit.>.

The generalized Lyapunov equation (<ref>) has received increased attention over the last decade. Results on low-rank approximability has emerged <cit.><cit.>, i.e., similarly to the standard Lyapunov equation one can in certain cases when the right-hand-side B is of low rank, r≪ n, expect the singular values of the solution to decay rapidly even for the generalized Lyapunov equation. Such low-rank approximability is useful for large-scale problems since algorithms can be adapted to exploit the low-rank format, reducing computational effort and storage requirement. Example of algorithms exploiting low-rank structures are a Bilinear ADI method <cit.>, specializations of Krylov methods for matrix equations <cit.>, as well as greedy low-rank methods <cit.>, and exploitations of the fixed-point iteration <cit.>.
Through the connection with bilinear control systems there is an extension of IRKA, known as bilinear iterative rational Krylov (BIRKA) <cit.>.
There are also methods based on Lyapunov and ADI-preconditioned GMRES and BICGSTAB <cit.>, and in general for problems with tensor product structure <cit.>.
We also mention that in the case when the correction Π has low operator-rank, typically if all the matrices N_i are of low rank, there is a specialization of the Sherman-Morrison-Woodbury formula to the linear matrix equation, see, e.g. <cit.>.




 §.§ Outline and summary of contributions

The paper is outlined as follows: In Section <ref> we present some existing theory and preliminary results which sets the context of the paper. This is followed up in Section <ref> where we prove that the alternating linear scheme (ALS) presented by Kressner and Sirković in <cit.> computes search directions which at each step fulfill a first order necessary condition for being 2̋-optimal, in a certain sense. Moreover, we show equivalence between ALS and BIRKA. In Section <ref> we make an analogue to the fixed-point iteration, showing that it minimizes an upper bound of the energy-norm,  before we in Section <ref> present a residual-based generalized rational-Krylov-type subspace adapted for solving the generalized Lyapunov equation. While the main focus of this paper is  on solution methods to the generalized Lyapunov equation, in Section <ref> we draw parallels with rational Krylov subspaces for the standard Lyapunov equation. We end the paper with some numerical experiments presented in Section <ref>, and conclusions and outlooks in Section <ref>.






§ PRELIMINARIES



 §.§ Generalized matrix equations and approximations

In this paper we are concerned with iterative methods for computing approximative solutions to the generalized Lyapunov equation (<ref>). In this section we recall some basic definitions and results that will be used in later in the paper.
Regarding the generalized Lyapunov equation there is one special class of problems where more can be said, and that is the symmetric case.

We call equation (<ref>) a symmetric generalized Lyapunov equation if A = A^T and N_i=N_i^T for i=1,…,m.


In general we will think of X̂_k∈^n× n as an approximation of the solution to (<ref>), where k is an iteration count. The goal is to find an X̂_k such that X-X̂_k is small for some norm, where X is the exact solution to (<ref>). However, in many situations there is no precise definition, and no need for one, of what is meant by approximation. Nevertheless, in order to discuss projection methods and make the results precise, 
we make the following (standard) definition of a Galerkin approximation.


Let _k⊆^n be an n_k≤ n dimensional subspace for k=0,1,…, and let _k∈^n× n_k be a matrix containing an orthogonal basis of _k.
We call X̂_k the Galerkin approximation to (<ref>), in _k, if X̂_k = _k Y_k _k^T and Y_k is determined by the condition

    _k^T((X̂_k) + Π(X̂_k)  + BB^T)_k = 0,

known as the projected problem.

It is possible that the projected problem does not have a solution, or that the solution is not unique, and in such case the definition of (the) Galerkin approximation in Definition <ref> is nonsensical. For the generalized Lyapunov equation there are certain sufficient conditions for the Galerkin approximation to exist and be unique, e.g., the criteria in <cit.><cit.> or <cit.>.
In fact, for the symmetric generalized Lyapunov equation with our assumptions on the spectral radius, it is possible to assure the existence of the Galerkin approximation. Since the property <cit.> is preserved by orthogonal projections in the symmetric case, cf. Proposition <ref>. 
However, we will not delve deeper into this topic, and further on simply assume that the projected problem has a unique solution.

Related to the (Galerkin) approximation we also have the (Galerkin) residual.

Let _k and _k be as in Definition <ref> and, assuming that it exists, let X̂_k be the Galerkin approximation in _k to (<ref>). We define the Galerkin residual as 

    _k := (X̂_k) + Π(X̂_k) + BB^T.


By using this definition the projected problem is commonly expressed as _k^T_k_k=0, which is also known as the Galerkin orthogonality condition and is a characteristic feature of the Galerkin residual.





It is possible to define a generic residual, which we with a slight abuse of notation also call _k. We define it as: for any matrix X̂_k ∈^n× n, _k := (X̂_k) + Π(X̂_k) + BB^T. Some of the results and arguments presented are valid for a (generic) residual and others, more specialized, only for the Galerkin residual, or another residual specific to some approximation scheme. However, we will always connect the residual to a specific approximation and the Galerkin residual will always be referenced as such.


We end this section with a specialization of a classical result from linear algebra, sometimes called the residual equation.


Let  and Π be as in (<ref>)-(<ref>), and let X∈^n× n be the solution to (<ref>). Moreover, let X̂_k ∈^n× n be a given matrix, let _k be the (generic) residual, and define the error X^e_k := X - X̂_k. Then

    (X^e_k) + Π(X^e_k) + _k = 0.












The result in Proposition <ref> is useful since it may serve as a starting point to derive iterative algorithms for computing approximations to (<ref>). 

Hence, one strategy for computing updates to the current iterate is to compute, or approximate, X^e_k, which Proposition <ref> allows us to do based on the known, or computable, quantities , Π and _k.

The result for the Lyapunov equation was presented already by Smith in <cit.>, and for generalized matrix equations cf. <cit.>, and <cit.>.





 §.§ Bilinear systems

We recall some control theoretic concepts for bilinear control systems of the form 

    Σ{ẋ(t)    = Ax(t) + ∑_i=1^m N_i x(t) w_i(t) + Bu(t)
    
    y(t)    = C x(t),
    .

with A,N_i ∈^n× n, B∈^n× r and C∈^r× n and control inputs u(t)∈^r and w(t)∈^m.


We highlight that the bilinear system Σ in (<ref>) differs from the usual notation used in the  literature, see e.g.,  <cit.>. The formulation (<ref>) is convenient since it allows for m≠ r. However, the system Σ can be put on the usual form by considering the input vector [w(t)^T u(t)^T]^T, adding m zero-columns to the beginning of B, i.e., considering [0, B], and considering the matrices N_i = 0 for i=m+1,m+2,…,m+r. Thus it is reasonable to consider the same Gramians. The system Σ can also be compared to systems from applications, e.g., <cit.>.

As in  <cit.>, for a MIMO bilinear system (<ref>), we define the controllability and observability Gramian as follows.

Let A∈^n× n be a stable matrix, N_1,…,N_m∈^n× n, B∈^n× r, and C∈^r× n, and consider the bilinear system (<ref>).
Moreover, let P_1(t_1) := e^At_1B, P_j(t_1,…,t_j) := e^At_j[N_1P_j-1, …, N_mP_j-1] for j=2,3,…, Q_1(t_1) := Ce^At_1, and Q_j(t_1,…,t_j) := [N_1^TQ_j-1^T,…,N_m^TQ_j-1^T]^T e^At_j for j=2,3,…. We define the controllability and observability Gramian respectively as

    P    := ∑_j=1^∞∫_0^∞…∫_0^∞ P_j P_j^T dt_1… dt_j
    
    Q    := ∑_j=1^∞∫_0^∞…∫_0^∞ Q_j^T Q_j dt_1… dt_j.





It is possible that the generalized Gramians from Definition <ref> do not exist; sufficient conditions are given in, e.g., <cit.>. However, if the Gramians exist we also know that they satisfy the following matrix equations 

    AP + PA^T + ∑_i=1^m N_i P N_i^T +BB^T    = 0
    
    A^TQ + QA + ∑_i=1^m N_i^T Q N_i +C^T C    = 0.











In relation to the generalized controllability and observability Gramians, one might also define a generalized cross Gramian similar to the SISO case discussed in <cit.>. 
Consider the generalized Lyapunov equation (<ref>), and an approximation X̂_k with related error X_k^e and residual _k.
One can easily verify that for the auxiliary system 

    Σ^e = {ẋ(t)    = Ax(t) + ∑_i=1^m N_i x(t) w_i(t) + B__ku(t)
    
    y(t)    = C__k x(t),
    .

where B__k = US^1/2, and C__k = S^1/2V^T, with _k = USV^T being a singular value decomposition of _k, the associated cross Gramian coincides with the error X_k^e. In the special case of a symmetric system, we additionally have the following result.






Consider a symmetric generalized Lyapunov equation. Let X̂_k be an approximation such that the residual _k = _k^T ≽ 0. Then one can choose B__k=C__k^T and the error X^e_k is the controllability and observability Gramian of the auxiliary bilinear system.


For what follows, let us recall the definition of the ℋ_2-norm for bilinear systems that was introduced by Zhang and Lam in <cit.>.


We define the 2̋-norm for the bilinear system (<ref>) as
  
    Σ_2̋^2:=(∑_j=1^∞∫_0^∞…∫_0^∞∑_ℓ_1,⋯,ℓ_j=1^mg_j^(ℓ_1,…,
      ℓ_j)(g_j^(ℓ_1,…,ℓ_j))^T ds_1 ⋯ds_j ) ,

  with g_j^(ℓ_1,…,ℓ_j)(s_1,…,s_k)
  :=Ce^As_jN_ℓ_1e^As_j-1 N_ℓ_2⋯ e^As_1b_ℓ_j.
 
It has been shown in <cit.> that if the Gramians P and Q from Definition <ref> exist, then it holds that

    Σ_2̋^2 = (CPC^T) = (B^T Q B).









§ ALS AND 2̋-OPTIMAL MODEL REDUCTION FOR BILINEAR SYSTEMS


In this section, we discuss a low rank approximation method proposed by Kressner and Sirković in <cit.>. We show that several results can be generalized from the case of standard Lyapunov equations to the more general form (<ref>). Moreover, we show that in the symmetric case the method allows for an interpretation in terms of 2̋-optimal model reduction for bilinear control systems. With this in mind, we assume that A=A^T ≺ 0 and N_i=N_i^T for i=1,…,m. If additionally we have that ρ(^-1Π)<1, the operator ℳ(X):=-(X)-Π(X) allows us to define a weighted inner product via

    ⟨·,·⟩_ℳ   ^n× n×^n × n→
    ⟨ X,Y ⟩_ℳ   = ⟨ X,ℳ(Y) ⟩ = (X^T ℳ(Y) ),

with a corresponding induced ℳ-norm, also known as energy norm,

    X^2 = ⟨ X,X ⟩_ℳ.




 §.§ ALS for the generalized Lyapunov equation


In <cit.>, it is suggested to construct iterative approximations X̂_k by rank-1 updates that are locally optimal with respect to the ℳ-norm. To be more precise, assume that X is a solution to the symmetric Lyapunov equation (<ref>), i.e.,

    AX + XA + ∑_i=1^m N_i X N_i + BB^T=0.

Given an approximation X̂_k, we consider the minimization problem 

    min_v,w∈ℝ^n X-X̂_k - v w^T_ℳ^2    = ⟨ X-X̂_k - v w^T, X-X̂_k - v w^T ⟩_ℳ.

Since the minimization involves the constant term X-X̂_k^2_ℳ, it suffices to focus on

    J(v,w):= ⟨ vw^T , vw^T ⟩_ℳ - 2 ( wv^T _k),

where _k is the current residual, i.e.,

    _k:= A X̂_k + X̂_k A + ∑_i=1^m N_i X̂_k N_i + BB^T.

Locally optimal vectors v_k and w_k are then (approximately) determined via an alternating linear scheme (ALS). The main step is to fix one of the two vectors, e.g., v and then minimize the strictly convex objective function in order to obtain an update for w. A pseudocode of the algorithm presented in <cit.> is given in Algorithm <ref>.

The ALS-based approach for computing new subspace extensions can also, in terms of Proposition <ref>, be seen as searching for an approximation to X^e_k of the form v_kw_k^T by iterating ((v_kw_k^T) + Π(v_kw_k^T) + _k)w_k = 0 when determining v_k and v_k^T((v_kw_k^T) + Π(v_kw_k^T) + _k) = 0 when determining w_k. This is to say that the error is approximated by a rank-1 matrix, and at convergence this would result in the new residual, _k+1, being left-orthogonal to v_k and right-orthogonal to w_k.
In the symmetric case, local minimizers of (<ref>) are necessarily symmetric positive semidefinite. This yields the following extension of <cit.>.


  Consider the symmetric generalized Lyapunov equation (<ref>) and assume that A≺ 0, ρ(^-1Π)<1, and _k=_k^T≽0. Let J be as in (<ref>).
 Then every local minimum (v_*,w_*) to J is such that v_*w_*^T is symmetric positive semidefinite.


The proof, naturally, follows along the lines of <cit.>, and hence without loss of generality we assume that v_*≠ 0 , w_*≠ 0, and v_*=w_*. Thus v_*w_*^T is positive semidefinite if and only if v_*=w_*. The proof is by contradiction and we assume that v_*≠ w_*. Then from the strict convexity of J we have that 
    2 J(v_*,w_*)    < J(v_*,v_*) + J(w_*,w_*).

  Simplifying the left-hand-side we get
  
    2 J(v_*,w_*) = 
         
     -2((v_*w_*^T)^T(v_*w_*^T)) -2((v_*w_*^T)^TΠ(v_*w_*^T)) - 4((v_*w_*^T)^T_k)
     
       =
     -2 v_*^T(v_*w_*^T)v_* -2 v_*^TΠ(v_*w_*^T)w_* - 4 v_*^T_k w_*,

 and similarly the right-hand-side gives
 
    J(v_*,v_*) + J(w_*,w_*) =
        - v_*^T(v_*v_*^T)v_* - v_*^TΠ(v_*v_*^T)v_* - 2 v_*^T_k v_*
     
       - w_*^T(w_*w_*^T)w_* - w_*^TΠ(w_*w_*^T)w_* - 2 w_*^T_k w_*.

 Collecting the terms involving the  -operator we observe that 
 
    -2 v_*^T(v_*w_*^T)w_* + v_*^T(v_*v_*^T)v_* + w_*^T(w_*w_*^T)w_*    =
    
     2(v_*^Tv_*)(w_*^TAw_* - v_*^TAv_*) + 2w_*^Tw_*(v_*^TAv_* - w_*^TAw_*)    = 0,

 and thus the inequality reduces to
 
    -2 v_*^TΠ(v_*w_*^T)w_* - 4 v_*^T_k w_* <
          - v_*^TΠ(v_*v_*^T)v_* - 2 v_*^T_k v_* - w_*^TΠ(w_*w_*^T)w_* - 2 w_*^T_k w_*.

 We rewrite the inequality  as
  
    -2 v_*^TΠ(v_*w_*^T)w_* + v_*^TΠ(v_*v_*^T)v_* + w_*^TΠ(w_*w_*^T)w_* 
     < -2 (v_*-w_*)^T _k (v_*-w_*).

 Now if we can show that 
 
    -2 v_*^TΠ(v_*w_*^T)w_* + v_*^TΠ(v_*v_*^T)v_* + w_*^TΠ(w_*w_*^T)w_* ≥ 0,

 this would imply that -2 (v_*-w_*)^T _k (v_*-w_*)>0 in contradiction to the positive semidefiniteness of _k.
The argument is hence concluded by considering the expression -2 v_*^TΠ(v_*w_*^T)w_* + v_*^TΠ(v_*v_*^T)v_* + w_*^TΠ(w_*w_*^T)w_*. We can without loss of generality consider m=1, hence only one N-matrix, since the following argument can be applied to all terms in the sum independently. Observe that
 
    -2 v_*^T N v_*w_*^T N w_* + v_*^TNv_*v_*^TNv_* + w_*^TNw_*w_*^TNw_* = 
     (v_*^TNv_* - w_*^TNw_*)^2 ≥ 0,

 which shows the desired inequality, and thus concludes the proof.
 
Algorithm <ref> and the argument in Lemma <ref> is based on a residual. However, if one considers X̂_k = 0, then _k = BB^T, and hence the result is applicable directly to any symmetric generalized Lyapunov equation. The focus on the residual in the previous results is natural since it leads to the following extension of <cit.> to the case of symmetric generalized Lyapunov equations.

  Consider the symmetric generalized Lyapunov equation (<ref>) with the additional assumptions that A≺ 0 and ρ(^-1Π)<1. Moreover, consider the sequence of approximations constructed as
  
    X̂_0    = 0
    X̂_k+1   = X̂_k + v_k+1v_k+1^T,      k = 0,1,…,
 
where v_k+1 is a locally optimal vector computed with ALS (Algorithm <ref>).
 Then _k=_k^T ≽ 0 for all k≥ 0.


We show the assertion by induction. It trivially holds that _0=_0^T≽ 0. Now assume that this is the case for some _k. 
From Lemma <ref> the local minimizers of (<ref>) are symmetric and hence X̂_k+1 as in (<ref>) is reasonable as defined. Moreover, since X̂_k+1 and the operators in (<ref>) are symmetric we know that _k+1 is symmetric. Thus what is left to show is that _k+1≽0.

We have that _k+1≽0 if and only if y^T_k+1y≥0 for all y∈^n. Hence take an arbitrary y∈^n and consider y^T_k+1y. We derive properties similar to <cit.>: 

Since (v_k+1,v_k+1) is a (local) minimizer of J(v,w), it also follows that v_k+1 is a (global) minimizer of the (convex) cost function

    J_w(v):=J(v,w)=⟨ vw^T,vw^T⟩ _ℳ - 2(wv^T_k),

where w=v_k+1. Note that the gradient ∇_vJ_w of J_w with respect to v is given by

    (∇_vJ_w)_i =2⟨ e_iw^T,vw^T ⟩_ℳ- 2e_i^T _k w .

Due to the optimality of v_k+1 with respect to J_v_k+1, first order optimality conditions then imply that

    -Av_k+1v_k+1^Tv_k+1-v_k+1v_k+1^TAv_k+1-∑_i=1^mN_iv_k+1v_k+1^TN_iv_k+1=_k v_k+1 .

Striking this equality with v_k+1^T  from the left implies that

    2  v_k+1^TA v_k+1v_k+1^2 = - v_k+1^T _k v_k+1 - ∑_i=1^m ( v_k+1^T N_i   v_k+1)^2 .

Based on this we can write the residual as


    y^T_k+1y = y^T_k y + y^T((v_k+1v_k+1^T) + Π(v_k+1v_k+1^T))y
    
       =
     y^T_k y + y^T(Av_k+1v_k+1^T + v_k+1v_k+1^TA + ∑_i=1^m N_iv_k+1v_k+1^T N_i)y
    
       =
     y^T_k y + ∑_i=1^m y^TN_iv_k+1v_k+1^T N_i y + 1/v_k+1^2(
    
          +  y^T(-_k v_k+1v_k+1^T - ( v_k+1^TA v_k+1)  v_k+1v_k+1^T - ∑_i=1^m ( v_k+1^T N_i v_k+1)N_iv_k+1v_k+1^T)y
    
          +  y^T(-v_k+1  v_k+1^T_k - v_k+1v_k+1^T ( v_k+1^TAv_k+1)  - ∑_i=1^m ( v_k+1^T N_iv_k+1)v_k+1v_k+1^T N_i)y).

Here, the third equality follows from (<ref>) and its transpose by exploiting the symmetry of the involved matrices. We rearrange, identify the term -2( v_k+1^T A  v_k+1)  v_k+1 v_k+1^T and insert (<ref>) to get


    y^T_k+1y
    =
     y^T_k y
    
          + 1/v_k+1^2y^T(-_k v_k+1v_k+1^T -v_k+1v_k+1^T_k + 1/v_k+1^2v_k+1^T_kv_k+1 v_k+1 v_k+1^T )y
    
          +1/v_k+1^2 y^T (∑_i=1^m N_i v_k+1v_k+1^TN_i v_k+1^2 + 1/v_k+1^2∑_i=1^m ( v_k+1^TN_i v_k+1)^2 v_k+1v_k+1^T) y
    
          +1/v_k+1^2 y^T( - ∑_i=1^m ( v_k+1^T N_i v_k+1)N_iv_k+1v_k+1^T - ∑_i=1^m ( v_k+1^T N_iv_k+1)v_k+1v_k+1^T N_i)y [0]
    
       =
    
     y^T_k y 
    
          + 1/v_k+1^2y^T(-_k v_k+1v_k+1^T -v_k+1v_k+1^T_k + 1/v_k+1^2v_k+1^T_kv_k+1 v_k+1 v_k+1^T )y
    
          +  1/v_k+1^2(∑_i=1^m (y^TN_iv_k+1)^2v_k+1^2 +  1/v_k+1^2 ( v_k+1^TN_i v_k+1)^2 ( v_k+1^Ty)^2 
        - 2 ( v_k+1^T N_i v_k+1) (y^TN_i v_k+1)( v_k+1^T y ) ) [0]
    
       =
    y^T_ky 
    
         + 1/v_k+1^2( -2(y^T_k v_k+1)(v_k+1^Ty) + 1/v_k+1^2( v_k+1^T_kv_k+1)( v_k+1^Ty)^2  )
    
          + 1/v_k+1^2∑_i=1^m( v_k+1(y^TN_i v_k+1) -  1/v_k+1( v_k+1^TN_i v_k+1)( v_k+1^T y))^2 [0]
    
       =
    (y - v_k+1v_k+1^Ty/v_k+1^2 )^T_k(y - v_k+1v_k+1^Ty/v_k+1^2 )
    
          + 1/v_k+1^2∑_i=1^m( v_k+1(y^TN_i v_k+1) -  1/v_k+1( v_k+1^TN_i v_k+1)( v_k+1^T y))^2
    
       ≥ 0.

This asserts the inductive step and hence concludes the proof.



The iteration (<ref>) produces an increasing sequence of approximations 0=X̂_0 ≼X̂_1≼…≼ X.




 §.§ 2̋-optimal model reduction for symmetric state space systems


For the standard Lyapunov equation it has been shown, in <cit.>, that minimization of the energy norm induced by the Lyapunov operator (see <cit.>) is related to 2̋-optimal model reduction for linear control systems. As it turns out, a similar conclusion can be drawn for the minimization of the cost functional (<ref>) and 2̋-optimal model reduction for symmetric bilinear control systems. In this regard, let us briefly summarize the most important concepts from bilinear model reduction. Given a bilinear system Σ as in (<ref>) with dim(Σ)=n, the goal of model reduction is to construct a surrogate model Σ of the form

    Σ{ẋ(t)    = Ax(t) + ∑_i=1^m N_i x(t) w_i(t) + Bu(t)
    y(t)    = Cx(t),
    .

with A,N_i ∈^k× k, B∈^n× k, C∈^r× k and control inputs u(t)∈^r and w(t)∈^m. In particular, the reduced system should satisfy k≪ n and y(t)≈ y(t) in some norm. For the bilinear 2̋-norm defined in Definition <ref>, in  <cit.>, the authors have suggested an algorithm, BIRKA, that iteratively tries to compute a reduced model satisfying first order necessary conditions for 2̋-optimality. A corresponding pseudocode is given in Algorithm <ref>. 









To establish the connection we introduce the following generalizations of the operator ℳ:

    ℳ   ℝ^n× r→ℝ^n× r,   ℳ(X) := - A X -  X  Â -∑_i=1^mN_i X   N̂_i, 
    ℳ   ℝ^r× r→ℝ^r× r,   ℳ(X) := - Â  X -  X  Â -∑_i=1^m N̂_i X N̂_i,

where Â=V^T A V, N̂_i=V^TN_i V for i=1,…,m, and V∈^n× r is orthogonal. Our first result is concerned with the invertibility of the operators ℳ and ℳ, respectively. 


 If σ(ℳ)=-σ(+Π)⊂ℂ_+ then  σ(ℳ) ⊂ℂ_+ and σ(ℳ) ⊂ℂ_+ . 



  

  Note that σ(ℳ) is determined by the eigenvalues of the matrix 
  
    𝐌 := -I⊗ A  - Â⊗ I - ∑_i=1^m N̂_i ⊗  N_i .

  Similarly, we obtain  σ(ℳ) by computing the eigenvalues of the matrix 
  
    𝐌 := -I⊗ A - A   ⊗ I - ∑_i=1^m  N_i⊗  N_i .

Since A and N_i are assumed to be symmetric, we conclude that 𝐌=𝐌^T≻ 0. Let us then define the orthogonal matrix 𝐕=  I⊗ V. It follows that 𝐌=𝐕^T 𝐌𝐕 and, consequently, 𝐌 =𝐌^T≻ 0. A similar argument with 𝐕=V⊗ V can be applied to show the second assertion.





Given a reduced bilinear system, we naturally obtain an approximate solution to the generalized Lyapunov equation. Moreover, the error with respect to the ℳ inner product is given by the 2̋-norms of the original and reduced system, respectively.

  Let Σ denote a bilinear system (<ref>) with A=A^T≺ 0,N_i=N_i^T, i=1,…,m and B=C^T be given. Assume that ρ(^-1Π)<1. Given an orthogonal  V∈^n × k,k< n, define a reduced bilinear system Σ̂ via Â=V^T A V, N̂_i=V^TN_i V and B̂=V^TB=Ĉ^T. Let X be the solution to (<ref>), and let X̂ be the solution to ÂX̂+X̂Â+∑_i=1^m N̂_i X̂N̂_i +B̂B̂^T=0. Then 
  
    X-VX̂V^T _ℳ^2 = Σ_2̋^2  - Σ̂_2̋^2.



  By assumption it holds that the controllability Gramian X exists and the operator ℳ is invertible. Moreover, for the reduced system we obtain
  
    ℳ(X̂)    = -ÂX̂ - X̂Â - ∑_i=1^m N̂_i X̂N̂_i 
    = -V^T ( AVX̂V^T + VX̂V^T A + ∑_i=1^m N_i VX̂V^T N_i ) V 
       = V^T ℳ(VX̂V^T)V .

Hence, we obtain 

    ⟨ X-VX̂V^T , X -    VX̂ V^T ⟩_ℳ 
       = ⟨ X-VX̂V^T, ℳ(X)-ℳ(VX̂V^T) ⟩
        = ⟨ X-VX̂V^T, ℳ(X) ⟩ - ⟨ X-VX̂V^T, ℳ(VX̂V^T) ⟩
        = ⟨ X-VX̂V^T, BB^T ⟩ - ⟨ℳ(X),VX̂V^T ⟩ + ⟨ VX̂V^T, ℳ(VX̂V^T)⟩
        = (B^TXB) -(B^TVX̂V^TB) 
         - (B^T VX̂V^TB)+(X̂V^T ℳ(VX̂V^T)V) 
        = Σ _2̋^2 - 2 Σ̂_2̋^2 + (X̂ℳ(X̂) ) 
       = Σ _2̋^2 -  Σ̂_2̋^2.



Extending the results from <cit.>, we obtain a lower bound for the previous terms by means of the 2̋-norm of the error system Σ-Σ.


  Let Σ denote a bilinear system (<ref>) with A=A^T≺ 0,N_i=N_i^T, i=1,…,m and B=C^T be given. Assume that ρ(^-1Π)<1. Given an orthogonal  V∈^n × k,k< n, define a reduced bilinear system Σ via Â=V^T A V, N̂_i=V^TN_i V and B̂=V^TB=Ĉ^T. Then, it holds
  
    Σ- Σ̂_2̋^2 ≤Σ_2̋^2  - Σ̂_2̋^2,

  with equality if Σ is locally 2̋-optimal.



The proof follows by arguments similar to those used in <cit.>.
  By definition of the 2̋-norm for bilinear systems we have 
  
    Σ- Σ_2̋^2 = ([   B^T -B̂^T ] X_e [   B; -B̂ ]) ,

where X_e= [   X   Y; Y^T  X̂ ] is the solution of 

    [ A 0; 0 Â ] X_e + X_e [ A 0; 0 Â ] + ∑_i=1^m [  N_i    0;    0 N̂_i ] X_e [  N_i    0;    0 N̂_i ] + [  B; B̂ ][  B^T B̂^T ]=0.

Analyzing the block structure of X_e we thus find the equivalent expression 

    Σ- Σ̂_2̋^2    = (B^TXB-B̂^T X̂B̂ - 2 [ B^T,-B̂^T ][  Y; X̂ ]B̂ )
       = Σ_2̋^2 - Σ̂^2_2̋-2( [ B^T,-B̂^T ][  Y; X̂ ]B̂ ).

We claim that ( [ B^T,-B̂^T ][  Y; X̂ ]B̂ )≥ 0 which then shows the first assertion. In fact, Y and X̂ are the solutions of

    AY + YÂ + ∑_i=1^m N_i Y N̂_i + BB̂^T   =0 
    ÂX̂ + X̂Â + ∑_i=1^m N̂_i X̂N̂_i + B̂B̂^T   =0.

With the operators introduced in (<ref>) and (<ref>), we obtain 

    (B^T Y B̂) -   (B̂^T X̂B̂) = vec(BB̂^T)^T vec(Y) -  vec(B̂B̂^T)^T vec(X̂) 
       =  vec(BB̂^T)^T 𝐌^-1vec(BB̂^T) -  vec(B̂B̂^T)^T 𝐌^-1vec(B̂B̂^T) 
       = vec(BB̂^T)^T( 𝐌^-1 - 𝐕 ( 𝐕^T 𝐌𝐕)^-1𝐕)vec(BB̂^T).

As in <cit.>, it follows that the previous expression contains the Schur complement of 𝐌^-1 in 𝐒= [ 𝐕^T 𝐌𝐕    𝐕^T;      𝐕   𝐌^-1 ] which can be shown to be positive semidefinite. We omit the details here and instead refer to <cit.>. 

Assume now that Σ is locally 2̋-optimal. From <cit.>, we have the following first-order necessary optimality conditions

    Y^T Z + X̂Ẑ   = 0,    
      Z^T N_iY + X̂N̂_iẐ=0,   i=1,…,m, 
     
      Z^TB + ẐB̂   =0,     CY+ĈX̂ =0,

where Y,X̂ are as before and Z,Ẑ satisfy 

    A^T Z + ZÂ + ∑_i=1^m N_i^T Z N̂_i -C^T Ĉ   =0, 
    Â^T Ẑ + ẐÂ + ∑_i=1^m N̂^T_i ẐN̂_i + Ĉ^T Ĉ   =0.

From the symmetry of A,Â,N_i and N̂_i as well as the fact that B=C^T and B̂=Ĉ^T, we conclude that Ẑ=X̂ and Z=-Y. Hence, from the optimality conditions, we obtain 

    0=Z^TB+ẐB̂=-Y^TB+X̂B̂

which in particular implies that ( [ B^T,-B̂^T ][  Y; X̂ ]B̂ )= 0. This shows the second assertion.

As a consequence from Proposition <ref> and Proposition <ref>, we obtain the following result.
 

Let Σ denote a bilinear system (<ref>) with A=A^T≺ 0,N_i=N_i^T, i=1,…,m and B=C^T be given. Assume that ρ(^-1Π)<1. Given an orthogonal  V∈^n × k,k< n, define a reduced bilinear system Σ via Â=V^T A V, N̂_i=V^TN_i V and B̂=V^TB=Ĉ^T. Assume that X̂ solves ÂX̂+X̂Â+∑_i=1^m N̂_i X̂N̂_i^T +B̂B̂^T=0. If Σ is locally 2̋-optimal, then VX̂V^T is locally optimal with respect to the ℳ-norm.
 





 §.§ Equivalence of ALS and rank-1 BIRKA




So far we have shown that having a subspace producing a locally 2̋-optimal model reduction implies having a subspace for which the Galerkin approximation is locally optimal in the ℳ-norm.
In this part we, algorithmically, establish an equivalence between BIRKA and ALS. More precisely, for the symmetric case the equivalence is between BIRKA applied with the target model reduction subspace of dimension 1 for (<ref>), and ALS applied to (<ref>). The proof is based on the following lemmas.





Consider using BIRKA (Algorithm <ref>) with k=1, i.e., having initial guesses and output being vectors. Then Ã∈ is a scalar and hence we can take Λ̃= Ã and R=1 in Step <ref>. Thus B̂ = B̃, Ĉ = C̃, N̂_1 = Ñ_1, …, N̂_m = Ñ_m, and hence Steps <ref>-<ref> are redundant. Moreover, since Ṽ and W̃ are vectors, Step <ref>, is redundant.


By direct computation.


When speaking about redundant steps and operations we mean that the entities assigned in that step is exactly equal to another, existing, entity. In such a situation the algorithm can be rewritten, by simply changing the notation, in a way that skips the redundant step and still produces the same result.


Consider a symmetric generalized Lyapunov equation and let v,w∈^n be two given vectors.
Let v_birka,w_birka∈^n be the approximations obtained by applying BIRKA (Algorithm <ref>) to (<ref>) with C = B^T and initial guesses v and w. If v=w, then v_birka = w_birka.


The proof is by induction, and it suffices to show that if Ṽ = W̃ at the beginning of a loop, the same holds at the end of the loop. Thus assume Ṽ = W̃. Then Ñ_i = (W̃^T Ṽ)^-1W̃^T N_i Ṽ = Ṽ^T N_i Ṽ/V^2 = Ṽ^T N_i^T Ṽ/V^2 = N̂_i^T for i=1,…,m, and C̃ = C Ṽ = B^T W̃ = B̃^T. By Lemma <ref> we do not need to consider Steps <ref>-<ref>. We can now conclude that Step <ref> and Step <ref> are equal, and thus at the end of the iteration we still have Ṽ = W̃.




Consider a symmetric generalized Lyapunov equation and let v,w∈^n be two given vectors.
Let v_als,w_als∈^n  be the approximations obtained by applying the ALS algorithm (Algorithm <ref>) to (<ref>) with initial guesses v and w.
If v=w, then v_als = w_als.


Similar to the proof of Lemma <ref> it is enough to show that if v = w at the beginning of a loop then it also holds at the end of the loop. Hence we assume that v = w. Then Â_1 = Â_2 follows by direct calculations. 
Moreover, by assumption _k=_k^T. Thus Step <ref> and Step <ref> are equal, and hence at the end of the iteration we still have that v=w.



Consider a symmetric generalized Lyapunov equation and let v∈^n be a given vector. Let v_birka∈^n be the approximation obtained by applying BIRKA (Algorithm <ref>) to (<ref>) with C= B^T and initial guess v. Moreover, let v_als∈^n  be the approximation obtained by applying the ALS algorithm (Algorithm <ref>) to (<ref>) with initial guess v.
Then v_birka= v_als.


Firstly, Lemma <ref> and Lemma <ref> makes it reasonable to assess the algorithms with only a single initial guess as well as a single output. Moreover, Step <ref> in BIRKA as well as Steps <ref>-<ref> in ALS are redundant. Furthermore, we know from Lemma <ref> that in this situation Steps <ref>, <ref>, and <ref> of BIRKA are also redundant. Hence we need to compare the procedure consisting of Steps <ref> and <ref> from BIRKA, with the procedure consisting of Steps <ref>, <ref>, and <ref> from ALS. It can be observed that the computations are equivalent and thus the asserted equality holds if they stop after an equal amount of iterations. We hence consider the stopping criteria and note that they are the same, since 
(v^TA^Tv + v^TAv)/2v^2 = v^TAv/v^2 = Ã∈.



Theorem <ref> is applicable with ALS changed to BIRKA, using subspaces of dimension 1.









§ FIXED-POINT ITERATION AND APPROXIMATIVE ℳ-NORM MINIMIZATION

In this section we present a similar result as in the previous section, but for the fixed-point iteration. Instead of (locally) minimizing the error in the ℳ-norm with rank-1 updates, the fixed-point iteration minimizes an upper bound for the ℳ-norm, but with no rank-constraint on the minimizer.
As in the previous section, we assume that A=A^T ≺ 0, and N_i=N_i^T for i=1,…,m, as well as ρ(^-1Π)<1.

We remind ourselves that with these assumptions the symmetric generalized Lyapunov equation (<ref>) has a unique solution X, and specifically it is symmetric and positive definite, see, e.g., <cit.>.

Recall the fixed-point iteration
for the generalized Lyapunov equation (<ref>),

    (X̂_k+1) = -Π(X̂_k) - BB^T,      k=0,1,…,

with X̂_0 = 0.
The iteration has been presented as a convergent splitting in, e.g., <cit.>, <cit.>, and <cit.>.
We want to relate iteration (<ref>) to the ℳ-norm minimization problem. Hence consider the problem

    min_Δ
    Δ=Δ^T≽ 0X - X̂_k - Δ^2_ℳ.

A reason to restrict the minimization to symmetric positive semidefinite matrices is that we know that the solution X=X^T≽ 0. Hence an iteration started with X_0=0 creates a sequence of symmetric positive semidefinite approximations, ordered as X̂_k+1≽X̂_k, and converging to the symmetric positive definite solution, cf. Section <ref> and <cit.>.
Naturally Proposition <ref> gives us the solution in just one step. However, the computation is as difficult as the original problem and hence we consider minimizing an upper bound.
Similar as before we disregard the constant term X-X̂_k^2_ℳ in the minimization and consider

    min_Δ
    Δ=Δ^T≽ 0⟨Δ, Δ⟩_ℳ   - 2( Δ^T _k)
       =
    min_Δ
    Δ=Δ^T≽ 0( Δ^T(-(Δ)-Π(Δ))) - 2( Δ^T _k)
    
       =
    min_Δ
    Δ=Δ^T≽ 0( Δ^T(-(Δ)-2_k)) -( Δ^TΠ(Δ))
    
       ≤min_Δ
    Δ=Δ^T≽ 0( Δ^T(-(Δ)-2_k)),

where the inequality comes from that Δ^T and Π(Δ) are positive semidefinite matrices, and hence the trace is non-negative <cit.>.
We want to say that the last step is minimized by Δ=-^-1(_k), which is true if the residual is symmetric and positive semidefinite. We also want to show that this iteration is equivalent to the fixed-point iteration. The argument is closed by the following result.


Consider the symmetric generalized Lyapunov equation (<ref>) with the additional assumptions that A≺ 0 and ρ(^-1Π)<1. Moreover, consider the sequence of approximations constructed as
  
    X̂_0    = 0
    X̂_k+1   = X̂_k - ^-1(_k),      k = 0,1,…,
 
where _k is the residual associated with X̂_k.
Then X̂_k=X̂_k^T ≽ 0 and _k=_k^T ≽ 0, for all k≥ 0. Furthermore, the iteration (<ref>) is equivalent to (<ref>).



The first part of the proof is by induction and similar to that of Theorem <ref>. It trivially holds that X_0=X_0^T≽0 and _0=_0^T≽ 0. Now assume that this is the case for some k.
Then Δ:=-^-1(_k) is symmetric and positive semidefinite, and hence X̂_k+1 is symmetric and positive semidefinite. Moreover, since X̂_k+1 is symmetric and because of the structure of the operators in (<ref>) we know that _k+1 is symmetric. Thus what is left to show for the first part is _k+1≽0.
We have that _k+1≽0 if and only if y^T_k+1y≥0 for all y∈^n. Hence take an arbitrary y∈^n and consider

    y^T _k+1y
       =
    
    
    y^T _k y + y^T((Δ) + Π(Δ))y
    = y^T(Π(Δ))y ≥ 0,

where we in the first equality use the linearity of the operators and in the second the definition of Δ. The last inequality holds since Δ is symmetric and positive semidefinite and Π is a symmetric operator.

For the second part of the proof, what is left to show is that (<ref>) and (<ref>) are equivalent. This comes directly from adding and subtracting (X̂_k) from the right-hand-side of (<ref>), and using the linearity of .



The fixed-point iteration (<ref>) produces an increasing sequence of approximations 0=X̂_0 ≼X̂_1≼…≼ X.




























One could consider creating a subspace iteration from (<ref>), by computing a few singular vectors of ^-1(_k) and add to the basis. The method seems to have nice convergence properties per iteration in the symmetric case, but not in the non-symmetric case. However, the (naïve) computations are prohibitively expensive. See <cit.> for a computationally more efficient way of exploiting the fixed-point iteration.






§ A RESIDUAL-BASED RATIONAL KRYLOV GENERALIZATION

Given Proposition <ref> and the discussion in Sections <ref> and <ref> it seems as constructing an iteration based on the current residual can be a useful technique when solving the generalized Lyapunov equation.
Moreover, in <cit.> it is suggested that, so called, preconditioned residuals can be used in expanding the search space. It is further suggested that, for the Lyapunov equation, one such preconditioner could be a one-step-ADI preconditioner P^-1_ADI = (A-σ I)^-1(A-σ I)^-1, for a suitable choice of the shift.
This strategy can be, somewhat, motivated by Theorem <ref>.

In this section we present a method that can be seen as a generalization of the rational Krylov subspace method. We further motivate why it can be natural to see it as a generalization of the rational Krylov subspace method, by considering the standard Lyapunov equation.






 §.§ Suggested search space

For the generalized Lyapunov equation (<ref>), we suggest the following search space:

    _k := {B,(A-σ_1 I)^-1r_0,(A-σ_2 I)^-1r_1,…,(A-σ_k I)^-1r_k-1},

where r_k-1 is the most dominant left singular vector of the Galerkin residual _k-1 of _k-1, and {σ_i}_i=1^k is a sequence shift that needs to be chosen. In analogy with the discussion in <cit.>, we suggest that the shift is chosen iteratively, where the approximation error along the current direction is the largest. More precisely,

    σ_k := _σ∈[σ_min,σ_max]( r_k-1 - (A - σ I)_k-1(A_k-1 - σ I)^-1_k-1^T r_k-1),

where _k-1 is a matrix with orthogonal columns containing basis of _k-1, A_k-1 = _k-1^T A _k-1, and [σ_min,σ_max] is a search interval. Typically for a stable matrix A, then σ_min is the negative real part of the eigenvalue of A with largest real part, and correspondingly σ_max is the negative real part of the eigenvalue of A with smallest real part.
Equations (<ref>) and (<ref>) can be straightforwardly incorporated in a Galerkin method for the generalized Lyapunov equation, the pseudocode is presented in Algorithm <ref>.





In practice the computation of the left singular vector can typically be done approximatively in an iterative fashion. This would also remove the need of computing the approximative solution X̂_k in Step <ref> and the residual in Step <ref> explicitly since the matrix vector product can be implemented as _k v = A _k Y_k _k^T v + _k Y_k _k^T A^T v + ∑_i=1^m N_i _k Y_k _k^T N_i^T v + BB^Tv. Such computations may, however, introduce inexactness and present a difficulty in a convergence analysis.


The dynamic shift-search in Step <ref> can, heuristically, be changed for and analogue to <cit.>. In our case we suggest

    σ_k := _σ∈∂ S1/|𝔯_k-1(z)|,

where S is approximates the mirrored spectrum of A and  ∂ S is the boundary of S, and

    𝔯_k-1(z) := ∏_j=1^(_k-1) z-λ_j^(k-1)/∏_ℓ=1^k-1z-σ_ℓ,

with λ_j^(k-1) being the Ritz values of A_k-1. Typically S is approximated at each step using the convex hull of the Ritz values of A_k-1. It has been observed efficient in experiments since the maximization of (<ref>) is computationally faster compared to (<ref>). For a practical comparison of convergence properties, see Section <ref>.


The steps <ref>-<ref> in Algorithm <ref> can be changed for a tangential-direction approach according to <cit.>. One practical way, although a heuristic, is to do the shift search according to either (<ref>) or (<ref>), and then compute the principal direction(s) according to <cit.>, i.e., through a singular value decomposition of _k-1 - (A - σ_k I)_k-1(A_k-1 - σ_k I)^-1_k-1^T _k-1. It has been observed in experiments that such an approach tend to speed up the convergence, in terms of computation time,  since the computation of the residual is costly.


It is (sometimes) desirable to allow for complex conjugate shifts σ_k and σ_k, although, for reasons due to computations and model interpretation, one wants to keep the basis real. This is achievable also in the proposed setting by observing the following standard relation 
{(A-σ_k I)^-1r_k-1, (A-σ_kI)^-1r_k-1}
= 
{((A-σ_k I)^-1r_k-1), ((A-σ_k I)^-1r_k-1)}. 

Although it requires two shifts to be used together with the vector r_k-1.







 §.§ Analogies to the linear case

To give further motivation to the suggested subspace in (<ref>), we draw parallels with the (standard) rational Krylov subspace for the standard Lyapunov equation. The reasoning in this section can be compared to <cit.>.
We consider the equation

    (X) + bb^T = 0,

where  is defined by (<ref>) and b∈^n. Note that Definitions <ref> and <ref> are analogous for the (standard) Lyapunov equation (<ref>), but with Π=0. To prove the main result of this section we need the following lemma.




Let A∈^n× n and s_a∈ be any scalar such that (A-s_aI) is nonsingular. Moreover, let ∈^n× k, k≤ n, be orthogonal, i.e., ^T = I, and let ∈^n× n be such that ((A-s_a I)^-1)⊆().
Then = (A-s_aI)(^T A  - s_aI)^-1^T.


To prove the statement we consider the right-hand-side of the asserted equality,


    (A - s_aI)(^T A    - s_aI)^-1^T 
       =
    (A - s_aI)(^T A  - s_aI)^-1^T (A-s_aI)(A-s_aI)^-1
       = (A - s_aI)(^T A  - s_aI)^-1^T (A-s_aI)^T(A-s_aI)^-1,

where the first equality is just an insertion of an identity matrix, and the second equality follows from the assumption that ((A-s_aI)^-1)⊆(). By observing that (^T A  - s_aI)^-1^T (A-s_aI) = I the expression can be further simplified as


    (A - s_aI)(^T A  - s_aI)^-1^T    = (A - s_aI)^T(A-s_aI)^-1
       = (A - s_aI)(A-s_aI)^-1 = ,

where the second equality, again, follows from ((A-s_aI)^-1)⊆().




Let A∈^n× n, b∈^n, and let {s_i}_i=1^k+1 be a sequence of shifts such that A-s_i I is nonsingular, and define _k := {b,(A-s_1 I)^-1b,…,(A-s_kI)^-1b}, and _k+1 analogously. Let _k be a basis of _k, _k+1 a basis of _k+1, and let v_k+1∈^n be such that _k+1 = [    _k v_k+1 ]. Moreover, let _k∈^n× n be the Galerkin residual with respect to (<ref>).
Then each column of (A-s_k+1I)^-1_k is in (V_k+1), i.e., ((A-s_k+1I)^-1_k)⊆(V_k+1). Furthermore, if ((A-s_k+1I)^-1_k)⊆(V_k), then _k = 0.


We start by proving the first claim. From <cit.> we have that _k = UJU^T with

    J = [ 0 1; 1 0 ]     and      
    U = [ u_1 u_2 ],

with u_1 := V_kY_kH_k^-Te_kh_k+1,k and u_2 := v_k+1s_k+1-(I-V_kV_k^T)Av_k+1.
From the definition of u_1 one can see that u_1 = V_kα for some α∈^n, hence u_1∈_m, and thus (A-s_m+1I)^-1u_1∈_m+1. 
Similarly, from the definition of u_2 one can see that u_2 = (s_k+1I-A)v_k+1 + V_kV_k^T A v_k+1 = (s_k+1I-A)v_k+1 + V_kβ, for some β∈^n. Thus (A-s_k+1I)^-1u_2 = -v_k+1 + (A-s_k+1I)^-1V_kβ∈_k+1.
Hence ((A-s_k+1I)^-1_k)⊆_k+1, which proves the first claim.








To prove the second claim we assume that ((A-s_m+1I)^-1_m)⊆_k=(_k). Under this assumption we can use Lemma <ref> and the fact that _k = _k^T to get

    _k
       = (A - s_k+1I)_m(_k^TA_k - s_k+1I)^-1_k^T_k
       = (A - s_k+1I)_m(_k^TA_k - s_k+1I)^-1_k^T_k_k(_k^TA^T_k-s_k+1I)^-1_k^T(A^T - s_k+1I)
       = 0,

since _k is the Galerkin residual and thus _k^T_k_k=0.


The main observation and conclusion from Theorem <ref> can be phrased as follows: Consider the two spaces _k := {b,(A-s_1 I)^-1b,…,(A-s_kI)^-1b} and _k := {_-1,(A-s_1 I)^-1_0,…,(A-s_kI)^-1_k-1}, where _-1 = b and _i is the Galerkin residual in space _i, i=0,1,…,k-1 .
Then for for all relevant cases, i.e., _i≠ 0 for i=-1,0,…,k-1, we have that _k = _k. Hence the suggested subspace in (<ref>) can be seen as a natural generalization of a Rational Krylov subspace for linear matrix equations.


The result in Theorem <ref> generalizes naturally to the case B∈^n× r, where the arguments have to be done with block matrices. As noted in <cit.> the changes when going to blocks are mostly technical. Results needed to generalize <cit.> are found in, e.g., <cit.>, and the block case is implemented in the code available at Simoncini's webpage[<http://www.dm.unibo.it/ simoncin/software.html>].







§ NUMERICAL EXAMPLES

In this section we numerically compare different methods discussed in the paper. All algorithms are treated in a subspace fashion and we compare practically achieved approximation properties as a function of subspace dimension.
For small and moderate sized problems there are algorithms for computing the full solution, although costly, see <cit.>, cf. <cit.>. Nevertheless this allows for inspection of the actual relative error, i.e., 

    X-X̂_k/X,
 where X̂_k is the approximation and X is the exact solution. Moreover, it also allows comparison with the (in the Frobenius norm) optimal low-rank approximation based on the SVD.

We summarize some of the implementation details. Specifically, BIRKA is implemented as described in Algorithm <ref>, with a maximum allowed number of iterations equal to 100. Convergence tolerance is implemented as relative norm difference of the vector of sorted eigenvalues and was set to 10^-3. Each subspace is computed separately since the algorithm does not have a natural way to extend the subspace dimension based on a smaller subspace.
The complete method based on ALS is a subspace method along the lines of <cit.>, rather than an iteratively updated method as described in Theorem <ref>. Moreover, because of the structure of the generalized Lyapunov equation, the solution is symmetric even if the coefficient matrices are not. Hence we use a symmetric version of ALS even for the non-symmetric examples. The convergence tolerance is implemented analogously to BIRKA, but here it is only a scalar each time. The maximum allowed number of iterations in ALS was set to 20, and the tolerance was set to 10^-2.
Regarding rational-Krylov-type methods there are a plethora of variants to choose from. We compare the following methods, which we give short labels for the legends further down:


  * A: 𝒦_k as in (<ref>), according to Algorithm <ref>

  * B: Algorithm <ref> but with tangential directions according to Remark <ref>, though with shifts according to (<ref>)

  * C: Algorithm <ref> but with shifts according to (<ref>)

  * D: Algorithm <ref> but with tangential directions according to Remark <ref> and shifts according to (<ref>)

  * E: Standard Rational Krylov. More precisely, 𝒦_k similar to (<ref>) but instead of using the singular vector of the residual, r_k-1, we use the right-hand-side B in both (<ref>) and (<ref>)

  * F: A: 𝒦_k as in (<ref>), but with on-beforehand-prescribed shifts given as the recycling of mirrored eigenvalues from a size-10-BIRKA (convergence tolerance set to 10^-3). Mirrored eigenvalues are potentially complex, with positive real part and taken in ascending order according to the real parts.

For C, D, and F the shifts may be complex-valued, and the complex arithmetic is avoided by creating the space in accordance with Remark <ref>.
For the shift-computing versions of the rational-Krylov-type algorithms, the maximum and minimum eigenvalues  were computed and the mirrored shift-search-boundaries were perturbed with the factors 1.01 and 0.99 for the upper and lower boundaries respectively, as to slightly enlarge the region.
Orthogonalization of the basis is implemented very similar to how the SVD-based orthogonalization is done in . More precisely, it is using  built-in QR factorization and keep vectors only if the corresponding diagonal element in R is large enough.
Implementation for the methods A-F is available online.[<https://people.kth.se/ eringh/software/res_rat_Kry_type/>] 
The simulations were done in  R2018a (9.4.0.813654) on a computer with four 1.6 GHz processors and 16 GB of RAM.

We test the algorithms on three different problems. In all problems we are computing an approximation to an associated controllability Gramian to a bilinear control system, as in (<ref>). The examples all have stable Lyapunov operators. The first example is symmetric, the second is non-symmetric but symmetrizable, and the third example non-symmetric. 




 §.§ Heat equation


The first example is motivated by an optimal control problem for the selective cooling of steel profiles, see <cit.>. In this example, the state variable y models the evolution of a temperature and is described by a two-dimensional heat equation,

    ∂/∂ t w(x,y,t)    = Δ w(x,y,t)              (x,y,t)∈(0,1)×(0,1)×(0,T)
    
    w(x,y,0)   =w_0(x,y)        x∈(0,1)
    
    -∂/∂ x w(0,y,t)   = 0.5 (w(0,y,t) - 1) u(t)        (y,t) ∈ (0,1)×(0,T)
    
    w(1,y,t)   =0        (y,t) ∈ (0,1)×(0,T)
    
    w(x,0,t)   =0        (x,t) ∈ (0,1)×(0,T)
    
    w(x,1,t)   =0        (x,t) ∈ (0,1)×(0,T),

where a control u(t) enters bilinearly  from the left boundary through a Robin boundary condition, and the other spatial boundaries satisfy homogeneous Dirichlet conditions. The control can be interpreted as the spraying intensity of a cooling fluid. The system is discretized in space using centred finite difference, which yields a bilinear system with
A∈^5041× 5041, B∈^5041, m=1, and N_1=N∈^5041× 5041. 

It can be further noted that, A=A^T≺0 and N=N^T, and hence the theory of 2̋-optimality and the definition of the ℳ-norm, (<ref>), from above is applicable.



We compare different methods discussed in the paper, and we compare both the relative residual norm, as well as the relative error. For readability the plots have been split, hence in Figure <ref> we compare across different methods, and in Figure <ref> we compare between different flavours the rational-Krylov-type methods.
It can be observed, see Figure <ref>, that for this problem BIRKA has extremely good performance, even outperforming the SVD in relative residual norm. Nevertheless, the larger BIRKA subspaces are rather costly to compute. In comparison ALS shows good performance compared to the rational-Krylov-type subspace, and is rather cheap to compute.
When comparing the different Rational-Krylov-type methods, see Figure <ref>, we see that standard rational-Krylov (E) has the problem that it reaches an invariant subspace and is unable to extend larger than dimension 13. However, stagnation has already been observed. The methods A, C, and F are similar, although B and D are only slightly worse in the error per subspace dimension comparison and practically sometimes even faster to compute.

Since the ℳ-norm is defined for this problem we compare the relative error also in this norm, see Figure <ref>. The trend is similar as in the Frobenius norm, although it can be observed that in general the error is smaller and BIRKA has best performance, even compared to the SVD (consider that this is relative error, however, in another norm).




 §.§ 1D Fokker-Planck

The second example is from quantum physics, where a one-dimensional Fokker-Planck equation describes the evolution of a probability density function, ρ, of a particle affected by a potential. Parts of the potential can be manipulated by a so-called optical tweezer, which constitutes our control. For more details see <cit.>.
More precisely, we consider

    ∂/∂ tρ(x,t)     = ν∂^2/∂ x^2ρ(x,t)
    + ∂/∂ x(ρ(x,t)∂/∂ xV(x,t)) 
            (x,t)∈(-6,6)×(0,T)
    ρ(x,0)   =ρ_0(x)        x∈(-6,6)
    ν∂/∂ xρ(x,t)   =-ρ(x,t)∂/∂ xV(x,t)       (x, t) ∈{-6,6}× (0,T),

where the potential is V(x,t) = W(x) + α(x) u(t), with the ground (fixed) potential being W(x) = (((0.5x^2-15)x^2 + 199)x^2 + 28x + 50)/200, and α(x) is an approximately linear control shape function, see <cit.> for further details. In a weighted inner product, the dynamics can be described by self-adjoint operators. However, here we employ an upwinding type finite difference scheme with 5000 grid point, leading to a non-symmetric system. As has been pointed out in <cit.>, the system matrix A is not asymptotically stable due to a simple zero eigenvalue associated with the stationary probability distribution. Using a projection-based decoupling, it is however possible to work with an asymptotically stable system of dimension n=4999.  Similar to the first example, our control variable is a scalar and, consequently, we only obtain a single bilinear coupling matrix N_1=N.  Since the system is non-symmetric, the operator ℳ is generally indefinite and hence we make no comparisons in the ℳ-norm.



The plots in Figures <ref> and <ref> are analogous to the plots in Figures <ref> and <ref> respectively. However, for this example the direct solver stagnated at a relative residual of about 10^-8, which can be seen in the stagnation of the SVD approximation in the left of Figure <ref>. As a result, the comparisons of relative error performance, the right of Figures <ref> and <ref>, has an artificial stagnation around 10^-8. At this level the convergence stagnates since it measures the discrepancy between the method-approximations and the inexact reference solution, rather than the true error of the method-approximations. Nevertheless we believe the comparisons to be fair more or less upto to point of stagnation, which is further justified by the relative residual plots showing similar behaviour. However, the relative residual also indicating stagnation around 10^-8 for the other methods as well, although not quite as clear as for the SVD.

From Figure <ref> we see the BIRKA performs well for this problem. However, the subspaces of dimension 28 and 29 did not converge in a 100 iterations and hence these, as well as iteration 30, are left out of the plots. This illustrates a drawback with the method. The performance difference between ALS and the rational-Krylov-type method is slightly smaller compared to the previous example. Among the rational-Krylov-type methods A, B, and F seems to have similar performance, whereas C is clearly worse. Method E has the same problem as in the previous example, and method D also ends up with an insufficient subspace.












 §.§ Burgers equation


In the third example we consider an approximation to the one-dimensional viscous Burgers equation

    ∂/∂ t w(x,t) + w(x,t) ∂/∂ x w(x,t)    = ν∂^2/∂ x^2w(x,t)               (x,t)∈(0,1)×(0,T)
    
    w(x,0)   =w_0(x)        x∈(0,1)
    
    w(0,t)   =u(t)        t ∈ (0,T)
    
    w(1,t)   =0        t ∈ (0,T),





where ν=0.1 is constant. The solution w(x,t) can be interpreted as a velocity and the equation occurs in, e.g., modelling of gas or traffic flow. A control input u(t) is applied at the left boundary. The problem is discretized in space using centered finite differences with 71 uniformly distributed grid points. Using a second order Carleman bilinearization, we obtain a bilinear control system approximation with A,N∈^5112× 5112 and B∈^5112, see <cit.> for further details. 

Note that in this case A is an asymptotically stable but non-symmetric matrix. In order to ensure the positive semidefiniteness of the Gramian, we scale the control matrices N and B with a factor α=0.25. We emphasize that the control law is scaled proportionally with 1/α such that the dynamics remain unchanged, for further discussion see <cit.>.

The comparison is similar to the previous examples and the Figures <ref> and <ref> are analogous to the respective Figures <ref> and <ref>. The problem is difficult in the sense that the singular values of the solution decays slowly. Moreover, the direct method stagnates at a relative residual norm of 5·10^-6. This is, however, less visible compared to the previous example since in general the convergence is slower.

For this example the performance of BIRKA is not significantly better than other methods, which is not surprising since the theoretical justifications for the method are not valid. ALS shows faster convergence in relative residual norm but slower convergence in relative error, as well as indications of stagnation. However, the theoretical justifications for ALS are also not valid for this problem and the result is in line with the results in <cit.>. Regarding the rational-Krylov-type methods it seems as if method D and B has good performance, and method E is not working for this problem either.






§ CONCLUSIONS AND OUTLOOKS

In this paper we have studied iterative methods for computing approximations to the generalized Lyapunov equation. 
The methods have been studied from an energy-norm perspective, as well as a model-reduction perspective, with connections made in between. 
Common for all methods studied is that they use the current residual in the iterations. Computing the residual can in itself be costly for a truly large scale problem, although approximate dominant directions can be computed in an iterative fashion, resulting in an inner-outer-type iteration. However, more research is needed to understand the consequences of such inexact subspaces. Moreover, we have proposed a rational-Krylov-type subspace for solving the generalized Lyapunov equation. Simulations indicate competitive performance, at least in the non-symmetric case where optimality statements for the other methods are no longer valid. Simulations show that methods A and F do good or decently good for all three examples. The ALS iteration, as well as results from the literature, cf. <cit.>, 
seems to indicate that subspaces of the type (A-σ I -μ N)^-1B could be useful. Although we have not been able to exploit this efficiently. More research is needed to understand the theoretical aspects of the suggested, and related, spaces.









§ ACKNOWLEDGMENT

This research started when the second author visited the first author at the Karl-Franzens-Universität in Graz; the kind hospitality was greatly appreciated. The visit was made possible due to the generous support from the European Model Reduction Network (COST action TD1307, STSM grant 38025).










spmpsci

















