### Author:  [@qrebjock](https://discourse.julialang.org/u/qrebjock)
#### Post date: [March 20, 2021, 3:14pm UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/1 "2021-03-20T15:14:43Z")
Hello,
It seems like the function `qr` doesn’t take keywords argument to return an economy-sized factorization.
In this thread [qr(A; full=false) no longer available? · Issue #27397 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/27397) they say that the internal representation is efficient so it shouldn’t be a problem. But when it comes to do matrix multiplications with the Q factor I get weird performances.
Here is an example:
```julia
function fun1()
A = randn(10000, 500)
B = randn(10000, 500)
Q, R = qr(A)
return Q * B
end
function fun2()
A = randn(10000, 500)
B = randn(500, 500)
Q, R = qr(A)
return Q * B
end
```
The function `fun2` should use the thin QR when evaluating the product `Q * B`.
But when we time the functions
```julia
using BenchmarkTools
@btime fun1();
@btime fun2();
```
we find that `fun2` is as slow as `fun1`, or even slower because there are more allocations.
So my question would be: how to take advantage of the thin QR decomposition to get better performances?
Side note: If I try to do the multiplication with `Q[:, 1:500] * B` it is even slower (actually it takes forever – not sure what’s happening here).
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [March 20, 2021, 3:33pm UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/2 "2021-03-20T15:33:16Z")
What is the `typeof(Q) `? Maybe you need to do `Matrix(Q) ` to get an ordinary matrix which is fast to multiply by. The Q returned by qr might be lazy.
Possibly related
[https://github.com/JuliaLang/julia/issues/38972](https://github.com/JuliaLang/julia/issues/38972)
[https://github.com/JuliaLang/julia/issues/37102](https://github.com/JuliaLang/julia/issues/37102)
---
### Author:  [@ettersi](https://discourse.julialang.org/u/ettersi)
#### Post date: [March 21, 2021, 2:50am UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/3 "2021-03-21T02:50:06Z")
`Q*Bthin` internally zero-extends the `n x n` matrix `Bthin` to an `m x n` matrix `Bfull` and then computes `Q*Bfull`, see
>
That’s why you do not get any speedup in your `fun2()`.
According to the documentation for `qr()`, you can get the thin `Q` factor as a dense matrix through `Matrix(Q)`, so `Matrix(Q)*Bthin` avoids the zero-padding, but it seems that the dense matrix algebra is less efficient than working with the compressed Householder representation:
```julia
julia> A = rand(10000,500)
Bfull = rand(10000,500)
Bthin = rand(500,500)
Q,R = qr(A)
@time Q*Bfull
@time Matrix(Q)*Bthin
;
0.184003 seconds (4 allocations: 40.054 MiB)
0.300148 seconds (8 allocations: 116.349 MiB, 5.89% gc time)
```
---
### Author:  [@qrebjock](https://discourse.julialang.org/u/qrebjock)
#### Post date: [March 21, 2021, 8:59pm UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/4 "2021-03-21T20:59:19Z")
Thank you for your replies.
The type of `Q` is `LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}`. And indeed, `Matrix(Q)` returns a `10000 x 500` matrix. I understand better the internal representations now.
I have the same results as you, ettersi. When I perform the operation `Matrix(Q) * Bthin` it is as slow, or even slower than `Q * Bfull`.
Do you think it is something that we can hope to improve? That is, is it an issue with the way Julia builds `Matrix(Q)`? Or there is some inherent limitation of the `QR` decomposition that makes the thin decomposition not much faster than the full one?
---
### Author:  [@ettersi](https://discourse.julialang.org/u/ettersi)
#### Post date: [March 22, 2021, 1:17am UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/5 "2021-03-22T01:17:04Z")
On second thought, it is not too hard to see why `Matrix(Q) * Bthin` is slower than `Q * Bthin`: the `QRCompactWYQ` data type represents `Q` as a sequence of Householder reflectors,
Q = \prod\_{k = 1}^n (I - 2 \, u\_k u\_k^T) .
Therefore:
- `Q * Bthin = Q * [Bthin; zeros(m-n,n)]` requires `O(m*n)` operations.
- `Matrix(Q) * Bthin` requires `O(m*n^2)` operations.
This probably also explains why Julia (or LAPACK for that matter) does not provide an “optimised” code for `Q*Bthin`. The runtime of this operation would still be `O(m*n)`, i.e. it would have the same asymptotic complexity as `Q*Bpadded` and “only” a slightly better prefactor. Specifically, I believe `Q*Bpadded` requires roughly `5*m*n` additions and multiplications while an optimal implementation of `Q*Bthin` would require `(2m + 3n)*n` additions and multiplications.
---
### Author:  [@qrebjock](https://discourse.julialang.org/u/qrebjock)
#### Post date: [March 22, 2021, 9:02am UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/6 "2021-03-22T09:02:26Z")
I see, thank you that’s very helpful!
It seems like that computing the product `Q * Bthin` is already mostly optimized.
I am just wondering: is the zero padding allocating memory? Or it is something performed lazily and internally Julia knows that there is just a bunch of zeros?
---
### Author:  [@ettersi](https://discourse.julialang.org/u/ettersi)
#### Post date: [March 22, 2021, 3:41pm UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/7 "2021-03-22T15:41:53Z")
No, the zero padding allocates memory.
---
### Author:  [@qrebjock](https://discourse.julialang.org/u/qrebjock)
#### Post date: [March 22, 2021, 4:58pm UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/8 "2021-03-22T16:58:14Z")
Ok, I think this explains why it doesn’t scale linearly in `m` (considering `A` is `m x n`). It’s slightly more than linear.
When I am trying the same code on MATLAB it does scale linearly if I use the thin QR decomposition.
Is there any chance to rewrite things in a way that would be linear? Maybe with low-level LAPACK functions? I am looking into it but it’s not so clear.
---
### Author:  [@ettersi](https://discourse.julialang.org/u/ettersi)
#### Post date: [March 25, 2021, 12:52am UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/9 "2021-03-25T00:52:46Z")
I suspect that Julia failing to be linear in `m` is actually a sign of how optimised the code is. Specifically, I suspect that the nonlinearity is due to slower memory accesses if you operate on a larger amount of data. If so, the only way to achieve linear scaling would be to slow down the fast processing of small amounts of data which of course would make no sense.
This is just a guess. We might be able to come to more meaningful conclusions if you share the hard data which lead to your conclusion that Matlab scales linearly but Julia doesn’t.
---
### Author:  [@qrebjock](https://discourse.julialang.org/u/qrebjock)
#### Post date: [April 5, 2021, 4:46pm UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/10 "2021-04-05T16:46:39Z")
Ah I see what you mean. Yes it wouldn’t make sense. Actually I benchmarked a little more rigorously and here are the results I get.
In Julia:
```julia
using LinearAlgebra
using Plots
function fun(A, B)
Q = qr(A).Q
C = Q * B
end
function benchmark(n, sizes, num_repeats)
num_sizes = length(sizes)
times = zeros(num_sizes)
for i = 1:num_sizes
s = 0
for repeat = 1:num_repeats
A = randn(sizes[i], n)
B = randn(n, n)
s += (@timed fun(A, B)).time
end
times[i] = s / num_repeats
end
return times
end
n = 500
num_repeats = 5
num_sizes = 20
sizes = Int.(floor.(exp10.(range(3, stop=4, length=num_sizes))))
times = benchmark(n, sizes, num_repeats)
plot(sizes, times, xaxis=:log, yaxis=:log, label="")
plot!(sizes, sizes ./ 1e5, xaxis=:log, yaxis=:log, label="")
```
The plot I get:
In MATLAB:
```matlab
n = 500;
num_repeats = 5;
num_sizes = 20;
sizes = floor(logspace(3, 4, num_sizes));
times = zeros(num_sizes, 1);
for i = 1:num_sizes
disp(sizes(i))
s = 0;
for repeat = 1:num_repeats
A = randn(sizes(i), n);
B = randn(n, n);
tic
Q = qr(A, 0);
C = Q * B;
s = s + toc;
end
times(i) = s / num_repeats;
end
hold on
set(gca, 'XScale', 'log', 'YScale', 'log')
loglog(sizes, times)
loglog(sizes, sizes / 1e5)
```
The plot I get:
As you see both scale more or less linearly. But MATLAB is much faster.
There are maybe a few things to change in my code for the benchmark to be more fair. Let me know if it’s the case.
---
### Author:  [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)
#### Post date: [April 6, 2021, 6:34am UTC](https://discourse.julialang.org/t/economy-sized-qr-decomposition/57608/11 "2021-04-06T06:34:59Z")
What do you get from
```julia
julia> LinearAlgebra.BLAS.vendor()
:openblas64
```
?
Matlab uses MKL as BLAS, which is usually much faster than openblas. You can use MKL by installing MKL.jl