Can one outperform Cramer's rule for the inversion of a 3 by 3 matrix

Question

I know that for general matrices, Cramer's rule is far from ideal for the numerical computation of the matrix inverse. However, can it be outperformed in the case of a $3 \times 3$ matrix? One particular advantage is that the algorithm is non-branching. For example, the transpose of the inverse can be implemented as follows, assuming the matrix $\bf{a}$ is invertible and the matrices $\bf{a}$ and $\bf{b}$ are stored as contiguous arrays.

void inverse3(double* a, double* b) {
    // cofactors
    b[0] = a[4]*a[8] - a[5]*a[7];
    b[1] = a[5]*a[6] - a[3]*a[8];
    b[2] = a[3]*a[7] - a[4]*a[6];
    b[3] = a[7]*a[2] - a[8]*a[1];
    b[4] = a[8]*a[0] - a[6]*a[2];
    b[5] = a[6]*a[1] - a[7]*a[0];
    b[6] = a[1]*a[5] - a[2]*a[4];
    b[7] = a[2]*a[3] - a[0]*a[5];
    b[8] = a[0]*a[4] - a[1]*a[3];
    // determinant
    double det = b[0]*a[0] + b[1]*a[1] + b[2]*a[2];
    // inverse
    b[0] /= det;
    b[1] /= det;
    b[2] /= det;
    b[3] /= det;
    b[4] /= det;
    b[5] /= det;
    b[6] /= det;
    b[7] /= det;
    b[8] /= det;
}

Can this be improved in terms of efficiency or numerical robustness?

Answer 1

I don't know if another factorization would be faster overall, but you can (almost) certainly speed the above code up by inverting the determinant once and multiplying by the inverse 9 times rather than performing the division over and over. It's possible the compiler will figure this out for you, but not guaranteed (look at the assembly).

Edited to add: For the record, a couple of attempts with the Intel 13 compiler give different results depending on what it knows about the divisor. If it comes from command line input, then the compiler dutifully divides each one (vectorized on my Sandy Bridge machine).

Answer 2

In terms of numerical robustness, yes, for sure. For instance, det could overflow in your code; a method based on a Givens QR factorization won't instead. I guess also that there could be significant cancellation in your cofactor computations.

I am afraid that you can't have perfect stability and minimal number of operations at the same time; you'll have to choose in which of the two directions you wish to improve.