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

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

[Alvin Ang]

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]

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

Answer 1

[Bruno Luong]

% 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

[John D'Errico]

Thanks. This is simpler than mine, but has the same effect, since the rot90(X,2) does the same thing as my pair of matrix multiplies. I had a feeling there was an alternate solution like that to be found.

[Alvin Ang]

Hi Bruno & John, thanks for the replies. However my ldl function can't seem to work on my A matrix.

Error using ldl

Matrix must have real diagonal.

Here is my code that displays that error:

%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(A)

[Bruno Luong]

"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]

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

Answer 2

[John D'Errico]

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.

https://en.wikipedia.org/wiki/Idempotence

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.

[Alvin Ang]

Thank you so much for your advice!.