Conditioning of matrix factorizations and square roots

Question

For my application, I need factors $\tilde C$ so that

$$ \tilde C{}^T \tilde C = C^TMC $$ where $C$ is a long matrix, i.e. $C$ has much more columns than rows, and $M$ is a small symmetric positive definite (mass) matrix. Since my code assembles $MC$, I will also have to use some inverses of $M$.

I have found that with the Cholesky factorization $M=L^TL$ and $\tilde C = LC$ the relative difference in the Frobenius norm $$ \|LM^{-1}MC - L^{-T}MC\|~/~\|MC\| $$ is of order $10^{-4}$.

Question: Is there a better factorization for this purpose?

I basically use the scipy.inv and scipy.cholesky routines of Python's scipy module.

Answer 1

If I understand your question correctly you are looking for a factorisation $$\tilde{C}^T \tilde{C} = C^TMC,$$ where $M$ is symmetric positive definite. Moreover you have at hand $M$ and $B=MC$.

For $M$ you can compute Cholesky, LDLT, and eigenvalue decompositions: $$M = L_1 L_1^T = L_2 D_2 L_2^T = Q \Lambda Q^T,$$ where $L$ are lower triangular, $Q$ orthogonal, $D_2$ and $\Lambda$ diagonal and positive.

Accordingly you can define and compute $\tilde C$ as

1. $\tilde C = L_1^TC$, to be computed as $\tilde C = L_1^{-1} B$, where $L_1^{-1}\cdot$ is a back-substitution on the columns of $B$

2. $\tilde C = D^{\frac 12} L_2^TC$, to be computed as $\tilde C = D^{-\frac12}L_2^{-1} B$. Here we have back-substitution followed by row-scaling.

3. $\tilde C = Q\Lambda^{\frac12}Q^T C$, to be computed as $\tilde{C} = Q\Lambda^{-\frac12}Q^T B$. Different approches are possible, but maybe computing $M^{-\frac12}B$, with $M^{-\frac12} = Q\Lambda^{-\frac12}Q^T$ is advisable.

It is hard to guess which approach is best, without more information on $M$ and it's properties, but I would first try approach 2., which could be more stable, for ill-conditioned $M$.