Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

We have been told in our homework to implement the following program:

function A = Cholesky(A) % Cholesky Decomposition of a Matrix A, so that the resulting matrix L gives A=L*L

N = size(A,1);

for n=1:N
if abs(A(n,n)) error('* ERROR * Cholesky Decomp doesn't exist'); end; A(n:N,n) = A(n:N,n) - A(n:N,1:n-1) * A(n,1:n-1)';
A(n:N,n) = A(n:N,n) / sqrt(A(n,n));
end;

Could anyone please explain what the line

if abs(A(n,n))<eps

means. I can sort of see that it is saying that the absolute values of the diagonal values can't be zero however what impact would this have on my decomposition. Thanks!!

share|improve this question
    
You will run into complex numbers as matrix entries of $L$ when you implement the algorithm without the checking condition. –  user17762 Jan 9 '11 at 4:34

1 Answer 1

up vote 1 down vote accepted

In one of the following lines, one divides by $A_{n,n}$. So $A_{n,n}$ must be non-zero. Now, computations in the computer are not accurate so that the correct way to test $A_{n,n} = 0$ is $|A_{n,n}| < \epsilon$. In fact, when $0<|A_{n,n}| < \epsilon$ the entire procedure will be very sensitive to the exact value of $A_{n,n}$, and so the result of the computation won't be trustworthy.

The Cholesky decomposition is usually defined only for (symmetric) positive-definite matrices, since in that case (contrary to the non-negative case) the decomposition is unique. The algorithm you describe uses some form of Gaussian elimination, so that $A_{n,n}$ does stand for an eignevalue.

share|improve this answer
    
thanks @Yuval Filmus Does the second part of your answer mean this is not the usual way to do a Cholesky decomposition in matlab? –  Lucy Marshall Jan 9 '11 at 4:40
1  
matlab probably uses one of the standard libraries like LAPACK. The "unblocked" version of the code (netlib.org/lapack/double/dpotf2.f) seems pretty similar to your code. The optimized routine (netlib.org/lapack/double/dpotrf.f) is more complicated. –  Yuval Filmus Jan 9 '11 at 6:00

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.