Cholesky factorization: How to change code to produce A = L' * D * L instead of A = L * D * L' .

6 views (last 30 days)
Alvin Ang
Alvin Ang on 18 Dec 2021
Commented: Alvin Ang on 19 Dec 2021
Hi guys i want to obtain A = L' * D * L with L: unit lower triangular.
However, the code i found was for A = L * D * L':
% LDL factorization of a symmetrical matrix A.
%
% [L, D] = LDL_GOLUB(A) returns lower triangular matrix L and diagonal
% matrix D for symmetric matrix A. It should hold that:
% L*D*L' = A.
%
% The implementation is a direct rewrite from the book.
% No attempts were made to optimize the function.
%
% Example:
% n = 4; % size of the generated matrix
% x = randi(5,n,n); % square matrix
% A = x*x'; % symmetrical matrix
% [L, D] = ldl_golub(A);
% obtained = L*D*L';
% assert(all(all(abs(A - obtained) < eps(200))))
% From Golub 1996:
% Matrix Computations,
% algorithm 4.1.2 (in revision from 2013, it is algorithm 4.1.1)
function [L, D] = ldl_golub(A)
n = length(A);
v = zeros(n,1);
for j = 1:n
for i = 1:j - 1
v(i) = A(j, i)*A(i, i);
end
A(j, j) = A(j, j) - A(j, 1:j - 1) * v(1:j - 1);
A(j + 1:n, j) = (A(j + 1:n,j) - A(j + 1:n, 1:j - 1) * v(1:j - 1))/A(j, j);
end
L = tril(A,-1)+eye(n);
D = diag(diag(A));
  2 Comments
Bruno Luong
Bruno Luong on 19 Dec 2021
Hint: if you want to achieve that you need to maje the reverse loop.
for j = n-1:-1:1
for i = j+1:n
...
end
end

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 18 Dec 2021
% Generate random symetric A
A=randn(5); A=A'+A;
disp(A)
-3.1881 1.0281 0.2981 -1.6865 -0.6830 1.0281 2.5594 0.6286 -0.0361 2.1537 0.2981 0.6286 -0.2223 0.7094 -0.7746 -1.6865 -0.0361 0.7094 0.5899 2.0677 -0.6830 2.1537 -0.7746 2.0677 -3.1272
[L,D]=ldl(A);
disp(L*D*L')
-3.1881 1.0281 0.2981 -1.6865 -0.6830 1.0281 2.5594 0.6286 -0.0361 2.1537 0.2981 0.6286 -0.2223 0.7094 -0.7746 -1.6865 -0.0361 0.7094 0.5899 2.0677 -0.6830 2.1537 -0.7746 2.0677 -3.1272
[L,D]=ldl(rot90(A,2)); L=rot90(L,2)'; D=rot90(D,2);
disp(L'*D*L)
-3.1881 1.0281 0.2981 -1.6865 -0.6830 1.0281 2.5594 0.6286 -0.0361 2.1537 0.2981 0.6286 -0.2223 0.7094 -0.7746 -1.6865 -0.0361 0.7094 0.5899 2.0677 -0.6830 2.1537 -0.7746 2.0677 -3.1272
  4 Comments
Alvin Ang
Alvin Ang on 19 Dec 2021
Understood, thank you! I have achieved A = L' * D * L using a self written function together with your advice:
[L,D]=ldlt(rot90(A,2));
L=rot90(L,2)'; D=rot90(D,2);
disp(L'*D*L)
function [L,D]=ldlt(A)
%
% Figure out the size of A.
%
n=size(A,1);
%
% The main loop. See Golub and Van Loan for details.
%
L=zeros(n,n);
for j=1:n,
if (j > 1),
v(1:j-1)=L(j,1:j-1).*d(1:j-1);
v(j)=A(j,j)-L(j,1:j-1)*v(1:j-1)';
d(j)=v(j);
if (j < n),
L(j+1:n,j)=(A(j+1:n,j)-L(j+1:n,1:j-1)*v(1:j-1)')/v(j);
end;
else
v(1)=A(1,1);
d(1)=v(1);
L(2:n,1)=A(2:n,1)/v(1);
end;
end;
%
% Put d into a matrix.
%
D=diag(d);
%
% Put ones on the diagonal of L.
%
L=L+eye(n);
end

Sign in to comment.

More Answers (1)

John D'Errico
John D'Errico on 18 Dec 2021
Edited: John D'Errico on 18 Dec 2021
The trivial solution is...
You want to compute a U*D*U' factorization, where U is UPPER triangular. MATLAB already offers LDL, which returns a LOWER triangular L.
As a simple example, I'll create the test matrix A.
A = rand(5) + eye(5); A = A + A'
A = 5×5
2.1794 1.2047 1.1618 0.8876 1.3451 1.2047 2.1395 0.8523 1.4961 1.0057 1.1618 0.8523 3.7328 0.3476 0.9420 0.8876 1.4961 0.3476 2.4837 1.2139 1.3451 1.0057 0.9420 1.2139 3.4335
We can use LDL. Apply it to a modified version of A.
Q = flip(eye(size(A)));
[L,D] = ldl(Q*A*Q);
Next, form these matrices as products, using Q, L, and D.
U = Q*L*Q;
Dhat = Q*D*Q;
Is U upper triangular? Of course. As long as L was lower triangular, then U is upper triangular, as you want.
U
U = 5×5
1.0000 0.4047 0.2273 0.2005 0.3918 0 1.0000 0.1636 0.5551 0.2929 0 0 1.0000 0.0071 0.2744 0 0 0 1.0000 0.3535 0 0 0 0 1.0000
Do these matices have the desired property? Yes.
U*Dhat*U'
ans = 5×5
2.1794 1.2047 1.1618 0.8876 1.3451 1.2047 2.1395 0.8523 1.4961 1.0057 1.1618 0.8523 3.7328 0.3476 0.9420 0.8876 1.4961 0.3476 2.4837 1.2139 1.3451 1.0057 0.9420 1.2139 3.4335
norm(U*Dhat*U' - A)
ans = 5.2687e-16
So zero to within floating point trash.
No code mods were required. If you have LDL, then you have a simple way to compute a UDU factorization. All of this works because the matrix Q=Q' is idempotent, so Q*Q equals the identity matrix. That means if we have
A = L*D*L'
then just pre and post multiply by Q.
Q*A*Q = Q*L*D*L'*Q
Next, insert extra pairs of Q*Q between the inner terms. Remember that Q is idempotent, so that changes nothing.
Q*A*Q = Q*L*Q*Q*D*Q*Q*L'*Q
Next, group some terms as...
Q*A*Q = (Q*L*Q)*(Q*D*Q)*(Q*L'*Q)
Finally, recognize that if L is lower triangular, then Q*L*Q is UPPER triangular.
Oh. Is this homework, and you really need to write code, with loops and all that crap? I hope not, but then this may not help you to do your homework assignment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by