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.

Sign up for GitHub

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

Jump to bottom

Fix multiplication of Triangular and Diagonal matrix #33334

Merged
merged 5 commits into from
Oct 10, 2019
Merged

Fix multiplication of Triangular and Diagonal matrix #33334

merged 5 commits into from
Oct 10, 2019

Conversation

goggle
Copy link
Contributor

@goggle goggle commented Sep 20, 2019

This PR intends to fix some examples from #30094 (comment), more precisely the multiplication of a triangular matrix with a diagonal matrix and vice versa.

There are more broken examples involving triangular and diagonal matrices, but I first wanted to check if this is done right.

Here are some examples that still need to be addressed:

julia> using LinearAlgebra

julia> D = Diagonal([0.5, 0.5, 0.5]);

julia> T = LowerTriangular(ones(Int, (3,3)));

julia> T'*D
ERROR: InexactError: Int64(0.5)
Stacktrace:
 [1] Int64 at ./float.jl:709 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:798 [inlined]
 [4] setindex! at ./multidimensional.jl:490 [inlined]
 [5] macro expansion at ./broadcast.jl:909 [inlined]
 [6] macro expansion at ./simdloop.jl:77 [inlined]
 [7] copyto! at ./broadcast.jl:908 [inlined]
 [8] copyto! at ./broadcast.jl:863 [inlined]
 [9] materialize! at ./broadcast.jl:822 [inlined]
 [10] rmul!(::Array{Int64,2}, ::Diagonal{Float64,Array{Float64,1}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:177
 [11] rmul! at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:187 [inlined]
 [12] *(::Adjoint{Int64,LowerTriangular{Int64,Array{Int64,2}}}, ::Diagonal{Float64,Array{Float64,1}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:219
 [13] top-level scope at REPL[4]:1

julia> transpose(T)*D
ERROR: InexactError: Int64(0.5)
Stacktrace:
 [1] Int64 at ./float.jl:709 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:798 [inlined]
 [4] setindex! at ./multidimensional.jl:490 [inlined]
 [5] macro expansion at ./broadcast.jl:909 [inlined]
 [6] macro expansion at ./simdloop.jl:77 [inlined]
 [7] copyto! at ./broadcast.jl:908 [inlined]
 [8] copyto! at ./broadcast.jl:863 [inlined]
 [9] materialize! at ./broadcast.jl:822 [inlined]
 [10] rmul!(::Array{Int64,2}, ::Diagonal{Float64,Array{Float64,1}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:177
 [11] rmul! at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:187 [inlined]
 [12] *(::Transpose{Int64,LowerTriangular{Int64,Array{Int64,2}}}, ::Diagonal{Float64,Array{Float64,1}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:228
 [13] top-level scope at REPL[5]:1

julia> D*T'
ERROR: InexactError: Int64(0.5)
Stacktrace:
 [1] Int64 at ./float.jl:709 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:798 [inlined]
 [4] setindex! at ./multidimensional.jl:490 [inlined]
 [5] macro expansion at ./broadcast.jl:909 [inlined]
 [6] macro expansion at ./simdloop.jl:77 [inlined]
 [7] copyto! at ./broadcast.jl:908 [inlined]
 [8] copyto! at ./broadcast.jl:863 [inlined]
 [9] materialize! at ./broadcast.jl:822 [inlined]
 [10] lmul!(::Diagonal{Float64,Array{Float64,1}}, ::Array{Int64,2}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:183
 [11] *(::Diagonal{Float64,Array{Float64,1}}, ::Adjoint{Int64,LowerTriangular{Int64,Array{Int64,2}}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:237
 [12] top-level scope at REPL[6]:1

julia> D*transpose(T)
ERROR: InexactError: Int64(0.5)
Stacktrace:
 [1] Int64 at ./float.jl:709 [inlined]
 [2] convert at ./number.jl:7 [inlined]
 [3] setindex! at ./array.jl:798 [inlined]
 [4] setindex!(::UpperTriangular{Int64,Array{Int64,2}}, ::Float64, ::Int64, ::Int64) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/triangular.jl:229
 [5] macro expansion at ./abstractarray.jl:1104 [inlined]
 [6] macro expansion at ./simdloop.jl:77 [inlined]
 [7] copyto! at ./broadcast.jl:908 [inlined]
 [8] copyto!(::UpperTriangular{Int64,Array{Int64,2}}, ::Base.Broadcast.Broadcasted{LinearAlgebra.StructuredMatrixStyle{UpperTriangular{Int64,Array{Int64,2}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64}},typeof(*),Tuple{Array{Float64,1},UpperTriangular{Int64,Array{Int64,2}}}}) at /home/alex/Projects/github/goggle/julia/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/structuredbroadcast.jl:175
 [9] *(::Diagonal{Float64,Array{Float64,1}}, ::Transpose{Int64,LowerTriangular{Int64,Array{Int64,2}}}) at ./broadcast.jl:822
 [10] top-level scope at REPL[8]:1

And in this example, I think we should expect in both cases a UpperTriangular as a result, right?

julia> T = LowerTriangular(ones(Float64, (3,3)));

julia> T'*D
3×3 UpperTriangular{Float64,Array{Float64,2}}:
 0.5  0.5  0.5
  ⋅   0.5  0.5
  ⋅    ⋅   0.5

julia> D*T'
3×3 Array{Float64,2}:
 0.5  0.5  0.5
 0.0  0.5  0.5
 0.0  0.0  0.5

@andreasnoack
Copy link
Member

The approach looks good to me.

@goggle
Copy link
Contributor Author

goggle commented Sep 20, 2019

I have now adopted the remaining multiplication methods in a similar fashion.

There are still some issues:

julia> using LinearAlgebra

julia> D = Diagonal([0.5, 0.5, 0.5]);

julia> L = LowerTriangular(ones(Int, (3,3)));

julia> L'*D
3×3 Array{Float64,2}:
 0.5  0.5  0.5
 0.0  0.5  0.5
 0.0  0.0  0.5

L'*D should be an UpperTriangular, right?

This would be solved, if similar(L') would return a UpperTriangular, which it currently doesn't:

julia> typeof(similar(L'))
Array{Int64,2}

So, should we add such similar methods? If so, where is the right place to do this?

@andreasnoack
Copy link
Member

So, should we add such similar methods? If so, where is the right place to do this?

Yes. I think it would make sense with such methods. I would think they should go in triangular.jl, if possible.

dkarrasch
dkarrasch reviewed Sep 20, 2019
View reviewed changes
stdlib/LinearAlgebra/src/triangular.jl
similar(A::Transpose{TI,TV}, ::Type{T}) where{T,TI,TV<:$t} = adjtypes[$t](similar(parent(parent(A)), T))
end
end

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I'm not mistaken, then writing this out explicitly yields 8 lines of clear code, as opposed to, well, what you have. 😛

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That's correct, I did not think about this...

@dkarrasch dkarrasch added the domain:linear algebra Linear algebra label Sep 20, 2019
@goggle
Copy link
Contributor Author

goggle commented Oct 9, 2019

I have modified the code acccording to @dkarrasch's suggestion, so this PR is ready to be reviewed.

@dkarrasch dkarrasch closed this Oct 10, 2019
@dkarrasch dkarrasch reopened this Oct 10, 2019
@dkarrasch dkarrasch merged commit 2f7f0bf into JuliaLang:master Oct 10, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@dkarrasch dkarrasch dkarrasch left review comments

Assignees
No one assigned
Labels
domain:linear algebra Linear algebra
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

None yet
3 participants
@goggle @andreasnoack @dkarrasch