The standard implementation of ldl factorization in Matlab does not work correctly

The ldl(P) applied to symmetric matrix P should return a lover-triangular matrix L and a diagonal matrix D, such that:

.

The attached file contain the 9x9 symmetric matrix P, which ldl doesn't works with correctly:

load('ldl_fail.mat');
P % original symmetric matrix

P = 9x9

1.0e-04 * 0.0029 -0.0000 -0.0001 -0.0000 -0.0049 0.0000 0.0000 -0.0014 0.0000 -0.0000 0.0029 -0.0000 0.0051 -0.0000 0.0003 0.0014 0.0000 0.0000 -0.0001 -0.0000 0.0034 0.0000 -0.0008 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0051 0.0000 0.0577 -0.0000 -0.0019 0.0025 0.0000 0.0000 -0.0049 -0.0000 -0.0008 -0.0000 0.0573 0.0000 0.0000 0.0024 -0.0000 0.0000 0.0003 0.0000 -0.0019 0.0000 0.1050 0.0001 -0.0000 0.0001 0.0000 0.0014 -0.0000 0.0025 0.0000 0.0001 0.0098 -0.0000 0.0000 -0.0014 0.0000 0.0000 0.0000 0.0024 -0.0000 -0.0000 0.0098 -0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0001 0.0000 -0.0000 0.0091

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

max(max(P-P'))

ans = 5.2940e-23

[L,D] = ldl(P); 
L % incorrect lover-triangular matrix factor

L = 9x9

-0.0863 -0.0000 -0.0435 -0.0000 1.0000 0 0 0 0 -0.0000 0.0884 -0.0000 1.0000 0 0 0 0 0 -0.0142 0.0000 1.0000 0 0 0 0 0 0 -0.0000 1.0000 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 0 0.0000 -0.0321 0.0000 0.1834 0.0000 1.0000 0 0 0 0.0000 0.0430 -0.0001 0.4812 0.0000 0.0000 1.0000 0 0 0.0421 0.0000 0.0141 0.0000 -0.4812 -0.0000 0.0000 1.0000 0 -0.0000 0.0006 -0.0000 0.0069 0.0001 0.0005 0.0000 -0.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

diag(D)'

ans = 1x9

1.0e-04 * 0.0573 0.0577 0.0033 0.0025 0.0025 0.1049 0.0091 0.0091 0.0091

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

L_ = chol(P)'; D_ = diag(L_).^2; L_ = tril(L,-1)+eye(size(L));
L_ % correct lover-triangular matrix factor (obtained with chol() function)

L_ = 9x9

1.0000 0 0 0 0 0 0 0 0 -0.0000 1.0000 0 0 0 0 0 0 0 -0.0142 0.0000 1.0000 0 0 0 0 0 0 -0.0000 1.0000 0 1.0000 0 0 0 0 0 1.0000 0 0 0 1.0000 0 0 0 0 0.0000 -0.0321 0.0000 0.1834 0.0000 1.0000 0 0 0 0.0000 0.0430 -0.0001 0.4812 0.0000 0.0000 1.0000 0 0 0.0421 0.0000 0.0141 0.0000 -0.4812 -0.0000 0.0000 1.0000 0 -0.0000 0.0006 -0.0000 0.0069 0.0001 0.0005 0.0000 -0.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D_'

ans = 1x9

1.0e-04 * 0.0029 0.0029 0.0033 0.0488 0.0487 0.1049 0.0091 0.0091 0.0091

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Answer 1

[Torsten]

From the documentation:

[L,D] = ldl(A) factorizes full matrix A into a permuted lower triangular matrix L and a block diagonal matrix D satisfying A = L*D*L'.

with stress on "permuted".