Does LDLt function from LinearAlgebra package only work with SymTriDiagonal matrices?

user_1477083 0
Viewed 74

I am studying LDL^t factorization of a symmetric matrix. My own code works correctly, however, when I want to use LDLt function from LinearAlgebra package, the following code does not work

using LinearAlgebra

A = [5.0 1 0 0 0; 1 5 1 0 0; 0 1 5 1 0; 0 0 1 5 1; 0 0 0 1 5]

ldltS = LDLt(A)

display(ldltS)

Julia reports the error "ERROR: LoadError: type Array has no field dv".

If instead of A as in the above code, I build A with SymTriDiagonal, then the code above works.

julia linear-algebra
2 Answers

your code works, try:

ldltS = LDLt(A);

However the documentation for LDLt reads:

Matrix factorization type of the LDLt factorization of a real SymTridiagonal matrix S

If you want to force it to work with regular matrices the problem is that show displays the object assuming it to be a SymTridiagonal matrix. See the getproperty method in LinearAlgebra\src\ldlt.jl:

function getproperty(F::LDLt, d::Symbol)      
    Fdata = getfield(F, :data)                
    if d === :d                               
        return Fdata.dv                       
    elseif d === :D                           
        return Diagonal(Fdata.dv)             
    elseif d === :L                           
        return UnitLowerTriangular(Fdata)     
    elseif d === :Lt                          
        return UnitUpperTriangular(Fdata)     
    else                                      
        return getfield(F, d)                 
    end                                       
end

The line:

return Diagonal(Fdata.dv)

could be replaced with

return Diagonal(Fdata)

You could just overwrite this function:

julia> function Base.getproperty(F::LDLt, d::Symbol)
           Fdata = getfield(F, :data)
           if d === :d
               return Fdata.dv
           elseif d === :D
               return Diagonal(Fdata)
           elseif d === :L
               return UnitLowerTriangular(Fdata)
           elseif d === :Lt
               return UnitUpperTriangular(Fdata)
           else
               return getfield(F, d)
           end
       end

now it works (with regards that algebra in LDLt assumes that this is a SymTridiagonal matrix) :

julia> ldltS
LDLt{Float64, Matrix{Float64}}
L factor:
5×5 UnitLowerTriangular{Float64, Matrix{Float64}}:
 1.0   ⋅    ⋅    ⋅    ⋅
 1.0  1.0   ⋅    ⋅    ⋅
 0.0  1.0  1.0   ⋅    ⋅
 0.0  0.0  1.0  1.0   ⋅
 0.0  0.0  0.0  1.0  1.0
D factor:
5×5 Diagonal{Float64, Vector{Float64}}:
 5.0   ⋅    ⋅    ⋅    ⋅
  ⋅   5.0   ⋅    ⋅    ⋅
  ⋅    ⋅   5.0   ⋅    ⋅
  ⋅    ⋅    ⋅   5.0   ⋅
  ⋅    ⋅    ⋅    ⋅   5.0
user_99577100

You can use LDLFactorizations.ldl:

julia> using LDLFactorizations

julia> M = [55 225 730; 225 979 3162; 730 3162 10216.0]
3×3 Matrix{Float64}:
  55.0   225.0    730.0
 225.0   979.0   3162.0
 730.0  3162.0  10216.0

julia> LDLT = ldl(M)
LDLFactorizations.LDLFactorization{Float64, Int64, Int64, Int64}(true, true, false, 3, [2, 3, -1], [2, 1, 0], [3, 3, 3], [1, 2, 3], [1, 2, 3], [1, 3, 4, 4], Int64[], Int64[], [2, 3, 3], [4.090909090909091, 13.272727272727273, 2.9999999999999942], [55.0, 58.54545454545462, 5.684341886080801e-13], [0.0, 0.0, 0.0], [1, 1, 2], 0.0, 0.0, 0.0, 3)

julia> LDLT.P # identity, so we can ignore
3-element Vector{Int64}:
 1
 2
 3

julia> (LDLT.L+I) * LDLT.D * LinearAlgebra.transpose(LDLT.L+I)
3×3 SparseArrays.SparseMatrixCSC{Float64, Int64} with 9 stored entries:
  55.0   225.0    730.0
 225.0   979.0   3162.0
 730.0  3162.0  10216.0
user_11001070
Related