MATLAB Answers

A QR complexity question

27 views (last 30 days)
David Tung
David Tung on 26 Jan 2021
Commented: David Tung on 27 Jan 2021
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 ?

  0 Comments

Sign in to comment.

Accepted Answer

Bruno Luong
Bruno Luong on 26 Jan 2021
Edited: Bruno Luong on 26 Jan 2021
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

  1 Comment

David Tung
David Tung on 27 Jan 2021
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.

Sign in to comment.

More Answers (1)

John D'Errico
John D'Errico on 26 Jan 2021
Edited: John D'Errico on 26 Jan 2021
(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.

  1 Comment

David Tung
David Tung on 27 Jan 2021
Amazing answer! Thanks.

Sign in to comment.

Tags

Products


Release

R2018a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!