# Optimizing a matrix computation algorithm

In the frame of the Kabsch algorithm to compute a spatial rotation matrix, I am implementing the following workflow:

• start from a given "cross-covariance" $3\times3$ matrix $A$ and compute its Singular Value Decomposition $A=U\Sigma V^T$; for this, compute the symmetric matrix $B=A^TA=V\Sigma^2V^T$ and its Eigendecomposition $B=V\Lambda V^T$; $\Lambda$ will give $\Sigma$;

• the Eigenvalues are obtained by solving the characteristic equation analytically, by means of trigonometric functions (as $B$ is symmetric, we are in the Casus Irreductibilis);

• from an Eigenvalue $\lambda$, an Eigenvector is obtained as three minors of the singular matrix $B-\lambda I$;

• when $V$ is known, $U$ can be computed as $AV\Sigma^{-1}$;

• the final answer is the matrix $R=UV^T=AV\Sigma^{-1}V^T$.

As the procedure will be invoked many times, I was wondering if any clever simplification was possible along the path from $A$ to $R$, for maximum efficiency. (I don't think that switching to an iterative approach can yield a faster solution, but I may be wrong.) I only address the $3\times3$ case, and I momentarily don't care about numerical stability nor degenerate cases.

Update:

As a variant, after $\Lambda=\Sigma^2$ and $V$ are found by diagonalizing $A^TA$, the Eigenvalues can be reused to diagonalize $AA^T=U\Sigma^2U^T$ and obtain $U$.

• The geometric least squares (GLS) plane of a set of points can also be computed using SVD (the normal of the least squares plane of the centered point set is the singular vector for the smallest singular value). But I have found it is much faster (approximately 30 times faster) to compute the GLS plane using Newton's method to solve a system of two quartics for the two unknowns defining the normal and an initial guess of the best ordinary least squares plane normal. So it might also be faster to use Newton's method to find this rotation matrix rather than the SVD. – J. Heller Dec 7 '16 at 3:36
• @J.Heller: I can believe a "30 times faster" against a general purpose/arbitrary dimension SVD solver, but not against an ad-hoc $3\times3$ algorithm. – Yves Daoust Dec 7 '16 at 10:57
• This timing comparison used Atlas optimized LAPACK for SVD, which I have found to perform very similarly to the Intel MKL. Finding the roots of a cubic equation to compute the SVD would definitely be much faster than this, but is it robust? – J. Heller Dec 7 '16 at 16:52
• @J.Heller I momentarily don't care about numerical stability nor degenerate cases. – Yves Daoust Dec 7 '16 at 18:05
• A timing comparison of symmetric 2x2 eigenvalue and eigenvector computation by solving a quadratic equation versus optimized LAPACK shows that the quadratic equation based method is 7-8 times faster than LAPACK. Do you have any idea how much faster symmetric 3x3 eigenvalue/eigenvector computation by solving a cubic is than using optimized LAPACK or Intel MKL? – J. Heller Dec 16 '16 at 17:24

The following paper might be of interest:

Computing the Singular Value Decomposition of 3×3 matrices with
minimal branching and elementary floating point operations
University of Wisconsin Department of Computer Science
Technical Report #1690


http://pages.cs.wisc.edu/~sifakis/project_pages/svd.html

http://pages.cs.wisc.edu/~sifakis/papers/SVD_TR1690.pdf

They compare their method to the method you use in Figure 4, getting modest time savings (84ns vs. 112ns on single core, 4.8ns vs. 7.0ns on 12 core).

The algorithm does look somewhat complicated. Someone has written an implementation that works on GPUs here: https://github.com/ericjang/svd3

• Also relevant: Müller et al., "A Robust Method to Extract the Rotational Part of Deformations", Proc. Motion in Games 2016 – Rahul Dec 7 '16 at 0:14
• Yep, the Sifakis paper is quite interesting. They seem to prefer iterative methods such as Jacobi. Anyway they compare SIMD optimized approaches to non optimized ones, which biases the comparisons. I was unaware of the paper by Müller et al. Also interesting (and using an original iterative approach). – Yves Daoust Dec 7 '16 at 0:41