# Is the "practical" complexity of linsolve direct solver O(n^2) ?

I recently timed the linsolve direct solver and I was kind of shocked to see that the solver seemed to be scaling quadratically even upto a 1000 dimensions. Specifically I ran the following code and it showed that the cubic coefficient of time taken $a_3 \approx \frac{1}{1000}a_2$. From my recollection the cubic coefficient is not supposed to be so small but more like $\frac{1}{10}a_2$. Why is the routine so fast? I am on a 5 year old macbook pro with only 2 core.

clear; close all; clc
pv = [ 10 20 30 40 50 60 70 80 90 100 200 300 400 500 600 700 800 900 1000];
tv = nan(length(pv), 1);
runs=3;
tv2 = nan(length(pv), 1);
for pi = 1:length(pv)
p = pv(pi);
v = randn(p, 1);
m = randn(p, p);
% Calculate time for mldivide.
tic; for i = 1:runs mldivide(m, v); end; tv(pi) = toc/runs * 1e9;
% Calculate time for linsolve.
tic; for i = 1:runs linsolve(m, v); end; tv2(pi) = toc/runs * 1e9;
end
plot(pv, tv, 'b');
hold on;
plot(pv, tv2, 'r');
ylabel('nano-seconds');
xlabel('dimension');


## 1 Answer

Your matrices are far too small to see the asymptoptic $O(n^3)$ running time behavior of the LU factorization used by linsolve. For very small matrices the overhead of computations surrounding the LU factorization will make it difficult to see the $O(n^3)$ growth. Furthermore, MATLAB will typically be making use of parallel routines for computing the LU factorization that only start to be efficient at the smaller matrix sizes you've experimented with. Try your demo again using a range of sizes from $p=500$ up to say $p=20000$ and see what happens.

• I had to go upto 40k on my university's computer to experience the cubic scaling. But the question still stands. Specifically I would greatly appreciate if someone could shed a light on why the algorithm only exhibits cubic runtime at 10k dimension andnot before Dec 17 '15 at 0:24
• The keyword in my answer is "asymptotic" These kinds of complexity results hold in the limit as the problem size $n$ goes to infinity. You shouldn't expect things to work in that way for small values of $n$. Dec 17 '15 at 0:28
• Is there a reference about the nuts and bolts of the "overhead surrounding LU" decomposition. the reference about guassian elimination on wikipedia gives a different relative factor. Dec 17 '15 at 0:29
• At the level of computer architecture, it's important understand that on contemporary hardware memory accesses are vastly slower than floating point operations. Modern processors attempt to adjust for this by using high speed "cache" memory to store data that is used frequently. For your problem with $p=500$, the entire matrix will fit into 2 megabytes of cache (each core of contemporary processor is likely to have this much cache) For matrices of this size and smaller nearly all of the time will be spent loading $O(n^{2})$ data into cache from memory. Dec 17 '15 at 0:36
• You really should take a course in modern computer architecture if you want to understand in more detail what's happening with parallel processing and cache in solving these systems of equations. The "overhead surrounding the LU factorization" includes things like loading the data from main memory and starting multiple processes to work on the various blocks of the matrix. Until the data has been brought in from memory and unless the matrix is large enough to make multithreading worth while, the multiple processor cores won't be working at full efficiency. Dec 17 '15 at 0:42