8
$\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?)

Improve this question
$\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, 2020 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, 2020 at 6:15
  • $\begingroup$ I think the spirit of the question still survives. I edited the question a little though. $\endgroup$
    – Grayscale
    Jun 12, 2020 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, 2020 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, 2020 at 8:01

3 Answers 3

Reset to default
9
$\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.

Improve this answer
$\endgroup$
1
  • 4
    $\begingroup$ The cited paper says the break even point is about $n=800$ $\endgroup$
    – vibe
    Jun 12, 2020 at 18:32
8
$\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

Improve this answer
$\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, 2021 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, 2021 at 16:00
4
$\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.

Improve this answer
$\endgroup$

Your Answer

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

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