BandedMatrices.jl icon indicating copy to clipboard operation
BandedMatrices.jl copied to clipboard

Type of factors in Cholesky decomposition

Open KlausC opened this issue 5 years ago • 1 comments

The types of the factors of a Cholesky decomposition of a banded matrix do not reflect the banded layout, but produce a full matrix in half of the cases. Would be nice to have that fixed.

julia> A = Symmetric(BandedMatrix(0=>ones(2)))
2×2 Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅ 
  ⋅   1.0

julia> C = cholesky(A)
Cholesky{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}
U factor:
2×2 UpperTriangular{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅ 
  ⋅   1.0

julia> typeof(C.L)
LowerTriangular{Float64,Array{Float64,2}} # that is NOK

julia> typeof(C.U)
UpperTriangular{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}

julia> B = Symmetric(BandedMatrix(0=>ones(2)), :L)
2×2 Symmetric{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅ 
  ⋅   1.0

julia> D = cholesky(B)
Cholesky{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}
L factor:
2×2 LowerTriangular{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅ 
  ⋅   1.0

julia> D.L
2×2 LowerTriangular{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅ 
  ⋅   1.0

julia> D.U
2×2 UpperTriangular{Float64,Array{Float64,2}}:  # That is NOK
 1.0  0.0
  ⋅   1.0

KlausC avatar Sep 07 '20 11:09 KlausC

Note that

julia> @which getproperty(C, :U)
getproperty(C::Cholesky, d::Symbol) in LinearAlgebra at /Users/sheehanolver/Projects/julia-1.5/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/cholesky.jl:411

and following the code we see the real issue is

julia> C.factors'
2×2 Adjoint{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅ 
  ⋅   1.0

julia> copy(C.factors')
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  1.0

This is easy to fix:

julia> Base.copy(Ac::Adjoint{<:Any,<:AbstractBandedMatrix}) = BandedMatrix(Ac)

julia> copy(C.factors')
2×2 BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}:
 1.0   ⋅ 
  ⋅   1.0

julia> C.U
2×2 UpperTriangular{Float64,BandedMatrix{Float64,Array{Float64,2},Base.OneTo{Int64}}}:
 1.0   ⋅ 
  ⋅   1.0

Are you able to make a PR for this? We should also fix copy(transpose(C.factors)) for completeness.

dlfivefifty avatar Sep 07 '20 13:09 dlfivefifty