1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Homework Help: Decomposing a positive semi-definite matrix

  1. Sep 28, 2012 #1
    1. The problem statement, all variables and given/known data

    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!
  2. jcsd
  3. Sep 29, 2012 #2

    Ray Vickson

    User Avatar
    Science Advisor
    Homework Helper

    If M is an nxn positive semi-definite matrix, then one can show that there exist an upper-triangular matrix U such that M = UT 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 = (mij) is:
    (1) If m11 = 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 m11 > 0. Set
    [tex] u_{11}= \sqrt{m_{11}} \text{ and } u_{1j} = \frac{m_{1j}}{u_{11}},\; j=2,\ldots, n[/tex]
    (2) For i ≥ 2: let
    [tex] u_{ii}^2 = m_{ii} - \sum_{p=1}^{i-1} u_{pi}^2.[/tex] 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 uij = 0 for all j ≥ i. Otherwise, set
    [tex] 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.[/tex]

    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 = VT 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.

    Last edited: Sep 29, 2012
  4. Sep 29, 2012 #3
    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).
  5. Sep 29, 2012 #4

    Ray Vickson

    User Avatar
    Science Advisor
    Homework Helper

    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?

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook