4
$\begingroup$

I have a problem where I need to diagonalize a large number of small Hermitian matrices. Typically the matrices are between 4 and 64 in size (skewed towards the low end) and the number of matrices is on the order of millions. This is then repeated many times, so I have been trying to figure out the best way to go about this. Most of the performance information I've been able to gather seems to focus on the case for matrices that are significantly larger.

So the question is: Is there an algorithm (or implementation thereof) that is particularly suited for diagonalizing small (Hermitian) matrices?

|cite|improve this question
$\endgroup$
  • 1
    $\begingroup$ Do you need to diagonalize the matrices or solve linear equations? Are the systems related in any way? If not, are they independent such that you could parallelize the diagonalization? $\endgroup$ – Bill Barth Jun 24 '12 at 3:49
  • 1
    $\begingroup$ They need to be diagonalized. The systems stem from a matrix-valued function, of which I need to do some integrals involving it's eigenvalues and eigenvectors. They are independent and parallelizing the computation is something I do. I've been toying with using the GPU for this, so knowing the algorithm itself would be necessary (as optimized implementations do not exist) so I could port it (to CUDA and/or OpenCL). $\endgroup$ – radioactivist Jun 24 '12 at 17:33
8
$\begingroup$

Let me add a few remarks explaining why there is no better method for "small" matrices in the range described (order 4 to 64) than the usual approach: tridiagonalize the Hermitian matrix with a unitary similarity transform (a direct step), then iteratively apply QR-like operations to approach a diagonal limit.

If the eigenvalues (diagonalized matrices) alone are needed, the reduction to tridiagonal form can be carried out with Householder transformations using $2n^3/3$ (complex) multiplications, or $4n^3/3$ multiplications if the eigenvectors (unitary similarity transforms) are also needed. This reduction step's complexity dominates that of the iterative steps, which is $O(n \log^3 n)$ "even in the case of high precision" (Reif,1999) if only eigenvalues are needed (increase by a factor of $n$ if eigenvectors too are needed). So in any case the usual combined complexity is $O(n^3)$ with a modest leading coefficient.

Finding eigenvalues is of course theoretically equivalent to getting roots of a characteristic polynomial, and their effective equivalence is outlined (with links to references) in the early portion of this previous Answer. The rest of that Answer explains why direct methods for orders 6 and above are impractical compared with the tried-and-true iterative approach(es).

For a 2x2 matrix it certainly makes sense to use the stable form of the quadratic formula to get the (real) eigenvalues directly. Probably a case can be made for using Cardano's formula for roots of a cubic to solve eigenvalues of a 3x3 matrix, or perhaps the formulas of Viète if we are to avoid complex arithmetic.

But there the "garden path" of direct methods comes to a sudden narrowing! Although the general quartic can be solved by radicals (Ferrari's method) and the general quintic polynomial by Bring radicals, these methods require solving "auxiliary polynomials" of lower degrees. Such cumbersomeness is compounded by complex arithmetic and monitoring for loss of significance in a floating point implementation.

The situation for order 6 and above is even more inauspicious because in place of the single argument Bring radical, exotic functions of at least two arguments are needed. Moreover the known methods that could be coded for variable degrees $n$ involve an element of trial-and-error to locate all roots.

In contrast the usual approach adapts to variable orders $n$ with no apparent difficulty and takes efficient advantage of the smallness of $n$ in both time and space.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for the great answer! The case you've made here against direct methods is quite convincing, so I guess any improvements would have to come from optimizing a specific implementation of a more usual method. $\endgroup$ – radioactivist Jun 24 '12 at 18:52
  • 1
    $\begingroup$ @radioactivist: It's been an active area of research, and as Arnold Neumaier says, the LAPACK folks stay on top of both the algorithm and platform developments. It's been a tricky enough computation not to want to implement from scratch (other than for educational purposes) for a long time. This older mathSE thread discusses the complete eigenvalue/eigenvector problem for slightly larger but sparse matrices, again mainly recommending LAPACK/BLAS. $\endgroup$ – hardmath Jun 24 '12 at 23:33
7
$\begingroup$

The most efficient package for platform adapted dense linear algebra (including eigensolvers) is the LAPACK library. http://en.wikipedia.org/wiki/LAPACK Every improvement in the state of the art is fairly soon reflected in this library, and a lot of effort is spent in tuning the kernels to optimal performance on all platforms.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ This is actually what I've been using for the most part (the Intel MKL specifically). However, given some projects such as Eigen claim better performance for this case (see here for example, which doesn't even use their fixed size code) I was wondering if there was some room for improvement. $\endgroup$ – radioactivist Jun 24 '12 at 19:32
  • $\begingroup$ @radioactivist: I think the eigen3 benchmarks use dynamically sized matrices for all the libraries, about which they wrote: "The reason why the values are typically low for small sizes, is that in this benchmark we deal with dynamic-size matrices which are relatively inefficient for small sizes." So their slight edge at those low sizes may have more to do with dynamic vs. fixed matrix manipulations, rather than the computational linear algebra. $\endgroup$ – hardmath Jun 25 '12 at 0:08
  • $\begingroup$ @hardmath I would think though that their fixed size code (size known at compile time) would be even faster at those low sizes (given this is what it is designed for) $\endgroup$ – radioactivist Jun 25 '12 at 0:19
  • $\begingroup$ @radioactivist: True, I expect it would be, but so would the code of the other libraries included in the benchmark. And the change might erase any slight advantage held by eigen3 at the small orders. $\endgroup$ – hardmath Jun 25 '12 at 1:17
  • $\begingroup$ @hardmath I don't think the other libraries have such optimizations for fixed size objects (as they can't do the optimizations at compile time) $\endgroup$ – radioactivist Jun 25 '12 at 1:53

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.