# Factorize Positive Definite Symmetric Matrix

Let's start from the assumption of disposal of a positive definite symmetric matrix of size $\ (N,N)$. For some reason I have to factorize this matrix: I am already aware of the Cholesky factorization method, but I am wondering if there is any procedure allowing for matrix factorization based on eigen-vectors and eigen-values.

I am not a mathematician, therefore my background knowledge on the topic is rather scarce and little. I need the factorization within a MATLAB code, and the Cholesky factorization is quite time-expensive, so that I'd more than appreciate moving towards a more practical (but effective) way to compute the lower triangular matrix.

Any support will be more than welcome.

• Since you said symmetric, I assume your matrix is real. Every positive definite symmetric matrix is diagonalizable in an orthonormal basis. Thus there exists $P$ invertible (and actually, $P$ can be taken to be an orthogonal matrix) such that $A=PDP^{-1}$ where $D$ is diagonal and the diagonal coefficients of $D$ are the eigenvalues of $A$, which are all positive. The same is true for positive definite complex matrices (they are not symmetric, they are self-adjoint or hermitian). Just replace orthogonal by unitary. – Julien Mar 5 '13 at 13:13
• "Cholesky factorization is quite time-expensive" Surely not more expensive than computing eigenvalues? – leonbloy Mar 5 '13 at 13:39
• @leonbloy ...unless for some reason the eigenvalues and/or eigenvectors are already computed. Could be. – adam W Mar 5 '13 at 13:53
• Is LU slow as well? One difference between LU and Cholesky is that scalar square roots are needed in the Cholesky. These would also be needed if using eigenvalues. – adam W Mar 5 '13 at 16:00
• LU is faster, but now I will have to account for the square roots as you suggest. – fpe Mar 5 '13 at 16:48

Edit: Another decomposition which involves an upper triangular matrix (therefore a lower triangular one if you apply it to the adjoint first) is the QR decomposition, which is based on Gram-Schmidt. Of course, there is also the general $LU$ decomposition, but Cholesky is twice as efficient. I don't know how the efficiency of QR compares to the one of Cholesky.

Now here is my previous answer, which did not answer the question...

Since you said symmetric, I assume your matrix is real.

Every positive definite symmetric (real, square) matrix is diagonalizable in an orthonormal basis.

Thus there exists $P$ invertible (and actually, $P$ can be taken to be an orthogonal matrix, wich means $P^tP=PP^t=I_n$) such that $$A=PDP^{-1}=PDP^t$$ where $D$ is diagonal and the diagonal coefficients of $D$ are the eigenvalues of $A$, which are all positive.

Note that $P$ is simply a matrix whose columns constitute an orthonormal basis made of eigenvectors.

The same is true for positive definite complex matrices (they are not necessarily symmetric, they are self-adjoint or hermitian). Just replace orthogonal by unitary (ie $P^*P=PP^*=I_n$) and $P^t$ by $P^*$.

• I am supposing then you meant to actually specify that $A = (PD^{\frac{1}{2}}P^t)(PD^{\frac{1}{2}}P^t)^t$? OP seems to have gotten the point. – adam W Mar 5 '13 at 14:00
• yes I already figured out what @julien meant – fpe Mar 5 '13 at 14:01
• @fpe Actually, I apologize: I read your question too quickly and thought you wanted any good factorization. The diagonalization does not yield the Cholesky decomposition, since nothing is upper or lower triangular here. Sorry about that. – Julien Mar 5 '13 at 14:09
• @adamW Actually, I had forgotten about Cholesky. ANd this does not help for Cholesky... – Julien Mar 5 '13 at 14:10
• @julien: you are rigt, your method does not lead to any lower or upper triangular matrix. I do have need for a lower triangular matrix, but without Cholesky – fpe Mar 5 '13 at 14:15

I wanted to comment on the LU stuff above, but i'm pretty new around here and i didn't know how.

LU is supposed to be less efficient than Choleski, by a factor 2.

As for the square root problem with Choleski, it is possible to compute $D$ and $U$ such that $R=\sqrt{D}U$ and $A=R^tR$ with a modified version of Choleski :

Input : A symmetric definite positive matrix $A\in \mathbb R^{n\times n}$.
Output : A triangular matrix $R$ with positive diagonal coefficients such that $A=R^T\cdot R$.

$U:=A$;
For{$\kappa := 1$ To $n$}{
$~~$For{$i:=1$ To $\kappa-1$}{
$~~~~$For{$j:=1$ To $i-1$}{
$~~~~~~$$u_{i,\kappa} := u_{i,\kappa} - \mu_{j,i} u_{j,\kappa}; ~~~~} ~~~~$$\mu_{i,\kappa} := \frac{u_{i,\kappa}}{u_{i,i}}$;
$~~$}
$~~$$u_{\kappa,\kappa} :=u_{\kappa,\kappa} -\sum_{i=1}^{\kappa-1} \mu_{i,\kappa} u_{i,\kappa}; ~~$$d_{\kappa} := \frac{1}{u_{\kappa,\kappa}}$;
}
Return $R=\sqrt{diag\left(\overrightarrow{d}\right)}\cdot Trig(U)$.

• What is $\overrightarrow{b}_{\kappa}$? – adam W Mar 5 '13 at 19:19
• a mistake on my end. I adapted an algorithm I wrote a long time ago to compute QR using choleski. Fixed now – imj Mar 5 '13 at 20:02
• Okay, I will code this later and test it, it looks like something I want to test out, but right now I wonder what is $Trig(U)$? Maybe when I code it up I will realize, but right now the only thing I can think of is some sort of scaling to normalize $U$. – adam W Mar 5 '13 at 23:17
• Trig(U) is the upper right part of U. Otherwise, since none of the lower left coefficients of $U$ are read/overwritten by the algorithm (they're redundant since the matrix is symetric), you'd end up with something that is not triangular :/ This truncation can be done at the beginning of the algorithm as well. – imj Mar 5 '13 at 23:50
• The algorithm is described in Accuracy and Stability of Numerical Algorithms (maths.manchester.ac.uk/~higham/asna/index.php), p197. The sqrt-free variant only gets a mention, I don't remember if I found an algorithm elsewhere or if I made it myself. – imj Mar 6 '13 at 0:06

Edit : Did not see the edit above, my answer provides comparaison between Choleski and QR:

Once you have $A= ( P D^{1/2} P^t) (P D^{1/2} P^t)^t$ you can perform $QR$ decomposition (using Householder for instance) on $B=(PD^{1/2}P^t)^t=QR$ Then since $A=B^tB$, you get :

$$A = R^tQ^tQR = R^tR$$ and thus your factorisation.

However if memory serves, Householder is more expensive than Choleski, and there isn't any (known) way to compute $QR$ faster than $4n^3/3$ operations (With Givens, Householder or Choleski), whereas Choleski costs $n^3/3$. So even if that technically works, I'm afraid it won't be of much help.

PS : using Gram Schmidt to compute QR is easier, but awful in terms of complexity, and numerical stability : $2n^3$ operations. The fact that computing $B^TB$ and using Choleski is one of the fastest way to compute $QR$ (but not great numerically) is a tell that QR probably isn't the answer.