no method for A*A' when A = LowerTriangular(Diagonal(...))

[tpapp]

julia> VERSION
v"0.7.1-pre.0"

julia> using LinearAlgebra

julia> A = LowerTriangular(Diagonal([1, 2]))
2×2 LowerTriangular{Int64,Diagonal{Int64,Array{Int64,1}}}:
 1  ⋅
 0  2

julia> A * A'
ERROR: MethodError: no method matching rmul!(::SparseArrays.SparseMatrixCSC{Int64,Int64}, ::Adjoint{Int64,LowerTriangular{Int64,Diagonal{Int64,Array{Int64,1}}}})

[tpapp]

BTW, I am afraid that having matrices like above comes from my own JuliaLang/julia#27823, eg

julia> using LinearAlgebra; cholesky(inv(Diagonal([1.0, 4.0]))).L
2×2 LowerTriangular{Float64,Diagonal{Float64,Array{Float64,1}}}:
 1.0   ⋅ 
 0.0  0.5

[tpapp]

I have been thinking about this issue, and broadly I can think of two solutions:

1. relax the requirement that getpropery(::Cholesky, ...) should always return an ::AbstractTriangular,

2. make Diagonal <: AbstractTriangular, as it technically is.

Also, can someone please add the linear algebra label to the issue?

[tpapp]

Friendly ping: I would be happy to make a PR fixing this, just need some guidance on the preferred solution. @andreasnoack, can you please help with this?

[andreasnoack]

I think it would be unfortunate to lose the diagonal structure which could easily happen if the Diagonal is hidden behind LowerTriangular. Hence, I think we should change getproperty to dispatch on the wrapped array type and return Diagonal when the wrapped type is Diagonal.

[dkarrasch]

@tpapp Do you still think we need this method, or was the issue kind of resolved with the diagonal cholesky decomposition?

[tpapp]

The original example still fails for me on master.

[dkarrasch]

The original example still fails for me on master.

Sure. But this is a very peculiar combination. Would you expect that one to occur in the wild?

[laborg]

@tpapp do you think this can be closed now?

julia> versioninfo()
Julia Version 1.8.0-DEV.1456
Commit 10ded70514* (2022-02-05 07:49 UTC)
...

julia> using LinearAlgebra

julia> A = LowerTriangular(Diagonal([1, 2]))
2×2 LowerTriangular{Int64, Diagonal{Int64, Vector{Int64}}}:
 1  ⋅
 0  2

julia> A * A'
2×2 Matrix{Int64}:
 1  0
 0  4

[tpapp]

@laborg: sorry for not responding earlier. Yes, I think this has been fixed somewhat accidentally by JuliaLang/julia#43127 (git bisect confirms this). I will make some a PR with a test so that this remains so.