Cholesky factorization: How to change code to produce A = L' * D * L instead of A = L * D * L' .
6 views (last 30 days)
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.
% 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);
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);
L = tril(A,-1)+eye(n);
D = diag(diag(A));
More Answers (1)
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'
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.
Do these matices have the desired property? Yes.
norm(U*Dhat*U' - A)
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.