Economy-sized QR decomposition

[qrebjock]

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:

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:

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.