Add methods for dividing triangular by diagonal matrices. #27999

tpapp · 2018-07-09T07:48:51Z

fredrikekre · 2018-07-09T08:46:26Z

stdlib/LinearAlgebra/test/diagonal.jl

@@ -455,4 +455,21 @@ end
 @test Transpose(x)*D*x == (Transpose(x)*D)*x == (Transpose(x)*Array(D))*x
 end

+@testset "Triangular division by Diagonal #27989" begin
+ for _ in 1:100


This loop seems unnecessary?

You mean that testing with multiple random values is unnecessary, and I should go for a single one?

The major cost is compilation, once compiled, additional random tests are relatively cheap.

Right, but we don't do this anywhere else. For testing correctness of the code path it would be more useful with a loop over e.g. (Float64, ComplexF64).

Thanks for the suggestion, this actually caught a bug.

fredrikekre · 2018-07-09T08:47:40Z

stdlib/LinearAlgebra/test/diagonal.jl

+ L = LowerTriangular(randn(K, K))
+ D = Diagonal(randn(K))
+ @test U / D isa UpperTriangular
+ @test U / D == UpperTriangular(Matrix(U) / Matrix(D))


Perhaps combine these two (and similar below) as

@test (U / D)::UpperTriangular == UpperTriangular(Matrix(U) / Matrix(D))

for the same effect?

@fredrikekre

suggested by @fredrikekre

For broadcasting, tranpose should be used instead of adjoint.

andreasnoack · 2018-07-10T09:54:29Z

stdlib/LinearAlgebra/src/diagonal.jl

@@ -356,6 +356,11 @@ rdiv!(A::AbstractMatrix{T}, transD::Transpose{<:Any,<:Diagonal{T}}) where {T} =
 (\)(A::Union{QR,QRCompactWY,QRPivoted}, B::Diagonal) =
 invoke(\, Tuple{Union{QR,QRCompactWY,QRPivoted}, AbstractVecOrMat}, A, B)

+(/)(U::UpperTriangular, D::Diagonal) = UpperTriangular(parent(U) ./ transpose(D.diag))


It would be great if you could define the in-place versions ldiv! and rdiv! instead and then adjust the methods in line 427-429 to use these.

Probably just something like

function rdiv!(U::UpperTriangular, D::Diagonal) broadcast!(/, parent(U), parent(U), diag(D)) return U end

although using broadcast here is actually cheating a bit since broadcasting and linear algebra operations have slightly different semantics in corner cases. However, that is an existing issue so we can clean that up in another PR.

Hm. The transpose is actually also not correct here since it is recursive so you probably want to reshape the vector instead.

Thanks for the review. I did the change you requested, but I am not convinced that going though mutating methods like rdiv! and ldiv! is the right approach.

Sure, it removes duplicate code, but at the same time the required promotion code makes the implementation convoluted (cf the rather simple approach with the commit I just made). Maybe I am doing the promotion wrong?

We'd need the rdiv! and ldiv! methods anyway and it would be confusing if the non-mutating versions didn't go through the mutating versions. Also, the promotion can be done in the more generic version so it's not that much code. Good idea to use broadcast of the type for that.

Removed previous direct implementations for upper/lower triangular, added ldiv!/rdiv! methods. Use reshape instead of transpose.

tpapp · 2018-07-11T08:54:01Z

I did not run the full test suite before committing, and now the (freebsd) tests show an ambiguity which I expect will crop up in the other CI runs, too.

Suggestions on how to handle this would be appreciated.

andreasnoack · 2018-07-11T09:22:49Z

You'd need to add definitions for

(/)(A::Transpose, D::Diagonal)
(/)(A::Adjoint, D::Diagonal)
(/)(x::RowVector, D::Diagonal)

where the latter should be defined as a deprecated method. The two former might be joined to AdjOrTrans but I'm not sure.

tpapp · 2018-07-11T09:48:30Z

Thanks, I am aware that defining the intersection is the standard solution for ambiguities, but I was wondering if there is a better method for streamlining the architecture.

I am wondering about rolling back to 985891d, replacing the transpose with reshape, and leaving the rdiv! and ldiv! for future work.

I am also skeptical about whether we need an ldiv! and rdiv! method at all, since eg ldiv! says that

The argument A should not be a matrix

and UpperTriangular and LowerTriangular are matrices.

andreasnoack · 2018-07-11T10:04:49Z

I'm also not too happy about all the need for disambiguating methods. Maybe just restrict (/)(A::AbstractMatrix, D::Diagonal) to (/)(A::StridedMatrix, D::Diagonal). That will make the AbstractMatrix case go through the fallback that calls \ but I think that is fine.

Regarding the formulation in the docstring for ldiv! then it is not precise but it doesn't change the fact that that \ and / should go through ldiv! and rdiv!. The reason for the formulation is that it should be possible to solve directly with the first argument in ldiv!. That typically means that a matrix should be in factored form but that is, of course, not necessary for a Diagonal.

andreasnoack · 2018-07-11T10:05:56Z

stdlib/LinearAlgebra/src/diagonal.jl

@@ -325,6 +325,11 @@ ldiv!(adjD::Adjoint{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
 ldiv!(transD::Transpose{<:Any,<:Diagonal{T}}, B::AbstractVecOrMat{T}) where {T} =
 (D = transD.parent; ldiv!(D, B))

+function ldiv!(D::Diagonal, A::Union{LowerTriangular,UpperTriangular})
+ broadcast!(/, parent(A), parent(A), diag(D))


I forgot the reshape here. Needs for be reshape(D.diag, 1, length(D.diag)) instead. (I recently learned that diag now copies.)

Happy to add this, but can you please explain why reshape is needed here? In any case it should be

reshape(D.diag, length(D.diag), 1)

for the ldiv!, but wouldn't using just the vector D.diag be equivalent?

I saw you were using / and thought this was rdiv!. This should then be

broadcast!(\, parent(A), D.diag, parent(A))

d\a and a/d is not always the same. E.g. if the elements are matrices.

andreasnoack · 2018-07-11T11:26:16Z

stdlib/LinearAlgebra/src/diagonal.jl

@@ -349,7 +349,7 @@ function rdiv!(A::AbstractMatrix{T}, D::Diagonal{T}) where {T}
 end

 function rdiv!(A::Union{LowerTriangular,UpperTriangular}, D::Diagonal)
- broadcast!(/, parent(A), parent(A), reshape(diag(D), 1, :))
+ broadcast!(/, parent(A), parent(A), permutedims(D.diag))


I didn't know we now have that method.

tpapp · 2018-07-12T07:43:28Z

Thanks for all the help!

Add methods for dividing triangular by diagonal matrices.

e6ececc

Fixes JuliaLang#27989.

fredrikekre reviewed Jul 9, 2018

View reviewed changes

tpapp added 2 commits July 9, 2018 11:14

Simplified tests.

f1ebb96

suggested by @fredrikekre

Test for element types instead, fix bug.

985891d

For broadcasting, tranpose should be used instead of adjoint.

fredrikekre approved these changes Jul 10, 2018

View reviewed changes

andreasnoack reviewed Jul 10, 2018

View reviewed changes

kshyatt added the domain:linear algebra Linear algebra label Jul 10, 2018

Matrix diagonal division through ldiv!/rdiv!.

452fb16

Removed previous direct implementations for upper/lower triangular, added ldiv!/rdiv! methods. Use reshape instead of transpose.

andreasnoack reviewed Jul 11, 2018

View reviewed changes

Fix ldiv! to \, use permutedims, restrict signature.

4ebd99d

andreasnoack reviewed Jul 11, 2018

View reviewed changes

andreasnoack merged commit 5f19151 into JuliaLang:master Jul 12, 2018

tpapp deleted the tp/fix-diagonal-div branch July 12, 2018 07:43

hbarthels mentioned this pull request Aug 28, 2019

Dividing diagonal by triangular matrix returns sparse matrix #33101

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add methods for dividing triangular by diagonal matrices. #27999

Add methods for dividing triangular by diagonal matrices. #27999

tpapp commented Jul 9, 2018

fredrikekre Jul 9, 2018

tpapp Jul 9, 2018

fredrikekre Jul 9, 2018

tpapp Jul 9, 2018

fredrikekre Jul 9, 2018

andreasnoack Jul 10, 2018

andreasnoack Jul 10, 2018 •

edited

Loading

andreasnoack Jul 10, 2018

tpapp Jul 11, 2018

andreasnoack Jul 11, 2018

tpapp commented Jul 11, 2018

andreasnoack commented Jul 11, 2018

tpapp commented Jul 11, 2018

andreasnoack commented Jul 11, 2018

andreasnoack Jul 11, 2018

tpapp Jul 11, 2018

andreasnoack Jul 11, 2018

andreasnoack Jul 11, 2018

tpapp commented Jul 12, 2018

Add methods for dividing triangular by diagonal matrices. #27999

Add methods for dividing triangular by diagonal matrices. #27999

Conversation

tpapp commented Jul 9, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andreasnoack Jul 10, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tpapp commented Jul 11, 2018

andreasnoack commented Jul 11, 2018

tpapp commented Jul 11, 2018

andreasnoack commented Jul 11, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tpapp commented Jul 12, 2018

andreasnoack Jul 10, 2018 •

edited

Loading