Left and right eigenspaces of the product of Gramians

Question

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?

Answer 1

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.

Answer 2

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).