1
$\begingroup$

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))

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!!

  • 0
    You will run into complex numbers as matrix entries of $L$ when you implement the algorithm without the checking condition.2011-01-09

1 Answers 1

2

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.

  • 1
    matlab probably uses one of the standard libraries like LAPACK. The "unblocked" version of the code (http://www.netlib.org/lapack/double/dpotf2.f) seems pretty similar to your code. The optimized routine (http://www.netlib.org/lapack/double/dpotrf.f) is more complicated.2011-01-09