# Extracting Q from qrfact allocation-free

**URL:** <https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083>\
**Category:** General Usage\
**Created:** [March 31, 2018, 5:09am UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083 "2018-03-31T05:09:52Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [March 31, 2018, 5:09am UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/1 "2018-03-31T05:09:52Z")

</div>

I am trying to extract the `Q` matrix from a QR decomposition without allocating. The following hack doesn’t work:

```julia
julia> using LinearAlgebra

julia> A = rand(3,2);

julia> qrf = qrfact(A);

julia> qrf.Q * Matrix{Float64}(I,2,2)
3×2 Array{Float64,2}:
 -0.554696 0.791229
 -0.64075 -0.603576
 -0.530803 -0.0982461

julia> Q = zeros(3,2);

julia> A_mul_B!(Q, qrf.Q, Matrix{Float64}(I,2,2))
┌ Warning: `A_mul_B!(C::AbstractMatrix, A::AbstractVecOrMat, B::AbstractVecOrMat)` is deprecated, use `mul!(C, A, B)` instead.
│ caller = top-level scope
└ @ Core :0
ERROR: DimensionMismatch("matrix A has dimensions (3,3), matrix B has dimensions (2,2)")
Stacktrace:
 [1] _generic_matmatmul!(::Array{Float64,2}, ::Char, ::Char, ::LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}, ::Array{Float64,2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\site\v0.7\LinearAlgebra\src\matmul.jl:518
 [2] generic_matmatmul!(::Array{Float64,2}, ::Char, ::Char, ::LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}, ::Array{Float64,2}) at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\site\v0.7\LinearAlgebra\src\matmul.jl:509
 [3] mul! at C:\cygwin\home\Administrator\buildbot\worker\package_win64\build\usr\share\julia\site\v0.7\LinearAlgebra\src\matmul.jl:156 [inlined]
 [4] A_mul_B!(::Array{Float64,2}, ::LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}, ::Array{Float64,2}) at .\deprecated.jl:57
 [5] top-level scope

```

Any idea how to do this in v0.7 and preferably v0.6?

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [March 31, 2018, 6:14am UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/2 "2018-03-31T06:14:39Z")

</div>

In my opinion,

```julia
copyto!(Q, CartesianIndices(Q), qrf.Q, CartesianIndices(Q))

```

should do what you want and not allocate, however

```julia
julia> using BenchmarkTools; @btime copyto!($Q, $(CartesianIndices(Q)), $(qrf.Q), $(CartesianIndices(Q)))
  5.764 μs (18 allocations: 1.88 KiB)
3×2 Array{Float64,2}:
 -0.945825 0.245839
 -0.323885 -0.668782
 -0.0226788 -0.701637

```

Honestly, I find `QRCompactWYQ` (`qrf.Q isa QRCompactWYQ`) extremely confusing:

```julia
julia> qrf.Q
3×3 LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}:
 -0.945825 0.245839 0.212082
 -0.323885 -0.668782 -0.669201
 -0.0226788 -0.701637 0.712174

julia> size(qrf.Q)
(3, 3)

julia> qrf.Q * zeros(2, 2) # fine, no problem at all
3×2 Array{Float64,2}:
 0.0 0.0
 0.0 0.0
 0.0 0.0

julia> copyto!(Q, qrf.Q) # oh no, can't do that
ERROR: BoundsError: attempt to access 3×2 Array{Float64,2} at index [Base.OneTo(3), Base.OneTo(3)]
Stacktrace:
 [1] throw_boundserror(::Array{Float64,2}, ::Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}) at ./abstractarray.jl:467
 [2] checkbounds at ./abstractarray.jl:432 [inlined]
 [3] copyto!(::Array{Float64,2}, ::LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}) at ./multidimensional.jl:1247
 [4] top-level scope

```

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [March 31, 2018, 7:05am UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/3 "2018-03-31T07:05:34Z")

</div>

Ya it’s confusing. I am actually interested in the least allocating way to extract `Q`. I think even `qrfact!(A)` allocates for `Q` and `R` so it only overwrites `A` with intermediate results? Anyways, since I am not interested in `R` and this is a hot part of my program, I will write my own QR-like function. Thanks for the help.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [March 31, 2018, 11:52am UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/4 "2018-03-31T11:52:33Z")

</div>

Note that for small matrices you could use `qr` from StaticArrays: [https://github.com/JuliaArrays/StaticArrays.jl/blob/d3c5662aa8a8a65ab391dbb3bc58441dac3a1461/src/qr.jl#L4-L10](https://github.com/JuliaArrays/StaticArrays.jl/blob/d3c5662aa8a8a65ab391dbb3bc58441dac3a1461/src/qr.jl#L4-L10). If `A` is larger, you could use LAPACK functions directly for BLAS float types:

- [https://docs.julialang.org/en/latest/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK.geqrt](https://docs.julialang.org/en/latest/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK.geqrt)!
- [https://docs.julialang.org/en/latest/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK.geqrf](https://docs.julialang.org/en/latest/stdlib/LinearAlgebra/#LinearAlgebra.LAPACK.geqrf)!

These allow you to pass in work buffers that can be reused, saving a bit on allocation, but you’ll still need to extract the `Q` matrix from the modified `A`.

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [March 31, 2018, 5:14pm UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/5 "2018-03-31T17:14:41Z")

</div>

You could file an issue at stdlib/LinearAlgebra for a R-less QR factorization and a method to obtain `Q` through the Householder transformations without allocating. The code for `ldiv!` reveals that for solving a system of equations Julia uses `view(lmul!(adjoint(A.Q), b)`.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [March 31, 2018, 10:17pm UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/6 "2018-03-31T22:17:03Z")

</div>

I don’t think it is possible to do an R-less QR factorization, since R pretty much comes for free. But extracting the Q matrix allocation-free should be possible. I did a proof of concept for a fully in-place QR factorization in a bit over a 100 lines, but the only problem is getting it to give the right result!

[https://github.com/JuliaLang/julia/issues/26673](https://github.com/JuliaLang/julia/issues/26673)

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [March 31, 2018, 11:27pm UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/7 "2018-03-31T23:27:45Z")

</div>

The LAPACK routines need some workspace which is allocated by the Julia wrappers. Is this too much allocation for you?

```julia
julia> A2=randn(1000,500); tmp=fill(0.0,size(A2,2));

julia> Q2=copy(A2);

julia> @time LAPACK.geqrf!(Q2,tmp);
  0.047948 seconds (7 allocations: 125.344 KiB)

julia> @time LAPACK.orgqr!(Q2,tmp);
  0.026197 seconds (6 allocations: 125.313 KiB)

julia> F=qrfact(A2);

julia> vecnorm(A2-Q2*F.R) # F[:R] in v0.6
5.634945506084882e-13

```

IIRC the size of the workspace depends on the preferred block size for level-3 BLAS, so it makes sense to conceal it from ordinary users.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [April 1, 2018, 12:26am UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/8 "2018-04-01T00:26:04Z")

</div>

> [@Ralph\_Smith](#):
>
> Is this too much allocation for you?

I am aiming for 0 allocations 🙂

Here is a proof of concept in about 160 lines of code with the API, [GitHub - mohamed82008/InplaceQR.jl: Provides a fully in-place, thin implementation of the Householder-based QR decompositon in Golub and Van Loan's book.](https://github.com/mohamed82008/InplaceQR.jl). A few allocations slipped through so it’s time for some optimization.

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [April 1, 2018, 2:32am UTC](https://discourse.julialang.org/t/extracting-q-from-qrfact-allocation-free/10083/9 "2018-04-01T02:32:54Z")

</div>

As mentioned in the issue and previously hinted by @Nosferican, `lmul!` does the trick.