`plyapc` returns a non-standard factorization

[JurajLieskovsky]

I am a bit perplexed by the design choice of plyapc computing the upper triangular fractor U of the solution X = U U'.

Cholesky factorization typically computes the lower triangular factor L of matrix A such that A = L L'. The lyapchol function in matlab, which is according to the documentation uses the SLICOT SB03OD routine, computes the upper triangular factor U, where X = U' U. I am not sure if the function is just a wrapper or a re-implementation but I assume the original SLICOT routine also uses the same factorization of the result. The factorization is then of course consistent with the typical Cholesky factorization as L = U'.

So far I haven't run into any issues with this, but plyapc is definitely the odd man out. Additionally if plyapc used the typical U' U' or L L', it could directly return a Cholesky type as a result.

[JurajLieskovsky]

Here's a bit of matlab code if anyone wants to double check

A = [-1 0 0; 0 -2 0; 0 0 -3];
B = [1 2 3; 2 3 1; 3 1 2];

X = lyap(A, B * B')

U = lyapchol(A, B)
X1 = U' * U

[andreasvarga]

I am a bit perplexed by the design choice of plyapc computing the upper triangular fractor U of the solution X = U U'.

plyapc computes X = U*U' to solve AX + XA' + BB' = 0 and X = U'*U to solve A'X + XA + B'B = 0.

The lyapchol function in matlab, which is according to the documentation uses the SLICOT SB03OD routine, computes the upper triangular factor U, where X = U' U. I am not sure if the function is just a wrapper or a re-implementation but I assume the original SLICOT routine also uses the same factorization of the result. The factorization is then of course consistent with the typical Cholesky factorization as L = U'.

This is not the case. plyapc uses the same interface as SB03OD, and the MATLAB interface is the odd man out.

Since you were perplexed by this choice, I will try to explain one of the reasons for it. The main application of this function is in order reduction of linear systems, where the two gramians for controllability P and observability Q satisfy AP + PA' + BB' = 0 and A'Q + QA + C'C = 0 and are computed as P = S*S' and Q = R'*R (see for example the gbalmr). The Hankel singular values, exceeding a certain threshold, are used to determine the reduced system order and are computed as the singular values of R*S, which is upper triangular provided S and R are upper triangular. This nice property is even preserved in the case of descriptor system representations, where the SVD computation involves a product R*E*S, with E upper triangular. In MATLAB you need to use R'*S, which is a full matrix.

[JurajLieskovsky]

Thank you for the very insightful explanation.