Cholesky

The Cholesky decomposition can be implemented in a number of different ways.

Outer Product Cholesky

OuterProductCholesky[A_]:= Module[
{m=Length[A],L=IdentityMatrix[m],B=A,a,b},
Do[
a=(Sqrt[B[[i,i]]]);
b=B[[i,i+1;;m]]/a;
B[[i+1;;m,i+1;;m]] -=Outer[Times,b,b];
L[[i,i]]=a;
L[[i+1;;m,i]]=b,
{i,1,m}];
L
]

m = 6;
A = RandomReal[{-1,1},{m,m}] + m IdentityMatrix[m];
A=A+Transpose[A];
L=OuterProductCholesky[A];
MatrixForm[L]
Norm[A-L.Transpose[L]]

Loopy Cholesky

LoopyCholesky[A_]:= Module[
{m=Length[A],L=ConstantArray[0,{m,m}],i,j},
Do[
Do[
L[[i,j]]=(1/(L[[j,j]]))(A[[i,j]]-Sum[L[[i,k]]L[[j,k]],{k,1,i-1}]),
{j,1,i-1}];
L[[i,i]]=(Sqrt[A[[i,i]]-Sum[((L[[i,k]])^2),{k,1,i-1}]]),
{i,1,m}];
L
]

m = 12;
A = RandomReal[{-1,1},{m,m}] + m IdentityMatrix[m];
A=A+Transpose[A];
L=LoopyCholesky[A];
Norm[A-L.Transpose[L]]

Time Trials

m = 12;
A = RandomReal[{-1,1},{m,m}] + m IdentityMatrix[m];
A=A+Transpose[A];
Timing[LoopyL=LoopyCholesky[A];]

Timing[OuterProductL=OuterProductCholesky[A];]
Norm[LoopyL-OuterProductL]

Built-In Commands

Built in commands may do more than you think.  They may also store stuff in compressed formats.

?LUDecomposition

A={
{1.0,2.0,3.0},
{3.,4.0,5.0},
{45.0,7.0,9.0}
};
{lu,p,c}=LUDecomposition[A]

m=3;
L=LowerTriangularize[lu,-1]+IdentityMatrix[m];
U=UpperTriangularize[lu];
L.U-A

?QRDecomposition

{Q,R}=QRDecomposition[A];
Norm[A-Transpose[Q].R]
Norm[Q.A-R]

Lots of other libraries do other things.  You always need to read the documentation!

gsl library web page
14.2 QR Decomposition

A general rectangular M-by-N matrix A has a QR decomposition into the product of an orthogonal M-by-M square matrix Q (where Q^T Q = I) and an M-by-N right-triangular matrix R,

     A = Q R

This decomposition can be used to convert the linear system A x = b into the triangular system R x = Q^T b, which can be solved by back-substitution. Another use of the QR decomposition is to compute an orthonormal basis for a set of vectors. The first N columns of Q form an orthonormal basis for the range of A, ran(A), when A has full column rank.
— Function: int gsl_linalg_QR_decomp (gsl_matrix * A, gsl_vector * tau)

    This function factorizes the M-by-N matrix A into the QR decomposition A = Q R. On output the diagonal and upper triangular part of the input matrix contain the matrix R. The vector tau and the columns of the lower triangular part of the matrix A contain the Householder coefficients and Householder vectors which encode the orthogonal matrix Q. The vector tau must be of length k=\min(M,N). The matrix Q is related to these components by, Q = Q_k ... Q_2 Q_1 where Q_i = I - \tau_i v_i v_i^T and v_i is the Householder vector v_i = (0,...,1,A(i+1,i),A(i+2,i),...,A(m,i)). This is the same storage scheme as used by lapack.

    The algorithm used to perform the decomposition is Householder QR (Golub & Van Loan, Matrix Computations, Algorithm 5.2.1).

?CholeskyDecomposition

SVD Images;

Pic=Import["/home/faculty/struther/Desktop/Link to 4610/Hiland-EileanDonan.jpg"]
Data=Reverse[N[Pic[[1,All,All,1]]]];

ListDensityPlot[Data]

Singula

ValCount=75;
Timing[{u,w,v}=SingularValueDecomposition[Data,ValCount];]
FilteredData=u.w.Transpose[v];
ListDensityPlot[FilteredData]

ValCount=4;
Timing[{u,w,v}=SingularValueDecomposition[Data,ValCount];]
FilteredData=u.w.Transpose[v];
ListDensityPlot[FilteredData]