RFC: Fast mul! with Q matrices from QR factorizations. #31581

nrontsis · 2019-04-02T11:49:32Z

Purpose

Uses LAPACK's gemqrt! and ormqr! when performing mul! with Q matrices, of QR decompositions, or their transpose. Fixes JuliaLang/LinearAlgebra.jl#612. Acknowledged duplicate of #31163.

Tasks list

Add mul! methods
Add tests
Discuss issues with lmul!/mul! on non-strided arrays.

Notes

No methods for mul! with a tranposed AbstractMatrix/AbstractVector input. This is because I could not find a way to perform this without extra memory allocation.

Use BLAS' gemqrt! and ormqr! when performing mul! with Q matrices or their transpose.

andreasnoack · 2019-04-05T06:53:50Z

stdlib/LinearAlgebra/src/qr.jl

+end
+
+function mul!(C::StridedVecOrMat{T}, Q::AbstractQ{T}, B::StridedVecOrMat{T}) where T<:BlasFloat
+    check_dimensions(C, 'N', 'N', Q, B)


Is it necessary to check the dimensions here? Both copyto! and two-argument lmul! have checks so I wouldn't think it would be needed here.

With StridedVecOrMats of BlasFloats, lmul! calls LAPACK.gemqrt! or LAPACK.ormqr! directly. So it depends on whether these assume correct sizes already.

On a different note, there is also function lmul!(A::QRPackedQ, B::AbstractVecOrMat). Should we then also have the proposed mul! methods here for abstract vectors or matrices?

LAPACK.gemqrt! and LAPACK.ormqr! seem to check dimensions, but copyto!(A, B) only checks that length(A) >= length(B) which is less strict than the checks I have.

Yes, thanks, I think that we should also support AbstractVecOrMat. On a related note, I cannot understand thought why rmul! requires StridedMatrix while lmul! does not, i.e.

function lmul!(A::QRPackedQ, B::AbstractVecOrMat) function rmul!(A::StridedMatrix, adjQ::Adjoint{<:Any,<:QRPackedQ})

If one wanted to have stricter checks, one could use copy! instead of copyto!, which requires equal axes. I couldn't find any use of copy! in qr.jl though. That sounds like there's no need for additional size checks then?

Hm, it occurred to me that for mutating *mul! functions you need memory to write into. So it could be that the destination arrays need to be strided, but for 3-arg mul!, the source array could be abstract. Does that make sense for a general rule?

@dkarrasch thanks I wasn't aware of copy!. I will use that and remove the additional size checks.

Well, I thought that since both the LAPACK functions and the generic Julian *mul! functions do size checks, it would suffice to leave it to them and to use the more "generous" copyto! (as suggested by @andreasnoack), independently of whether you go the LAPACK or the generic Julia route. Basically, leave everything as you have, but kick out size checks.

nrontsis · 2019-04-08T09:45:46Z

As suggested, I removed the size checks. Also, following @dkarrasch suggestion, I relaxed the type of the "input" of mul! (i.e. its second or third argument) to be AbstractVecOrMat.

@dkarrasch I must say though that I do not understand what you meant by the following:

for mutating *mul! functions you need memory to write into.

On a related note, I do not understand why Julia's version (non-LAPACK) of rmul! requires strided matrices while lmul! does not.

dkarrasch · 2019-04-08T10:15:30Z

What I meant was that for mutating versions, you require memory to write the result into. I assume that this is why you need a StridedVecOrMat for the output. For an input, it should be okay to have, say, a Range of appropriate size, it's just that you won't be able to call BLAS in that case, but then we still have the generic multiplication versions. From what I understand, strided arrays are arrays that have their data written out one way or another in the memory, as opposed to arrays which get their data, for instance, from function calls, like in ranges. For that reason, I guess both lmul! and rmul! should have a strided object as the one that is mutated, but no-one noticed because it's not tested? Could try to require the mutated arguments to be all strided, and then see if tests pass. Or wait until an expert tells what to do. 😃

dkarrasch · 2019-04-08T10:43:06Z

Consider the following example:

A = rand(4,4)
F = qr(A, Val(true))
F.Q * (1.0:4.0) # works fine
lmul!(F.Q, 1.0:4.0) # yields ERROR: setindex! not defined for UnitRange{Int64}
lmul!(F.Q, collect(1.0:4.0)) # works fine again, Float64.(1:4) gives a StridedVector

Same with F = qr(A), where F.Q::QRCompactWYQ. So I very much assume that the mutated vec or mat should always be strided, but the input can be abstract, unless we call LAPACK. @andreasnoack ?

nrontsis · 2019-04-08T16:22:49Z

Thanks alot @dkarrasch for the explanations!

Following the suggestions above I relaxed the definitions of the newly-added methods to be:

mul!(C::StridedVecOrMat{T}, Q::AbstractQ{T}, B::AbstractVecOrMat{T}) where {T}
mul!(C::StridedVecOrMat{T}, A::AbstractVecOrMat{T}, Q::AbstractQ{T}) where {T}
mul!(C::StridedVecOrMat{T}, adjQ::Adjoint{<:Any,<:AbstractQ{T}}, B::AbstractVecOrMat{T}) where {T}
mul!(C::StridedVecOrMat{T}, A::AbstractVecOrMat{T}, adjQ::Adjoint{<:Any,<:AbstractQ{T}}) where {T}

However, the LinearAlgebra/test/matmul.jl now fails on line 508 with the following stacktrace:

Test Failed at /Users/nrontsis/julia/stdlib/LinearAlgebra/test/ambiguous_exec.jl:4
  Expression: detect_ambiguities(LinearAlgebra; imported=true, recursive=true) == []
   Evaluated: Tuple{Method,Method}[(mul!(C::AbstractArray{T,2} where T, ...  in LinearAlgebra at /Users/nrontsis/julia/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/qr.jl:742)] == Any[]
ERROR: Error while loading expression starting at /Users/nrontsis/julia/stdlib/LinearAlgebra/test/ambiguous_exec.jl:4
caused by [exception 1]
There was an error during testing
method ambiguity: Test Failed at /Users/nrontsis/julia/stdlib/LinearAlgebra/test/matmul.jl:508

Any advice would be greatly appreciated.

dkarrasch · 2019-04-08T17:26:03Z

I have seen elsewhere in the code (or remember Andreas mention) that sometimes methods are split for vectors and matrices, like

mul!(C::StridedMatrix{T}, Q::AbstractQ{T}, B::AbstractMatrix{T}) where {T}
mul!(C::StridedVector{T}, Q::AbstractQ{T}, B::AbstractVector{T}) where {T}

etc., maybe that helps?

EDIT: You generalized the method signature in commit 3dd2729, and didn't mention the issue before. So that could be an indication that this is indeed the issue.

nrontsis · 2019-04-11T17:09:42Z

@dkarrasch thanks for the suggestion. Unfortunately it seems that splitting the methods to vectors and matrices did not resolve the method ambiguity.

However, restricting all the inputs to be strided solved it. This is obviously not ideal, as, according to the discussion above, the second or third input to mul! should be able to be non-strided.

Regardless, I restricted the method definitions to only allow for strided inputs, so as to have a reference commit at which the tests pass.

I have also updated the description of the issue to include a task list, hoping that this will facilitate the discussion.

dkarrasch · 2019-04-11T18:23:54Z

Good to know that there was a good reason for the strided requirement. After all, I think this should cover the regular use cases. Vectors that result from other multiplication processes will be likely strided anyway, and the Range example above was very artificial. I think this is a good solution. Probably, if you really want to, you can still pass a non-strided vector, it's just that it will go through the (slow) fallback.

nrontsis · 2019-04-11T20:08:43Z

Okay, in that case lets also restrict lmul! to strided arrays? If you agree with that then I will make the change and, in my point of view, I think we would be ready to merge to master?

dkarrasch · 2019-04-12T06:51:32Z

I'm not sure if we should add restrictions without immediate need. You could try to relax the rmul! method, because both the r/lmul! can be called without going through mul!.

nrontsis · 2019-04-12T14:29:53Z

Okay let's leave it as it is then. Do you think this is ready to be merged?

dkarrasch · 2019-04-12T14:38:28Z

I don't have authority to merge, but apparently I do have authority to approve, which I have just done.

nrontsis added 2 commits April 2, 2019 12:43

Fast mul! with Q matrices from QR factorizations.

cbf93c3

Use BLAS' gemqrt! and ormqr! when performing mul! with Q matrices or their transpose.

Fixing tests.

f656fca

ViralBShah requested review from andreasnoack and Sacha0 April 2, 2019 13:02

ViralBShah added the linear algebra Linear algebra label Apr 2, 2019

andreasnoack reviewed Apr 5, 2019

View reviewed changes

nrontsis added 2 commits April 8, 2019 10:38

Removing redundant size checks; relaxing mul! source type

8b7d133

Relax AbstractMatrix to AbstractVecOrMat

3dd2729

Relaxing eltype; removing ambiguous DimensionMismatch test.

12f92a8

Restrict mul! to strided inputs to avoid method ambiguity.

9ec2ba0

nrontsis changed the title ~~Fast mul! with Q matrices from QR factorizations.~~ RFC: Fast mul! with Q matrices from QR factorizations. Apr 12, 2019

dkarrasch approved these changes Apr 12, 2019

View reviewed changes

andreasnoack merged commit 4d5a901 into JuliaLang:master Apr 15, 2019

nrontsis mentioned this pull request Apr 15, 2019

fixed the issue with the performance of mul! on matrix type QRCompact… #31163

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

RFC: Fast mul! with Q matrices from QR factorizations. #31581

RFC: Fast mul! with Q matrices from QR factorizations. #31581

nrontsis commented Apr 2, 2019 •

edited

Loading

andreasnoack Apr 5, 2019

dkarrasch Apr 5, 2019

nrontsis Apr 5, 2019 •

edited

Loading

nrontsis Apr 5, 2019

dkarrasch Apr 5, 2019

dkarrasch Apr 5, 2019

nrontsis Apr 5, 2019

dkarrasch Apr 5, 2019

nrontsis commented Apr 8, 2019 •

edited

Loading

dkarrasch commented Apr 8, 2019

dkarrasch commented Apr 8, 2019

nrontsis commented Apr 8, 2019 •

edited

Loading

dkarrasch commented Apr 8, 2019 •

edited

Loading

nrontsis commented Apr 11, 2019 •

edited

Loading

dkarrasch commented Apr 11, 2019

nrontsis commented Apr 11, 2019

dkarrasch commented Apr 12, 2019

nrontsis commented Apr 12, 2019

dkarrasch commented Apr 12, 2019

RFC: Fast mul! with Q matrices from QR factorizations. #31581

RFC: Fast mul! with Q matrices from QR factorizations. #31581

Conversation

nrontsis commented Apr 2, 2019 • edited Loading

Purpose

Tasks list

Notes

andreasnoack Apr 5, 2019

Choose a reason for hiding this comment

dkarrasch Apr 5, 2019

Choose a reason for hiding this comment

nrontsis Apr 5, 2019 • edited Loading

Choose a reason for hiding this comment

nrontsis Apr 5, 2019

Choose a reason for hiding this comment

dkarrasch Apr 5, 2019

Choose a reason for hiding this comment

dkarrasch Apr 5, 2019

Choose a reason for hiding this comment

nrontsis Apr 5, 2019

Choose a reason for hiding this comment

dkarrasch Apr 5, 2019

Choose a reason for hiding this comment

nrontsis commented Apr 8, 2019 • edited Loading

dkarrasch commented Apr 8, 2019

dkarrasch commented Apr 8, 2019

nrontsis commented Apr 8, 2019 • edited Loading

dkarrasch commented Apr 8, 2019 • edited Loading

nrontsis commented Apr 11, 2019 • edited Loading

dkarrasch commented Apr 11, 2019

nrontsis commented Apr 11, 2019

dkarrasch commented Apr 12, 2019

nrontsis commented Apr 12, 2019

dkarrasch commented Apr 12, 2019

nrontsis commented Apr 2, 2019 •

edited

Loading

nrontsis Apr 5, 2019 •

edited

Loading

nrontsis commented Apr 8, 2019 •

edited

Loading

nrontsis commented Apr 8, 2019 •

edited

Loading

dkarrasch commented Apr 8, 2019 •

edited

Loading

nrontsis commented Apr 11, 2019 •

edited

Loading