Unexpected performance loss

[andreasvarga]

The computation of the orthogonal transformation matrix Q in the QR-decomposion is occasionally necessary. I tried the following two ways to recover Q from the computed factorization F = qr(A) of a matrix A of order n:

1. Q = F.Q*Matrix(I,n,n)
   and
2. Q = F.Q*I

and compared the timing results. Before the tests, the second way appeared to me being the more elegant approach for this computation.

The timing results have been obtained with Julia 1.3.1 (see Julia version information below) with the following code:

using LinearAlgebra
n = 250; 
A = rand(n,n); F = qr(A); 
@time Q = F.Q*Matrix(I,n,n);
@time Q1 = F.Q*I;
Q == Q1
Q ≈ Q1
norm(Q-Q1)

For n = 250 the results are:

julia> @time Q = F.Q*Matrix(I,n,n);
  0.003613 seconds (11 allocations: 1.014 MiB)

julia> @time Q1 = F.Q*I;
  4.305979 seconds (187.51 k allocations: 389.576 MiB, 0.79% gc time)

julia> Q == Q1
true

julia> Q ≈ Q1
true

julia> norm(Q-Q1)
0.0

For n = 500 the results are:

julia> @time Q = F.Q*Matrix(I,n,n);
  0.011218 seconds (11 allocations: 4.054 MiB)

julia> @time Q1 = F.Q*I;
 61.151209 seconds (750.01 k allocations: 2.908 GiB, 0.24% gc time)

julia> Q == Q1
false

julia> Q ≈ Q1
true

julia> norm(Q-Q1)
9.810156612825285e-15

For n = 1000 the results are:

julia> @time Q = F.Q*Matrix(I,n,n);
  0.065219 seconds (11 allocations: 16.213 MiB)

julia> @time Q1 = F.Q*I;
969.663585 seconds (3.00 M allocations: 22.717 GiB, 0.10% gc time)

julia> Q == Q1
false

julia> Q ≈ Q1
true

julia> norm(Q-Q1)
2.1475440222104054e-14

So, the "elegant" way is about 1200 times slower for n = 250, about 5450 times slower for n = 500 and about 14900 times slower for n = 1000.

I think the severe performace loss can be easily avoided (at least for this particular case), by explicitly addressing the case with the uniform scaling.

julia> versioninfo()
Julia Version 1.3.1
Commit 2d5741174c (2019-12-30 21:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-3770S CPU @ 3.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.1 (ORCJIT, ivybridge)
Environment:
  JULIA_EDITOR = "C:\Users\Andreas\AppData\Local\Programs\Microsoft VS Code\Code.exe"
  JULIA_NUM_THREADS =

[oscardssmith]

This problem seems to occur from missing routines for multiplying LinearAlgebra.QRCompactWYQ * Number. The fallback method for AbstractArray*Number uses broadcasting, so ends up indexing into the factorization which is relatively expensive. Defining

function (*)(A::LinearAlgebra.AbstractQ, I::UniformScaling)
    return Matrix(A) .* I.λ
end

fixes the problem. However, the easiest way around this problem is to just use Matrix(F.Q) instead of multiplying by I

[oscardssmith]

I've opened a PR that fixes it by adding this method.

[andreasvarga]

Nice fix. Thanks for the hint to use Matrix(F.Q).

[andreasvarga]

I think the above function not always does what I intended to compute, i.e., to explicitly determine the full transformation matrix Q from the QR-decompositions. The following example shows what I meant:

julia> using LinearAlgebra

julia> import Base: *

julia> function (*)(A::LinearAlgebra.AbstractQ, I::UniformScaling)
           return Matrix(A) .* I.λ
       end
* (generic function with 393 methods)

julia> A = rand(2,1)
2×1 Array{Float64,2}:
 0.8927617960163903
 0.6928997896223985

julia> F = qr(A)
LinearAlgebra.QRCompactWY{Float64,Array{Float64,2}}
Q factor:
2×2 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.789982  -0.61313 
 -0.61313    0.789982
R factor:
1×1 Array{Float64,2}:
 -1.130103421322657

julia> F.Q*I     # this is only the first column of Q
2×1 Array{Float64,2}:
 -0.7899823849497905
 -0.6131295388978103

julia> Matrix(F.Q)  # this is only the first column of Q
2×1 Array{Float64,2}:
 -0.7899823849497905
 -0.6131295388978103

julia> F.Q*Matrix(I,2,2)  # this is what I wanted to determine
2×2 Array{Float64,2}:
 -0.789982  -0.61313
 -0.61313    0.789982

I wonder if the functionality of Matrix(F.Q) is the intended one.

[ViralBShah]

julia> n=1000
1000

julia> A = rand(n,n); F = qr(A);

julia> @time Q = F.Q*Matrix(I,n,n);
  0.061705 seconds (7 allocations: 16.213 MiB, 9.72% gc time)

julia> @time Q1 = F.Q*I;
  0.052477 seconds (5 allocations: 15.259 MiB)