6
$\begingroup$

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?)

$\endgroup$
5
  • 7
    $\begingroup$ 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. $\endgroup$ – vibe Jun 12 '20 at 5:42
  • $\begingroup$ @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. $\endgroup$ – Grayscale Jun 12 '20 at 6:15
  • $\begingroup$ I think the spirit of the question still survives. I edited the question a little though. $\endgroup$ – Grayscale Jun 12 '20 at 6:28
  • 5
    $\begingroup$ 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. $\endgroup$ – Federico Poloni Jun 12 '20 at 6:39
  • 6
    $\begingroup$ 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. $\endgroup$ – Thijs Steel Jun 12 '20 at 8:01
8
$\begingroup$

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.

$\endgroup$
1
  • 4
    $\begingroup$ The cited paper says the break even point is about $n=800$ $\endgroup$ – vibe Jun 12 '20 at 18:32
6
$\begingroup$

LAPACK has been on the cutting edge for just about three decades, and probably still is for its niche. However, given given recent developments in libraries for the simpler BLAS-type matrix operations that LAPACK traditionally builds upon, it is perhaps conceivable that we could see the emergence of serious competitors to the traditional FORTRAN-based LAPACK in coming years.

In other words, libraries such as BLIS (C), MKL (C/C++) and Octavian.jl (Julia), which appear to beat traditional BLAS for the first time in decades, suggest that there is indeed performance still on the table even for such fundamental operations.

For comparison (via the Octavian.jl repo): BLAS-type library benchmarks

$\endgroup$
2
  • $\begingroup$ There are many libraries that are built over LAPACK. Would they be able to transition if they wanted to? $\endgroup$ – CuteCompute Jul 2 at 15:43
  • $\begingroup$ MKL seems to provides a complete LAPACK interface so that it can just be dropped in (so for instance you can use PETSc with MKL's BLAS and LAPACK, e.g. software.intel.com/content/www/us/en/develop/articles/…) -- the only catch being that MKL is apparently freeware but not necessarily open source. For future libraries it's hard to say! $\endgroup$ – cbk Jul 2 at 16:00
3
$\begingroup$

The reference LAPACK library prioritizes portability over performance, so it's missing many optimizations, such as vectorization and threading. (Additionally, it relies on BLAS for much of it's performance, so using the reference BLAS will significantly limit performance.)

However, most LAPACK users don't use the reference library, but various optimized BLAS+LAPACK libraries that provide the same API, but optimize various routines. Examples include Intel's MKL, IBM's ESSL, OpenBLAS, Cray's libsci, ARM's performance libraries, and CUDA's cuBLAS. This shared API allows a code to be written once while being able to use the best LA implementation available on a particular system.

Julia by default uses OpenBLAS, and SciPy appears to let the user control their choice of implementation.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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