Decomposing a positive semi-definite matrix

[xponential]

Homework Statement

if M is dxd positive semi-definite matrix, then:

M = L'(or L transpose) * L

where L is a matrix of dimensions rxd

How can I get L given M and the number of rows, r?

for example,

if M is 800x800 p.s.d, I want L such that it is a 30x800 matrix and so
M (800x800) = L' (800x30) * L (30x800)

If you know a MATLAB function that can do just that, it would be perfect!

 

[Ray Vickson]

Homework Statement

if M is dxd positive semi-definite matrix, then:

M = L'(or L transpose) * L

where L is a matrix of dimensions rxd

How can I get L given M and the number of rows, r?

for example,

if M is 800x800 p.s.d, I want L such that it is a 30x800 matrix and so
M (800x800) = L' (800x30) * L (30x800)

If you know a MATLAB function that can do just that, it would be perfect!

If M is an nxn positive semi-definite matrix, then one can show that there exist an upper-triangular matrix U such that M = U^T U. One obtains U by the algorithm for Cholesky factorization. Although many discussions of that algorithm assume M is positive-definite, if one is careful to put the argument in the correct form one can prove that it also applies to the positive-semidefnite case, with very slight modifications. I don't have access to Matlab, but the latest versions of Maple can do Cholesky factorization for positive-semidefinite matrices. Basically, the algorithm for a psd matrix M = (m_ij) is:
(1) If m₁₁ = 0, all elements in row 1 and column 1 must also be zero (otherwise, M could not be psd). If that holds, ignore row and column 1 and start again. Therefore, we may as well assume m₁₁ > 0. Set
$$u_{11}= \sqrt{m_{11}} \text{ and } u_{1j} = \frac{m_{1j}}{u_{11}},\; j=2,\ldots, n$$
(2) For i ≥ 2: let
$$u_{ii}^2 = m_{ii} - \sum_{p=1}^{i-1} u_{pi}^2.$$ If $u_{ii}^2 = 0$ we must have $m_{ij} - \sum_{p=1}^{i-1} u_{pi} u_{pj} = 0 \; \forall j > i.$ In that case, set u_ij = 0 for all j ≥ i. Otherwise, set
$$u_{ii} = \sqrt{u_{ii}^2}, \; u_{ij} = \frac{m_{ij} - \sum_{p=1}^{i-1} u_{pi} u_{pj}}{u_{ii}} \; \forall j > i.$$

So, in the psd case some rows of U will be all zeros (unlike in the pd case); we can omit the zero rows of U to get a new matrix V which is no longer upper-triangular, but still gives
M = V^T V.

This algorithm is so simple that I used to have it programmed into an HP hand-held calculator; it could decompose matrices up to about 10x10.

RGV

 

[xponential]

Thanks Ray .. I actually know about using Cholesky's in MATLAB .. however .. as in the example I provided .. I am required to find L such that it is (30x800) matrix instead of (800x800).

 

[Ray Vickson]

Have you TRIED cholesky? Perhaps U has 770 zero rows, which can be dropped. (I don't know if this is true; you would need to try it to see.)

Anyway, does your course notes or textbook say nothing about this?

RGV