1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Yes, you are right. Let me correct this. $\endgroup$ – Jan Jul 9 '14 at 8:08
1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ Just judging from some recent experience, the Cholesky factorization (1.) does a much better then the eigenvalue decomposition (3.). $\endgroup$ – Jan Sep 17 '14 at 7:06

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.