7
$\begingroup$

I solve the Lyapunov equations :

$$ A W_C E^T + E W_C A^T + B B^T = 0 $$ $$ A^T W_O E^T + E W_O A + C^T C = 0 $$

to obtain $ W_C $ and $W_O$. My aim is to get the left and right eigenspaces of $W_C W_O$. This is basically for model order reduction through balanced truncation.

Currently I am explicitly generating $W_C$ and $W_O$, and calculating the schur decomposition as

W_c = lyap(A,B*B',[],E);
W_o = lyap(A',C'*C,[],E);
W_j = W_c*W_o;
[Vs,Ts] = schur(W_j);
Vst = Vs';
V_lk = Vst(:,1:ordr);
V_rk = Vs(:,1:ordr);

How can I optimize the process of calculating the dominant eigenspaces? Links to some theory behind any techniques will be very helpful.

Additional Information : I am aware of ADI, and read quite a few of the papers based on it. However, there are 2 problems with it being applicable to my case.

  1. My original system is dense, a product of another model reduction. Isn't ADI specific to Sparse systems?
  2. My grammians will NOT be positive definite, because its not fully controllable/observable. They are in-fact indefinite (negative and positive) eigenvalues. Therefore Cholesky factor is not applicable. Is there a different factorization I could use with ADI?
$\endgroup$
  • $\begingroup$ To clarify, by "dominant" you mean the eigenvalues of largest magnitude? And you only want a few of the biggest ones, for which computing the Schur decomposition is overkill? $\endgroup$ – Daniel Shapero Jun 4 '13 at 14:08
  • $\begingroup$ Yes. To be exact, I want to obtain the vectors representing spaces spanned by the eigenvectors corresponding to some k largest eigenvalues. If not Schur, how best? Can you point me to some smarter methods? Thanks. $\endgroup$ – Milind R Jun 8 '13 at 11:41
  • 1
    $\begingroup$ The input to this reduction is a system originally of order ~10k, reduced to 200 by a different method. $\endgroup$ – Milind R Jun 17 '13 at 16:35
  • 3
    $\begingroup$ Regarding (2) @FedericoPoloni is right. Semi-definiteness is even welcomed in ADI, as one looks for low-rank factors. However, I am not sure about indefiniteness. Peter Benner does not give conditions for convergence but refers the reader to earlier papers by Wachspress on ADI. $\endgroup$ – Jan Jun 19 '13 at 7:25
  • 1
    $\begingroup$ Maybe I am making a silly mistake, but this should be a proof that ADI works also if the Gramian is semidefinite: (1) use Kalman decomposition $A=\begin{bmatrix}A_{11} & A_{12}\\\\ 0 & A_{22}\end{bmatrix}$, $B=\begin{bmatrix}B_1\\\\0\end{bmatrix}$; (2) notice that the second block never gets altered, so you could as well work on $(A_{11}, B_1)$ (3) This pair is controllable, so ADI works. As for indefinite, I don't know. $\endgroup$ – Federico Poloni Jun 19 '13 at 7:48
5
$\begingroup$

One of the best-performing methods for solving large-scale Lyapunov equations is ADI. It is an iterative algorithm that returns an approximate low-rank decomposition $X \approx VV^T$ of the solution $X$. In this case, you can work with this decomposition of both Gramians to reduce the eigenproblem to a smaller one.

I suggest you to start approaching this algorithm by reading the paper Numerical solution of large-scale Lyapunov equations, Riccati equations, and linear-quadratic optimal control problems, by Benner, Li and Penzl, the user manual for the Matlab library Lyapack, or one of the many talks on the webpage of Peter Benner.

$\endgroup$
  • $\begingroup$ @FredericoPoloni: Thanks, I have added more info to the question. $\endgroup$ – Milind R Jun 16 '13 at 4:20
  • $\begingroup$ Can you join chat.stackexchange.com/rooms/9436/lyapunov ? I have a few questions which would be difficult to ask in this format. $\endgroup$ – Milind R Jun 28 '13 at 8:38
0
$\begingroup$

In a 2012 paper, Simoncini, Szyld and Monsalve review several techniques for solving large scale Riccati equations but mainly focus on a Galerkin projection approach to Lyapunov equations.

However, as a basic assumption they have $A$ always Hurwitz, what would render the Gramians definite (if E is positive definite).

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.