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

Improve printing of LDLt factorization #32934

Merged
merged 2 commits into from
Aug 20, 2019

Conversation

goggle
Copy link
Contributor

@goggle goggle commented Aug 16, 2019

This PR improves the REPL printing of LDLt factorizations as proposed in #24588.

using LinearAlgebra
S = SymTridiagonal([1., 2., 3.], [-1., 1.])

Prior to this PR it looked like this:

julia> F = ldlt(S)
LDLt{Float64,SymTridiagonal{Float64,Array{Float64,1}}}([1.0 -1.0 0.0; -1.0 1.0 1.0; 0.0 1.0 2.0])

With this PR it now looks like this:

LDLt{Float64,SymTridiagonal{Float64,Array{Float64,1}}}
L factor:
3×3 UnitLowerTriangular{Float64,SymTridiagonal{Float64,Array{Float64,1}}}:
  1.0        
 -1.0  1.0    
  0.0  1.0  1.0
diagonal values:
3-element Array{Float64,1}:
 1.0
 1.0
 2.0

Instead of printing diagonal values, we could call it D factor and print a diagonal matrix. Let me know which method is preferred.

I have also added a new test file stdlib/LinearAlgebra/test/ldlt.jl to add a basic test for this.

@andreasnoack
Copy link
Member

Thanks for adding this. I think it would be a good idea to have the extractors (getproperty) defined since useeers would probably expect to have a way to extract the factors when they are printed likee that. Would you be up for adding a getproperty method as well and use it in the show method?

@andreasnoack andreasnoack added domain:display and printing Aesthetics and correctness of printed representations of objects. domain:linear algebra Linear algebra labels Aug 19, 2019
@goggle
Copy link
Contributor Author

goggle commented Aug 19, 2019

@andreasnoack Yes sure, I will try to do this. If I have questions, I will ask here.

@goggle
Copy link
Contributor Author

goggle commented Aug 19, 2019

Ok, so I have added a getproperty method for LDLt.

using LinearAlgebra
S = SymTridiagonal([1., 2., 3.], [-1., 1.]);

The factorization looks like this:

julia> F = ldlt(S)
LDLt{Float64,SymTridiagonal{Float64,Array{Float64,1}}}
L factor:
3×3 UnitLowerTriangular{Float64,SymTridiagonal{Float64,Array{Float64,1}}}:
  1.0        
 -1.0  1.0    
  0.0  1.0  1.0
D factor:
3×3 Diagonal{Float64,Array{Float64,1}}:
 1.0        
     1.0    
         2.0

I'm not sure if the Lt should also be shown on that output...

I'm also not sure about the types that should be returned. Currently it is like this:

Component Type
F.L UnitLowerTriangular
F.D Diagonal
F.Lt UnitUpperTriangular
F.d Array{T, 1}

If that should be changed, let me know.

Note that something like

julia> F.L*F.D*F.Lt

does not work:

ERROR: ArgumentError: cannot set off-diagonal entry (2, 1)
Stacktrace:
 [1] copyto!(::SymTridiagonal{Float64,Array{Float64,1}}, ::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{SymTridiagonal{Float64,Array{Float64,1}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64}},typeof(*),Tuple{SymTridiagonal{Float64,Array{Float64,1}},Array{Float64,2}}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/tridiag.jl:385
 [2] materialize! at ./broadcast.jl:817 [inlined]
 [3] rmul!(::SymTridiagonal{Float64,Array{Float64,1}}, ::Diagonal{Float64,Array{Float64,1}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/diagonal.jl:175
 [4] rmul!(::UnitLowerTriangular{Float64,SymTridiagonal{Float64,Array{Float64,1}}}, ::Diagonal{Float64,Array{Float64,1}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/diagonal.jl:187
 [5] * at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.3/LinearAlgebra/src/diagonal.jl:165 [inlined]
 [6] *(::UnitLowerTriangular{Float64,SymTridiagonal{Float64,Array{Float64,1}}}, ::Diagonal{Float64,Array{Float64,1}}, ::UnitUpperTriangular{Float64,SymTridiagonal{Float64,Array{Float64,1}}}) at ./operators.jl:529
 [7] top-level scope at REPL[8]:1

I would have assumed that the multiplication gets dispatched to a more generic method if there is no specialized method... What's the reason for this? Is it some bug?

@andreasnoack
Copy link
Member

Great work. The error is probably because we don't have specialized methods for UnitLowerTriangular{Float64,SymTridiagonal{Float64,Array{Float64,1}}}. Eventually, it would be good to have them but it can go into a separate PR.

This pull request was closed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
domain:display and printing Aesthetics and correctness of printed representations of objects. domain:linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants