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);
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;
plot(pv, tv, 'b');
hold on;
plot(pv, tv2, 'r');

a busy cat


1 Answer 1


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.

  • $\begingroup$ 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 $\endgroup$
    – Pushpendre
    Dec 17, 2015 at 0:24
  • 3
    $\begingroup$ 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$. $\endgroup$ Dec 17, 2015 at 0:28
  • $\begingroup$ 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. $\endgroup$
    – Pushpendre
    Dec 17, 2015 at 0:29
  • 1
    $\begingroup$ 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. $\endgroup$ Dec 17, 2015 at 0:36
  • 3
    $\begingroup$ 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. $\endgroup$ Dec 17, 2015 at 0:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.