Added function to compute Diagonal*F.Q' and faster methods for other variants. #41248

ArunS-tack · 2021-06-16T19:39:40Z

Requesting reviews from @dkarrasch @andreasnoack .

dkarrasch · 2021-06-17T08:50:08Z

I don't think this is the right approach. I think we should add

*(A::Diagonal, adjB::Adjoint{<:Any,<:Union{LinearAlgebra.QRCompactWYQ,LinearAlgebra.QRPackedQ}}) =
    rmul!(copyto!(Array{promote_type(eltype(A), eltype(adjB))}(undef, size(A)...), A), adjB)
*(A::BiTriSym, adjB::Adjoint{<:Any,<:Union{LinearAlgebra.QRCompactWYQ,LinearAlgebra.QRPackedQ}}) =
    rmul!(copyto!(Array{promote_type(eltype(A), eltype(adjB))}(undef, size(A)...), A), adjB)
*(adjA::Adjoint{<:Any,<:Union{LinearAlgebra.QRCompactWYQ,LinearAlgebra.QRPackedQ}}, B::Diagonal) =
    lmul!(adjA, copyto!(Array{promote_type(eltype(adjA), eltype(B))}(undef, size(B)...), B))
*(adjA::Adjoint{<:Any,<:Union{LinearAlgebra.QRCompactWYQ,LinearAlgebra.QRPackedQ}}, B::BiTriSym) =
    lmul!(adjA, copyto!(Array{promote_type(eltype(adjA), eltype(B))}(undef, size(B)...), B))

below

julia/stdlib/LinearAlgebra/src/special.jl

Lines 283 to 286 in 25efa99

    
           rmul!(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = 
        
               (B = adjB.parent; rmul!(full!(A), adjoint(B))) 
        
           *(A::AbstractTriangular, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = 
        
               (B = adjB.parent; *(copyto!(similar(parent(A)), A), adjoint(B)))

The methods that currently handle Q' * D treat Q like a quickly indexable matrix (in adjoint!), which is likely slow.

dkarrasch · 2021-06-17T09:46:09Z

Even simpler: we already have

julia/stdlib/LinearAlgebra/src/diagonal.jl

Line 292 in 6dfa690

    
           *(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; rmul!(Array(D), adjoint(Q)))

which should be

*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
    rmul!(Array{promote_type(eltype(D), eltype(adjQ))}(D), adjQ)

but perhaps we take the opportunity to make sure we also cover *(adjQ, D) and the other *diagonal types as suggested above.

ArunS-tack · 2021-06-17T10:08:13Z

Even simpler: we already have

julia/stdlib/LinearAlgebra/src/diagonal.jl

Line 292 in 6dfa690

*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) = (Q = adjQ.parent; rmul!(Array(D), adjoint(Q)))

which should be
*(D::Diagonal, adjQ::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}}) =
    rmul!(Array{promote_type(eltype(D), eltype(adjQ))}(D), adjQ)
but perhaps we take the opportunity to make sure we also cover *(adjQ, D) and the other *diagonal types as suggested above.

Nice! This one seems much simpler. But to make it clear, we are not implementing the function as made in the commit because it only covers a narrow number of cases and takes longer to compute the multiplication, right?

dkarrasch · 2021-06-17T10:11:41Z

The function as it is currently in the commit has a very narrow type signature, and it also doesn't mutate its first argument as it should. The idea is to provide [l/r]mul! with the correct eltype'd array in the * method.

ArunS-tack · 2021-06-17T10:33:28Z

function rmul!(A::StridedVecOrMat, adjB::Adjoint{<:Any,<:Union{QRCompactWYQ,QRPackedQ}})
    return rmul!(Array{promote_type(eltype(A), eltype(adjB))}(A), adjB)
end

I modified the code a bit, its also a lot faster now and most probably covers a wide range of type signature too. What are your views about this one? Also can we not use Union{Diagonal, BiTriSym} to cover it all in a single method/function?

dkarrasch · 2021-06-17T10:46:50Z

In your proposal, you again create a copy of A first and then multiply into it. This is not how *mul!s are supposed to work. And this also doesn't cover the case when the first factor is a Diagonal.

Also can we not use Union{Diagonal, BiTriSym} to cover it all in a single method/function?

This will likely lead to method ambiguities.

ArunS-tack · 2021-06-17T11:03:40Z

Got it! I will implement what you suggested above in diagonal.jl along with other methods then.

ArunS-tack · 2021-06-17T12:04:45Z

but perhaps we take the opportunity to make sure we also cover *(adjQ, D)

A method for that already exists I guess, its working fine.

dkarrasch · 2021-06-17T13:38:08Z

A method for that already exists I guess, its working fine.

Sure, but I guess it's slow.

ArunS-tack · 2021-06-17T13:51:32Z

well, actually it seems the methods already exist for all other types i.e BiTriSym and *(adjQ, D) as well for Diagonal.
Only one which doesn't exists is *(D,adjQ) for Diagonal. Shall I anyway add them in diagonal.jl?

ArunS-tack · 2021-06-17T13:54:46Z

A method for that already exists I guess, its working fine.

Sure, but I guess it's slow.

Yes, indeed. It is slow.

dkarrasch · 2021-06-17T13:58:44Z

Right. This can be fixed by my proposals above. I'll explain later why this is. Gotta leave now...

codecov · 2021-06-18T02:49:35Z

Codecov Report

Merging #41248 (6dfa690) into master (004cb25) will increase coverage by 1.82%.
The diff coverage is 84.73%.

❗ Current head 6dfa690 differs from pull request most recent head be45c69. Consider uploading reports for the commit be45c69 to get more accurate results

@@            Coverage Diff             @@
##           master   #41248      +/-   ##
==========================================
+ Coverage   86.15%   87.97%   +1.82%     
==========================================
  Files         350      351       +1     
  Lines       65991    74690    +8699     
==========================================
+ Hits        56854    65711    +8857     
+ Misses       9137     8979     -158

Impacted Files	Coverage Δ
base/atomics.jl	`84.00% <ø> (+0.66%)`	⬆️
base/baseext.jl	`94.73% <ø> (ø)`
base/compiler/typeutils.jl	`82.53% <ø> (+0.17%)`	⬆️
base/compiler/utilities.jl	`87.90% <ø> (+1.24%)`	⬆️
base/compiler/validation.jl	`88.99% <ø> (+0.64%)`	⬆️
base/complex.jl	`99.43% <ø> (-0.16%)`	⬇️
base/condition.jl	`78.00% <ø> (+4.08%)`	⬆️
base/coreio.jl	`68.75% <ø> (-6.25%)`	⬇️
base/cpuid.jl	`82.60% <ø> (ø)`
base/deepcopy.jl	`100.00% <ø> (ø)`
... and 512 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 083272b...be45c69. Read the comment docs.

dkarrasch

There are a few redundant lines, and then we need tests for the new BiTriSym methods. Would be also great to see some benchmarks, to make sure we do the right thing.

stdlib/LinearAlgebra/src/special.jl

stdlib/LinearAlgebra/src/diagonal.jl

Co-authored-by: Daniel Karrasch <[email protected]>

stdlib/LinearAlgebra/test/special.jl

Co-authored-by: Daniel Karrasch <[email protected]>

stdlib/LinearAlgebra/test/qr.jl

ArunS-tack · 2021-06-18T11:14:51Z

About Benchmarks, do we have to add them in some file or shall I just post it here as a comment?

Co-authored-by: Daniel Karrasch <[email protected]>

ArunS-tack · 2021-06-18T11:37:26Z

BENCHMARKS

julia> @benchmark Tri*A
BenchmarkTools.Trial: 
  memory estimate:  320 bytes
  allocs estimate:  2
  --------------
  minimum time:     542.112 ns (0.00% GC)
  median time:      549.644 ns (0.00% GC)
  mean time:        557.678 ns (0.56% GC)
  maximum time:     4.261 μs (86.09% GC)
  --------------
  samples:          10000
  evals/sample:     188

julia> @benchmark A*Tri

BenchmarkTools.Trial: 
  memory estimate:  320 bytes
  allocs estimate:  2
  --------------
  minimum time:     527.632 ns (0.00% GC)
  median time:      535.526 ns (0.00% GC)
  mean time:        543.155 ns (0.59% GC)
  maximum time:     3.431 μs (81.91% GC)
  --------------
  samples:          10000
  evals/sample:     190

julia> 

julia> @benchmark Bi*A
BenchmarkTools.Trial: 
  memory estimate:  320 bytes
  allocs estimate:  2
  --------------
  minimum time:     540.561 ns (0.00% GC)
  median time:      548.280 ns (0.00% GC)
  mean time:        555.362 ns (0.57% GC)
  maximum time:     3.101 μs (81.05% GC)
  --------------
  samples:          10000
  evals/sample:     189

julia> @benchmark A*Bi
BenchmarkTools.Trial: 
  memory estimate:  320 bytes
  allocs estimate:  2
  --------------
  minimum time:     528.511 ns (0.00% GC)
  median time:      535.963 ns (0.00% GC)
  mean time:        542.820 ns (0.58% GC)
  maximum time:     3.123 μs (81.92% GC)
  --------------
  samples:          10000
  evals/sample:     190

julia> @benchmark Sym*A
BenchmarkTools.Trial: 
  memory estimate:  320 bytes
  allocs estimate:  2
  --------------
  minimum time:     526.316 ns (0.00% GC)
  median time:      533.553 ns (0.00% GC)
  mean time:        541.513 ns (0.60% GC)
  maximum time:     3.133 μs (81.88% GC)
  --------------
  samples:          10000
  evals/sample:     190

julia> @benchmark A*Sym 
BenchmarkTools.Trial: 
  memory estimate:  320 bytes
  allocs estimate:  2
  --------------
  minimum time:     515.406 ns (0.00% GC)
  median time:      522.786 ns (0.00% GC)
  mean time:        530.500 ns (0.60% GC)
  maximum time:     3.139 μs (79.75% GC)
  --------------
  samples:          10000
  evals/sample:     192

dkarrasch · 2021-06-18T11:44:51Z

For comparison, what are the times on current master?

ArunS-tack · 2021-06-18T11:48:25Z

julia> @benchmark Bi*A

BenchmarkTools.Trial: 
  memory estimate:  5.31 KiB
  allocs estimate:  50
  --------------
  minimum time:     7.781 μs (0.00% GC)
  median time:      7.979 μs (0.00% GC)
  mean time:        8.122 μs (1.18% GC)
  maximum time:     250.542 μs (96.43% GC)
  --------------
  samples:          10000
  evals/sample:     4

julia> @benchmark A*Bi
BenchmarkTools.Trial: 
  memory estimate:  13.47 KiB
  allocs estimate:  131
  --------------
  minimum time:     20.250 μs (0.00% GC)
  median time:      20.833 μs (0.00% GC)
  mean time:        21.273 μs (1.76% GC)
  maximum time:     1.292 ms (98.03% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> @benchmark Sym*A
BenchmarkTools.Trial: 
  memory estimate:  5.20 KiB
  allocs estimate:  49
  --------------
  minimum time:     7.708 μs (0.00% GC)
  median time:      8.011 μs (0.00% GC)
  mean time:        8.119 μs (0.96% GC)
  maximum time:     273.562 μs (96.52% GC)
  --------------
  samples:          10000
  evals/sample:     4

julia> @benchmark A*Sym 
BenchmarkTools.Trial: 
  memory estimate:  13.36 KiB
  allocs estimate:  130
  --------------
  minimum time:     20.250 μs (0.00% GC)
  median time:      20.875 μs (0.00% GC)
  mean time:        21.530 μs (1.88% GC)
  maximum time:     1.479 ms (98.16% GC)
  --------------
  samples:          10000
  evals/sample:     1

ArunS-tack · 2021-06-18T17:03:54Z

Not sure why is it still showing 'Changes requested'? Afaik, I have resolved all the reviews.

…g#41248) Co-authored-by: Daniel Karrasch <[email protected]>

ArunS-tack added 6 commits June 17, 2021 00:50

added function.

5259712

added test case.

d92c146

Update qr.jl

3a348b6

removed white spaces.

6ed21b5

removed white spaces

958279b

Update qr.jl

bf5c9c6

dkarrasch added the domain:linear algebra Linear algebra label Jun 17, 2021

ArunS-tack added 3 commits June 17, 2021 19:45

Update diagonal.jl

8678398

Update qr.jl

bd6f306

Update special.jl

4d7335c

ArunS-tack changed the title ~~Added function to compute I(n)*F.Q' [element promotion of I(n)].~~ Added function to compute Diagonal*F.Q' and faster methods for other variants. Jun 17, 2021

Merge branch 'master' into bool

48525de

dkarrasch requested changes Jun 18, 2021

View reviewed changes

stdlib/LinearAlgebra/src/special.jl Outdated Show resolved Hide resolved

stdlib/LinearAlgebra/src/diagonal.jl Outdated Show resolved Hide resolved

ArunS-tack and others added 4 commits June 18, 2021 14:15

Update stdlib/LinearAlgebra/src/diagonal.jl

d4c4cd7

Co-authored-by: Daniel Karrasch <[email protected]>

Update stdlib/LinearAlgebra/src/special.jl

cc7fca7

Co-authored-by: Daniel Karrasch <[email protected]>

Merge branch 'JuliaLang:master' into bool

ed9caa2

testcases for BiTriSym.

f4c3a7d

dkarrasch reviewed Jun 18, 2021

View reviewed changes

stdlib/LinearAlgebra/test/special.jl Outdated Show resolved Hide resolved

Update stdlib/LinearAlgebra/test/special.jl

3f98e5a

Co-authored-by: Daniel Karrasch <[email protected]>

dkarrasch reviewed Jun 18, 2021

View reviewed changes

stdlib/LinearAlgebra/test/qr.jl Show resolved Hide resolved

Update stdlib/LinearAlgebra/test/qr.jl

be45c69

Co-authored-by: Daniel Karrasch <[email protected]>

ArunS-tack requested a review from dkarrasch June 18, 2021 17:10

dkarrasch approved these changes Jun 18, 2021

View reviewed changes

dkarrasch merged commit 3279e20 into JuliaLang:master Jun 18, 2021

ArunS-tack deleted the bool branch June 18, 2021 17:56

dkarrasch mentioned this pull request Jun 18, 2021

make AbstractArray{T}(...) constructor in LinearAlgebra more consistent #40831

Merged

johanmon pushed a commit to johanmon/julia that referenced this pull request Jul 5, 2021

Added multiplication functions for Ajoint{Q} and *diagonals (JuliaLan…

5155db1

…g#41248) Co-authored-by: Daniel Karrasch <[email protected]>

jessymilare pushed a commit to jessymilare/julia that referenced this pull request Sep 11, 2021

Added multiplication functions for Ajoint{Q} and *diagonals (JuliaLan…

c612015

…g#41248) Co-authored-by: Daniel Karrasch <[email protected]>

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Added function to compute Diagonal*F.Q' and faster methods for other variants. #41248

Added function to compute Diagonal*F.Q' and faster methods for other variants. #41248

ArunS-tack commented Jun 16, 2021 •

edited

Loading

dkarrasch commented Jun 17, 2021

dkarrasch commented Jun 17, 2021 •

edited

Loading

ArunS-tack commented Jun 17, 2021

dkarrasch commented Jun 17, 2021

ArunS-tack commented Jun 17, 2021 •

edited

Loading

dkarrasch commented Jun 17, 2021 •

edited

Loading

ArunS-tack commented Jun 17, 2021

ArunS-tack commented Jun 17, 2021

dkarrasch commented Jun 17, 2021

ArunS-tack commented Jun 17, 2021 •

edited

Loading

ArunS-tack commented Jun 17, 2021

dkarrasch commented Jun 17, 2021

codecov bot commented Jun 18, 2021 •

edited

Loading

dkarrasch left a comment

ArunS-tack commented Jun 18, 2021

ArunS-tack commented Jun 18, 2021

dkarrasch commented Jun 18, 2021

ArunS-tack commented Jun 18, 2021

ArunS-tack commented Jun 18, 2021

Added function to compute Diagonal*F.Q' and faster methods for other variants. #41248

Added function to compute Diagonal*F.Q' and faster methods for other variants. #41248

Conversation

ArunS-tack commented Jun 16, 2021 • edited Loading

dkarrasch commented Jun 17, 2021

dkarrasch commented Jun 17, 2021 • edited Loading

ArunS-tack commented Jun 17, 2021

dkarrasch commented Jun 17, 2021

ArunS-tack commented Jun 17, 2021 • edited Loading

dkarrasch commented Jun 17, 2021 • edited Loading

ArunS-tack commented Jun 17, 2021

ArunS-tack commented Jun 17, 2021

dkarrasch commented Jun 17, 2021

ArunS-tack commented Jun 17, 2021 • edited Loading

ArunS-tack commented Jun 17, 2021

dkarrasch commented Jun 17, 2021

codecov bot commented Jun 18, 2021 • edited Loading

Codecov Report

dkarrasch left a comment

Choose a reason for hiding this comment

ArunS-tack commented Jun 18, 2021

ArunS-tack commented Jun 18, 2021

dkarrasch commented Jun 18, 2021

ArunS-tack commented Jun 18, 2021

ArunS-tack commented Jun 18, 2021

ArunS-tack commented Jun 16, 2021 •

edited

Loading

dkarrasch commented Jun 17, 2021 •

edited

Loading

ArunS-tack commented Jun 17, 2021 •

edited

Loading

dkarrasch commented Jun 17, 2021 •

edited

Loading

ArunS-tack commented Jun 17, 2021 •

edited

Loading

codecov bot commented Jun 18, 2021 •

edited

Loading