BandedMatrices.jl
BandedMatrices.jl copied to clipboard
Type of factors in Cholesky decomposition
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
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.