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

20 vues (au cours des 30 derniers jours)
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 commentaires
Alvin Ang
Alvin Ang le 18 Déc 2021
Cholesky factorization
%Matrix A here is needed for factorization (A = L' * D * L)
%clear console; Remove all variables from memory; Close all figures
clc;clearvars;close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%All_Possible_Routes_On_Tree%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
TxRx = 4; %Number(Equal) of transmitters and Recievers %Number_Of_Levels_On_Tree
iter = 1; %Number of iterations/data signals
PTX = 10;
%Complex Gaussian Noise Signal, Noise
Cn = [4 1 0 2; 1 6 1 1; 0 1 5 2; 2 1 2 10];
[V,D] = eig(Cn);
squarerootCn = V*sqrt(D)*V';
Noise = 1/sqrt(2) * squarerootCn * (randn(TxRx,iter) + 1i*randn(TxRx,iter));
%Normal Distribution Matrix, H
mean = 0; variance = 1;
a = randn(TxRx) + 1i*randn(TxRx);
H = mean + sqrt(variance)*a;
%With PTX in block diagram
Hptx = sqrt(PTX/TxRx) * H;
%Random Signal (WO influence of Noise)
Integers = [0 1 2 3];
Constellations = [exp(1i*(pi/4)) exp(1i*(3.*pi/4)) exp(1i*(5.*pi/4)) exp(1i*(7.*pi/4))];
RandomIntegers = randi([0 3],iter,TxRx);
RandomConstellations = interp1(Integers, Constellations, RandomIntegers);
Sq = transpose(RandomConstellations);
%A = H' * Cn^-1 * H
A = H' * Cn^-1 * H
%A is non-negative definite if eigenvalues of A are non-negative.
eigenvalues_of_A = eig(A); %Checked: A is non-negative definite
%Cholesky factorization A = L' * D * L
[L, D] = ldl_golub(A); % but achieved A = L * D * L'
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));
end
Bruno Luong
Bruno Luong le 19 Déc 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

Connectez-vous pour commenter.

Réponse acceptée

Bruno Luong
Bruno Luong le 18 Déc 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 commentaires
Bruno Luong
Bruno Luong le 19 Déc 2021
Modifié(e) : Bruno Luong le 19 Déc 2021
"Matrix must have real diagonal"
What is not clear for you in this message? The LDL' factorization requires see a Hermitian matrix see the first sentence of doc, the Hermetian matrix has real values on diagonal, your matrix has complex value.
And this error has nothing to do with the question you asked.
Alvin Ang
Alvin Ang le 19 Déc 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

Connectez-vous pour commenter.

Plus de réponses (1)

John D'Errico
John D'Errico le 18 Déc 2021
Modifié(e) : John D'Errico le 18 Déc 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