Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Type of factors in Cholesky decomposition #202

Open
KlausC opened this issue Sep 7, 2020 · 1 comment
Open

Type of factors in Cholesky decomposition #202

KlausC opened this issue Sep 7, 2020 · 1 comment

Comments

@KlausC
Copy link
Contributor

KlausC commented Sep 7, 2020

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

@dlfivefifty
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants