Issue with low-rank factorizations

[julian-upc]

The positive_semidefinite example in test_real.jl appears to fail the invariant $$P A P^T \ = \ (L+I)D(L+I)^T.$$ The example is

A = [
      0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 4.0 0.0
      0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 5.0 0.0
      2.0 4.0 5.0 -2 4.0 1.0 2.0 2.0 2.0 0.0
      0.0 0.0 0.0 0.0 1.0 9.0 9.0 1.0 7.0 1.0
      0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
      1.0 3.0 2.0 1.0 4.0 3.0 1.0 0.0 0.0 7.0
      -3.0 8.0 0.0 0.0 0.0 0.0 -2.0 0.0 0.0 1.0
      0.0 0.0 0.0 5.0 7.0 9.0 0.0 2.0 7.0 1.0
      3.0 2.0 0.0 0.0 0.0 0.0 1.0 3.0 3.0 2.0
      0.0 0.0 0.0 0.0 -3 -4 0.0 0.0 0.0 0.0
    ]
ms = SparseMatrixCSC{Float64, Int64}(A*A')
10×10 SparseMatrixCSC{Float64, Int64} with 78 stored entries:
 16.0  20.0    8.0   28.0  4.0     ⋅      ⋅    28.0  12.0     ⋅ 
 20.0  25.0   10.0   35.0  5.0     ⋅      ⋅    35.0  15.0     ⋅ 
  8.0  10.0   78.0   47.0  2.0   43.0   22.0   45.0  28.0  -16.0
 28.0  35.0   47.0  214.0  7.0   47.0  -17.0  140.0  35.0  -39.0
  4.0   5.0    2.0    7.0  1.0     ⋅      ⋅     7.0   3.0     ⋅ 
   ⋅     ⋅    43.0   47.0   ⋅    90.0   26.0   67.0  24.0  -24.0
   ⋅     ⋅    22.0  -17.0   ⋅    26.0   78.0    1.0   7.0     ⋅ 
 28.0  35.0   45.0  140.0  7.0   67.0    1.0  209.0  29.0  -57.0
 12.0  15.0   28.0   35.0  3.0   24.0    7.0   29.0  36.0     ⋅ 
   ⋅     ⋅   -16.0  -39.0   ⋅   -24.0     ⋅   -57.0    ⋅    25.0

Now let's do the LDLT decomposition and check:

LDLT = ldlt(ms)

The left-hand side of the identity is

permute(ms, LDLT.P, LDLT.P)
10×10 SparseMatrixCSC{Float64, Int64} with 78 stored entries:
  25.0     ⋅   -24.0    ⋅     ⋅   -16.0  -39.0   ⋅   -57.0    ⋅ 
    ⋅    78.0   26.0    ⋅     ⋅    22.0  -17.0   ⋅     1.0   7.0
 -24.0   26.0   90.0    ⋅     ⋅    43.0   47.0   ⋅    67.0  24.0
    ⋅      ⋅      ⋅   16.0  20.0    8.0   28.0  4.0   28.0  12.0
    ⋅      ⋅      ⋅   20.0  25.0   10.0   35.0  5.0   35.0  15.0
 -16.0   22.0   43.0   8.0  10.0   78.0   47.0  2.0   45.0  28.0
 -39.0  -17.0   47.0  28.0  35.0   47.0  214.0  7.0  140.0  35.0
    ⋅      ⋅      ⋅    4.0   5.0    2.0    7.0  1.0    7.0   3.0
 -57.0    1.0   67.0  28.0  35.0   45.0  140.0  7.0  209.0  29.0
    ⋅     7.0   24.0  12.0  15.0   28.0   35.0  3.0   29.0  36.0

but the right-hand side is

(LDLT.L+I) * LDLT.D * (LDLT.L+I)'
10×10 SparseMatrixCSC{Float64, Int64} with 16 stored entries:
  25.0    ⋅   -24.0    ⋅     ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
    ⋅   78.0   26.0    ⋅     ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
 -24.0  26.0   90.0    ⋅     ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
    ⋅     ⋅      ⋅   16.0  20.0   ⋅    ⋅    ⋅    ⋅    ⋅ 
    ⋅     ⋅      ⋅   20.0  25.0   ⋅    ⋅    ⋅    ⋅    ⋅ 
    ⋅     ⋅      ⋅     ⋅     ⋅   0.0   ⋅    ⋅    ⋅    ⋅ 
    ⋅     ⋅      ⋅     ⋅     ⋅    ⋅   0.0   ⋅    ⋅    ⋅ 
    ⋅     ⋅      ⋅     ⋅     ⋅    ⋅    ⋅   0.0   ⋅    ⋅ 
    ⋅     ⋅      ⋅     ⋅     ⋅    ⋅    ⋅    ⋅   0.0   ⋅ 
    ⋅     ⋅      ⋅     ⋅     ⋅    ⋅    ⋅    ⋅    ⋅   0.0

We can see that they agree on the first 5 x 5 block, but disagree later.

This is probably due to the fact that D has too many zeros on the diagonal:

LDLT.D
LDLT.D
10×10 Diagonal{Float64, Vector{Float64}}:
 25.0    ⋅     ⋅        ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅   78.0    ⋅        ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅     ⋅   58.2933    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅     ⋅     ⋅      16.0   ⋅    ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅     ⋅     ⋅        ⋅   0.0   ⋅    ⋅    ⋅    ⋅    ⋅ 
   ⋅     ⋅     ⋅        ⋅    ⋅   0.0   ⋅    ⋅    ⋅    ⋅ 
   ⋅     ⋅     ⋅        ⋅    ⋅    ⋅   0.0   ⋅    ⋅    ⋅ 
   ⋅     ⋅     ⋅        ⋅    ⋅    ⋅    ⋅   0.0   ⋅    ⋅ 
   ⋅     ⋅     ⋅        ⋅    ⋅    ⋅    ⋅    ⋅   0.0   ⋅ 
   ⋅     ⋅     ⋅        ⋅    ⋅    ⋅    ⋅    ⋅    ⋅   0.0

[geoffroyleconte]

Hi! I assume you meant LDLT = ldl(ms)? ldlt is the LinearAlgebra function. I think this is to be expected because A is rank-deficient, so A*A' is too. That is why in the tests we used dynamic regularization: https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl/blob/7a8ee84565ba3e021bdf1997af19171739777280/test/test_real.jl#L223 Here is another example using the dynamic regularization: https://juliasmoothoptimizers.github.io/LDLFactorizations.jl/stable/tutorial/#Dynamic-Regularization

In the "ldl_mul!" testset we chose to solve systems with the right-hand side b in the image of the matrix to factorize: https://github.com/JuliaSmoothOptimizers/LDLFactorizations.jl/blob/7a8ee84565ba3e021bdf1997af19171739777280/test/test_real.jl#L340

Which example and testset were you referring to?

[julian-upc]

Hi Geoffroy,

Thank you and I appreciate the quick response!

I see, and sorry for the confusion between ldl and ldlt. Maybe this is not actually a bug, but just how ldl works.

However, in my application, solving for right-hand sides is not enough, and I really need the factorization property. (Briefly, I want to decompose $A = L\sqrt{D}(L\sqrt{D})^ =: CC^$ so that I can block-unitarize a possibly degenerate scalar product $(x,Ay) = ( C^x, C^y)$.)

Consider the smaller example

using LinearAlgebra, SparseArrays, LDLFactorizations

A = SparseMatrixCSC{Float64, Int64}([0 0 0 0; 0 5.0 -2.88675 -5.16389; 0 -2.88675 1.66667 2.98142; 0 -5.16389 2.98142 5.33333])
4×4 SparseMatrixCSC{Float64, Int64} with 9 stored entries:
  ⋅     ⋅         ⋅         ⋅ 
  ⋅    5.0      -2.88675  -5.16389
  ⋅   -2.88675   1.66667   2.98142
  ⋅   -5.16389   2.98142   5.33333

Since $A$ only has rank 3 and the first row and column are zero, ldlt fails:

ldlt(A)
ERROR: ZeroPivotException: factorization encountered one or more zero pivots. Consider switching to a pivoted LU factorization.

On the other hand, ldl succeeds, but gives back a zero matrix:

x=ldl(A);
x.L
4×4 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
  ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅ 
  ⋅    ⋅    ⋅    ⋅ 
x.D
4×4 Diagonal{Float64, Vector{Float64}}:
 0.0   ⋅         ⋅             ⋅ 
  ⋅   4.4e-323   ⋅             ⋅ 
  ⋅    ⋅        6.92037e-310   ⋅ 
  ⋅    ⋅         ⋅            6.92037e-310

This doesn't change even if I use ldl_analyze / ldl_factorize! as per the tutorial you linked to:

e = sqrt(eps(eltype(A)))
x.r1 = -e
x.r2 = e
x.tol = e
S=ldl_analyze(A)
S=ldl_factorize!(A, S)
S.L
4×4 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
  ⋅   6.92037e-310  6.92037e-310   ⋅ 
  ⋅    ⋅             ⋅             ⋅ 
  ⋅    ⋅             ⋅             ⋅ 
  ⋅    ⋅             ⋅             ⋅ 
S.D
4×4 Diagonal{Float64, Vector{Float64}}:
 0.0   ⋅             ⋅             ⋅ 
  ⋅   6.92037e-310   ⋅             ⋅ 
  ⋅    ⋅            6.92037e-310   ⋅ 
  ⋅    ⋅             ⋅            6.92037e-310

However, after permuting, we do get a factorization:

x=ldl(permute(A,[4,3,2,1],[4,3,2,1]));
x.L
4×4 SparseMatrixCSC{Float64, Int64} with 3 stored entries:
   ⋅          ⋅        ⋅    ⋅ 
  0.559017    ⋅        ⋅    ⋅ 
 -0.96823   -7.39535   ⋅    ⋅ 
   ⋅          ⋅        ⋅    ⋅ 
x.D
4×4 Diagonal{Float64, Vector{Float64}}:
 5.33333   ⋅            ⋅            ⋅ 
  ⋅       6.73026e-6    ⋅            ⋅ 
  ⋅        ⋅          -0.000201198   ⋅ 
  ⋅        ⋅            ⋅           0.0

However, because of the structure of L, the reconstructed matrix $LDL^T$ only has rank 2.

So right now I have two questions:

1. (related to your code) Is there a way to get a rank 3 LDL factorization of $A$? Maybe I used the dynamic regularization in a wrong way.
2. (related to my problem, on the off chance that you might know) Is there a pivot-aware Cholesky factorization available in julia?

Thanks a lot for your time, I appreciate it!

[julian-upc]

Sorry, I make an invalid comment here.

[geoffroyleconte]

Your matrix A is (almost) of rank 1 I think. Row 3 and 4 are multiple of row 2:

julia> [5.0, -2.88675, -5.16389] ./ [-5.16389, 2.98142, 5.33333]
3-element Vector{Float64}:
 -0.9682622983835829
 -0.9682466744034722
 -0.9682299801437376

julia> [5.0, -2.88675, -5.16389] ./ [-2.88675,   1.66667,   2.98142]
3-element Vector{Float64}:
 -1.732051615138131
 -1.7320465359069281
 -1.732023666574988

Also the LDL factorization is defined for symmetric positive/negative definite matrices, or more generally symmetric quasi-definite matrices. If your matrix to factorize is not symmetric quasi-definite, it might not have a LDL decomposition.

[dpo]

@julian-upc Are you looking to preserve semi-definiteness in the factors?

[julian-upc]

Oh I see. I wasn't aware that the factorization might not exist (not an expert on numerical linear algebra). I think for my application it is enough to consider the full-rank case.

Thank you all again for your time! Best, Julian