Work with Q in qr() is expensive.

[weiyang2048]

I have posted a question on this earlier on discourse

In the picture below, we can see that it is very efficient to generate a random unitary matrix Q using the qr() function in LinearAlgebra.

However, if we want to work with the generated Q, it appears to be very expensive in memory and time. With size 401 by 401, Q*1 takes 24 seconds and 1.5 GiB of memory, with size 500 by 500 it appears to take forever.

[antoine-levitt]

MWE showing the fallback used:

julia> @which qr(randn(2,2)).Q * 1
*(A::AbstractArray, B::Number) in Base at arraymath.jl:55

[carstenbauer]

Can confirm on 1.6 but seems to be fixed on master:

julia> for i in 1:100:500

           A = randn(i,i)
           println("Size = $i")
           printstyled("Time took for q : ", color=:green)
           @time Q = qr(A).Q
           printstyled("Time took for Q*1 : ", color=:red)
           @time Q = Q*1
           println("")
       end
Size = 1
Time took for q :   0.000008 seconds (5 allocations: 352 bytes)
Time took for Q*1 :   0.000010 seconds (2 allocations: 192 bytes)

Size = 101
Time took for q :   0.000418 seconds (8 allocations: 136.922 KiB)
Time took for Q*1 :   0.000158 seconds (4 allocations: 159.656 KiB)

Size = 201
Time took for q :   0.000888 seconds (8 allocations: 429.109 KiB)
Time took for Q*1 :   0.000433 seconds (4 allocations: 631.531 KiB)

Size = 301
Time took for q :   0.001847 seconds (8 allocations: 877.547 KiB)
Time took for Q*1 :   0.001077 seconds (4 allocations: 1.383 MiB)

Size = 401
Time took for q :   0.003631 seconds (8 allocations: 1.447 MiB)
Time took for Q*1 :   0.003361 seconds (4 allocations: 2.454 MiB)

julia> @which qr(randn(2,2)).Q * 1
*(A::LinearAlgebra.AbstractQ, b::Number) in LinearAlgebra at /Users/crstnbr/repos/julia/usr/share/julia/stdlib/v1.7/LinearAlgebra/src/qr.jl:618

julia> versioninfo()
Julia Version 1.7.0-DEV.915
Commit 2ba139c48e* (2021-04-13 04:20 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin20.3.0)
  CPU: Intel(R) Core(TM) i7-9750H CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
Environment:
  JULIA_EDITOR = subl
  JULIA_NUM_THREADS = 6
  JULIA_PKG_SERVER = eu-central.pkg.julialang.org

[antoine-levitt]

Was fixed in #39533

[oscardssmith]

The one thing that I think remains to be fixed is Q[a:b, c:d] which should first get Matrix(Q) to prevent Q from being calculated for each index.

[dkarrasch]

The one thing that I think remains to be fixed is Q[a:b, c:d] which should first get Matrix(Q) to prevent Q from being calculated for each index.

Just for the record. The following script

using BenchmarkTools, LinearAlgebra

Q = qr(randn(500, 300)).Q
@btime Q[50:200, 50:100];

yields

• 1.248 s (23105 allocations: 79.96 MiB) on v1.6.7
• 1.125 s (23105 allocations: 79.96 MiB) on v1.7.3
• 4.078 ms (6 allocations: 3.11 MiB) on v1.8.5
• 4.112 ms (6 allocations: 3.11 MiB) on v1.9
• 562.354 μs (6 allocations: 379.14 KiB) on master.