# Is LAPACK behind the cutting edge of dense linear algebra?

I have been digging into some numerical linear algebra lately, and reading in particular about how LAPACK solves symmetric eigenvalue problems. I noticed that the *stegr routine for computing eigenvalues and eigenvectors of a symmetric tridiagonal matrix (to which the more general cases are reduced) is $$\mathcal{O}(n^2)$$ in time, as stated in the references (interestingly this last URL has the suffix "#holygrail"). The LAPACK documentation (same as first link) states that this algorithm (the relatively robust representation algorithm) "is faster than all the other routines except in a few cases, and uses the least workspace."

However, I noticed on the Wikipedia article for tridiagonal matrices that there exist $$\mathcal{O}(n \log{n})$$ algorithms for computing the eigenvalues of a real symmetric tridiagonal matrix.

The reference given for the $$\mathcal{O}(n \log{n})$$ algorithm is from 2012. Based on this table, it looks like the last release of LAPACK was in 2000. Is it then the case that LAPACK is significantly supoptimal in its performance capabilities? If so, why do so many packages (like Julia's LinearAlgebra or Python's SciPy) use LAPACK under the hood?

(Note that I saw this post, which indicates that the $$\mathcal{O}(n \log{n})$$ algorithm for symmetric tridiagonal matrices does not help the overall asymptotic performance for dense matrices since the tridiagonalization step is more expensive. But I am still prompted by the above to wonder: are there known ways that LAPACK does not perform as well as it could in a significant way?)

• LAPACK 3.9.0 was released in November 2019, netlib.org/lapack/lapack-3.9.0.html. it is actively developed, and I doubt you will find a more state of the art dense linear algebra library available. – vibe Jun 12 '20 at 5:42
• @vibe Ah thanks for this clarification! I was misled by the lapack/releases page that I linked and apparently did not dig deep enough after seeing it. – Grayscale Jun 12 '20 at 6:15
• I think the spirit of the question still survives. I edited the question a little though. – Grayscale Jun 12 '20 at 6:28
• Another subtle issue is that in this context there are two targets for "performing well": speed and accuracy. If this new method does not help asymptotic performance, then the choice will be mainly driven by numerical accuracy. – Federico Poloni Jun 12 '20 at 6:39
• Additionally, you shouldn't look at the asymptotic performance. For very large problems, you would likely look at parallel implementations to run on supercomputers anyway. For moderate problem sizes, cache efficiency, vectorisation and even register efficiency matter way more. – Thijs Steel Jun 12 '20 at 8:01

When one says an algorithm is of order $$O(n)$$, that may mean that the complexity is given by: $$c + b*n$$. With every new element you add you increase in runtime (effectively). What mathematically minded people often forget is that these statements do not include how large the constants are. That of course carries over to $$O(n²)$$ and such. I can not answer your question in absolute terms of whether or not there might be special cases where other libraries are faster, but complexity orders are sadly only half of the story. When you design an algorithm you often have tradeoffs which allow you to push work from the linear term to the constant term, i.e. doing certain precomputations before doing element-wise computations. Optimizing for order will not always improve your runtime for a given problem. The break-even point may be arbitrarily large. The only proper way here is to measure.
• The cited paper says the break even point is about $n=800$ – vibe Jun 12 '20 at 18:32