# Copying a QR decomposition matrix Q is extremely slow

**URL:** <https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430>\
**Category:** Performance\
**Tags:** question, linearalgebra\
**Created:** [November 19, 2020, 12:51pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430 "2020-11-19T12:51:24Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [November 19, 2020, 12:51pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/1 "2020-11-19T12:51:24Z")

</div>

This is on julia 1.5.3

```julia
julia> f = qr(ones(100, 20), Val(true));

julia> copy(f.Q); @time copy(f.Q);
  0.057310 seconds (40.00 k allocations: 338.516 MiB, 10.85% gc time)

julia> copy(Array(f.Q)); @time copy(Array(f.Q));
  0.000074 seconds (7 allocations: 84.969 KiB)

julia> typeof(copy(f.Q)) == typeof(copy(Array(f.Q)))
true

```

The `copy` appears to be called in `svdvals` which slows to a crawl for large matrices. Is this a known issue?

It appears that `copy(f.Q)` relies on a default method for `AbstractArray`s, which might explain the difference in performance. Is this something that may be improved?

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [November 19, 2020, 1:29pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/2 "2020-11-19T13:29:06Z")

</div>

Upon further introspection it appears that indexing into `Q` is slow.

```julia
function getindex(Q::AbstractQ, i::Integer, j::Integer)
    x = zeros(eltype(Q), size(Q, 1))
    x[i] = 1
    y = zeros(eltype(Q), size(Q, 2))
    y[j] = 1
    return dot(x, lmul!(Q, y))
end

```

results in

```julia
julia> @btime $f.Q[1, 1];
  22.795 μs (6 allocations: 48.84 KiB)

```

As a consequence

```julia
julia> @time f.Q[:,:];
 30.861779 seconds (4.00 M allocations: 46.529 GiB, 2.50% gc time)

```

---

<div class="post-metadata">

**Author:** ![malacroi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/malacroi/32/19745_2.png) [@malacroi](https://discourse.julialang.org/u/malacroi)\
**Post date:** [November 19, 2020, 1:39pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/3 "2020-11-19T13:39:44Z")

</div>

I’m definitely not an expert here, but is that because the`q` isn’t present as a dense matrix until you copy it. It’s stored in packed form as a product of Householder reflectors. Do you actually a dense copy? You can multiply by it while still in packed form.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 19, 2020, 2:00pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/4 "2020-11-19T14:00:40Z")

</div>

@malacroi is correct. Extracting Q from the data structure you get from qr or qr! is expensive. IIRC, it’s O(N^3) work. If you just want to multiply D by Q’, for example, then `f.Q'*B` is what you want.

For example

```julia
julia> A=rand(1000,1000); AC=copy(A); AF=qr!(A);

julia> C=AF.Q*AF.R; norm(C-AC)
4.49007e-13

julia> D=AF.Q'*AC; RC=copy(AF.R); norm(D-RC)
1.27313e-13

```

The copy of AF.R is cheap because you don’t have to compute anything new. Not so if you do AF.Q.

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [November 19, 2020, 2:27pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/5 "2020-11-19T14:27:45Z")

</div>

The fast way to extract the `Q` factor as a dense matrix is just what you already do in your first example, namely `Array(f.Q)`. Indeed, `copy` will probably fall back on a general method that uses indexing, which is slow for the way `f.Q` is stored/represented.

If you want to make a copy of an array or some object that behaves as it, but is maybe not explicit (like a Q factor, or a transposed matrix), and you want the result to be a explicit `Array` that has fast indexing and its own memory, then `Array(object)` is a good way.

---

<div class="post-metadata">

**Author:** ![nalimilan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nalimilan/32/147_2.png) [@nalimilan](https://discourse.julialang.org/u/nalimilan)\
**Post date:** [November 19, 2020, 2:33pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/6 "2020-11-19T14:33:43Z")

</div>

Looks like we’re missing a specialized `copy(x::QRPackedQ) = Array(x)` method for performance. But maybe that method should even return a `QRPackedQ` object, since `copy` doesn’t have to return an `Array`. You could file an issue in GitHub.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 19, 2020, 2:37pm UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/7 "2020-11-19T14:37:50Z")

</div>

`C=Array(AF.Q)` is ok if you **must** have Q. It is far better than copy. The allocation burden seems a bit steep for the 8MB you need to store Q.

```julia
julia> A=rand(1000,1000); AF=qr!(A);

julia> @btime C=Array($AF.Q);
  31.310 ms (8 allocations: 22.89 MiB)

```

---

<div class="post-metadata">

**Author:** ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)\
**Post date:** [November 22, 2020, 8:58am UTC](https://discourse.julialang.org/t/copying-a-qr-decomposition-matrix-q-is-extremely-slow/50430/8 "2020-11-22T08:58:40Z")

</div>

Currently `svdvals` relies on a copy.

```julia
svdvals(A::AbstractMatrix{<:BlasFloat}) = svdvals!(copy(A))

```

I realize that I may evaluate `svdvals(Array(Q))` instead of `svdvals(Q)`, however is it reasonable to expect `svdvals` to perform this by itself? This is a different question from whether we should implement a specific method for `copy(::AbstractQ)`, as — from the suggestions above — `copy` might not be be the best choice for certain types of matrices.