A QR complexity question

It is a common understanding that the complexity of QR is O(n^3). But if I do the following simple tests,

>> tic; for k =1:100; a = randn(100); [q r] = qr(a); end; toc

Elapsed time is 0.082818 seconds.

>> tic; for k =1:100; a = randn(1000); [q r] = qr(a); end; toc

Elapsed time is 5.149743 seconds.

why the elapsed time is only increased by about 63 (instead of 1000) in my system ?

Answer 1

[Bruno Luong]

The O(n^3) is number of flops, which is not proportional to tic/toc.

• You time also RANDN
• Calling QR has overhead that is significant to aithmetic operations
• CPU caching is less effective for small size
• Multithreading/parallel is less effective for small size

[David Tung]

Multithreading probably is the major factor.

I want to use multiple tries to make sure the result is more stable, so I include randn in the loop. Thanks for the response.

Answer 2

[John D'Errico]

(Good points made by Bruno. This answer tries to dive a little more deeply than Bruno spent the time to do, but I really don't care if you just accept Bruno's answer, as it really did answer the heart of the question.)

To get a better answer as to how much time QR really takes requires much deeper investigation. And sometimes, there are incredibly strange interactions between various parts of your particular computer. There are several important things you need to do.

1. Look at a larger set of array sizes
2. Do not put the randn time into your time computation
3. Use timeit, as a far better measure of the time taken than tic/toc
4. Finally, plot AND model the results

I would do something like

N = [100:50:1000, 1100:100:2000, 2250:250:4000, 4500, 5000];
T = zeros(size(N));
for i = 1:numel(N)
    A = randn(N(i));
    T(i) = timeit(@() qr(A))
end

Now we would expect the resulting times to roughly follow a power model, of the form

T(N) = a + N^b

where b may be on the order of 3.

loglog(N',T','o')

Note that if this is a perfect relationship, we would expect to see a straight line on that plot.

But what happens is for smaller matrix sizes, you have other factors in there. A huge part of it is just function call overhead, thus startup costs, error checking, etc. Those startup costs often take a virtually constant time. But that confuses things.

polyfit(log(N),log(T),1)
ans =
       2.4965      -21.774

This suggests, using ALL of those times, that a reasonable model might have the QR time as O(N^2.5), at least based on my data. However, if we look only at the last 4 points from that curve, because there was a jog in the curve just below that...

k = 36:39;

fit(log(N(k))',log(T(k))','poly1')
ans = 
     Linear model Poly1:
     ans(x) = p1*x + p2
     Coefficients (with 95% confidence bounds):
       p1 =       3.203  (2.823, 3.583)
       p2 =      -27.55  (-30.73, -24.37)

We see a time behavior as possibly closer to O(N^3.2). Since that uses only 4 data points, it might be a little suspect. And the exponent bounds clearly contain your hoped for value of 3.

Next, you might want to recognize those jogs in the curve. They are pretty significant. They may be subtle things induced by my own computer, how it uses memory, or interactions from the blas, or perhaps all sort of things I have not even considered. Some of those bumps might be where MATLAB suddenly decided to multi-thread the problem. My computer has 8 real cores.

[David Tung]

Amazing answer! Thanks.