Menu
I am working on a project where I'm basically preforming PCA millions of times on sets of 20-100 points. Currently, we are using some legacy code that is using GNU's GSL linear algebra pack to do SVD on covariance matrix. This works, but is very slow.
I was wondering if there are any simple methods to do eigen decompositions on a 3x3 symmetric matrix, so that I can just put it on the GPU and let it run in parallel.
Since the matrices themselves are so small, I wasn't sure what kind of algorithm to use, because it seems like they were designed for large matrices or data sets. There's also the choice of doing a straight SVD on the data set, but I'm not sure what would be the best option.
I have to admit, I'm not stellar at Linear Algebra, especially when considering algorithm advantages. Any help would be greatly appreciated.
(I'm working in C++ right now)
397
For every real symmetric matrix A there exists an orthogonal matrix Q and a diagonal matrix dM such that A = (QT dM Q). This decomposition is called a spectral decomposition of A since Q consists of the eigenvectors of A and the diagonal elements of dM are corresponding eigenvalues.
Abstract. Real symmetric 3 × 3 matrices have 6 independent entries (3 diagonal elements and 3 off-diagonal elements) and they have 3 real eigenvalues (λ0,λ1,λ2). If 2 of these 3 eigenvalues are the same, i. e., λ0 = λ1 = λ2, then the number of independent entries of this matrix is reduced from 6 to 4.
Eigendecomposition is used to decompose a matrix into eigenvectors and eigenvalues which are eventually applied in methods used in machine learning, such as in the Principal Component Analysis method or PCA.
is real. That is, the eigenvalues of a symmetric matrix are always real. is orthogonal.
Using the characteristic polynomial works, but it tends to be somewhat numerically unstable (or at the very least inaccurate).
A standard algorithm to compute eigensystems for symmetric matrices is the QR method. For 3x3 matrices, a very slick implementation is possible by building the orthogonal transform out of rotations and representing them as a Quaternion. A (quite short!) implementation of this idea in C++, assuming you have a 3x3 matrix and a Quaternion class, can be found here. The algorithm should be fairly suitable for GPU implementation because it's iterative (and thus self-correcting), can make reasonably good use of fast low-dimensional vector math primitives when they're available and uses a fairly small number of vector registers (so it allows for more active threads).
59
Most methods are efficient for bigger matrices. For small ones the analytical method ist the quickest and simplest, but is in some cases inaccurate.
Joachim Kopp developed a optimized "hybrid" method for a 3x3 symmetric matrix, which relays on the analytical mathod, but falls back to QL algorithm.
An other solution for 3x3 symmetric matrices can be found here (symmetric tridiagonal QL algorithm).
27
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With
