6

Four by four real (well, floating point...) matrices are used in computer graphics to represent projections. Sometimes we need to compute their inverses. How many multiplications are required?

  • Valve's Source SDK implements the "pen&paper" algorithm: start with a 4x8 matrix whose left hand side is $A$ (the matrix to invert) and whose right hand side is $I_4$ (the 4x4 identity matrix); apply Gauss moves until the left hand side is $I_4$ and the right hand will be $A^{-1}$. This algorithm requires (by inspecting the code) $4 \cdot 8 + 4 \cdot 4 \cdot 8 = 160$ multiplications.
  • Solving $AX = I$ applying Cramer's rule for each element and computing each 3x3 determinant with Sarrus's rule improves on that: we need $6 \cdot 16$ multiplications for the cofactors, then $4$ more for the determinant of the 4x4 matrix, plus $16$ multiplications for the inverse of that, for a total of $6 \cdot 16 + 4 + 16 = 116$ multiplications.
  • Cleverly precalculating products of two elements as in this answer requires (again inspecting the code) $2\cdot 12 + 6 + 4 \cdot 16 = 94$ multiplications.

Can we do better? Can we prove we can't?

Community
  • 1
Jacopo Notarstefano
  • 2,297
  • In the third method (80), you can shave off four multiplications by multiplying each of c0..c5, s0..s5 by invdet instead of multiplying each of 16 entries by invdet. – vadim123 Aug 22 '13 at 20:45
  • 2
    The answer you link to is specialized to matrices with bottom row $0\ 0\ 0\ 1$, which allows simplification. The first two bullets invert a general $4 \times 4$ matrix. – Ross Millikan Aug 22 '13 at 20:47
  • @vadim123 I agree with you. This brings the total to 90 (80 was wrong, I was forgetting 6 multiplications required for invdet and also failing basic arithmetic...) – Jacopo Notarstefano Aug 22 '13 at 21:07
  • 1
    Keep in mind that more goes into the actual running time than simply the number of inverses. Additions, subtractions, and inversions add up. Register pressure can cause problems. Latency might be relevant. Some implementations might vectorize. And so forth. –  Aug 22 '13 at 21:07
  • @RossMillikan I don't agree with you. That code is precalculating some 2x2 determinants in such a way that the expansion of each cofactor requires 3 multiplications of matrix elements by those precalculated determinants. Maybe I'm misunderstanding something. Where do you claim the code is exploiting the structure of the matrix? – Jacopo Notarstefano Aug 22 '13 at 21:11
  • I just saw the comments that you could view the matrix as a block matrix and invert a $3 \times 3$ instead of a full $4 \times 4$, but didn't chase it. – Ross Millikan Aug 22 '13 at 21:17
  • The $4 \times 4$ matrices used in graphics are very special. The top-left $3 \times 3$ matrix is a rotation, and I would suspect that this makes things much easier. – bubba Aug 22 '13 at 23:52
  • Also, just out of curiosity, why do you want to minimize multiplications? I can see that you might want to minimize computation time, but that's a differemt question. – bubba Aug 22 '13 at 23:53
  • @bubba I was browsing Valve's code and thought the "pen&paper" algorithm was suboptimal. Number of multiplications was the first metric I thought of, but yes, it's clearly not the whole story. – Jacopo Notarstefano Aug 23 '13 at 19:37
  • So, can we assume that the top left 3x3 submtarix is a rotation? – bubba Aug 24 '13 at 10:22
  • @bubba No, I'm interested in the general problem. Most graphics library implement a general inverse of a 4x4 matrix. I think this means that sometimes the top left 3x3 submatrix is not a rotation. – Jacopo Notarstefano Aug 25 '13 at 16:54
  • You're right. Graphics libraries often implement a general $4 \times 4$ inverse. But they don't need to. Or, at least, they could get a significant speed-up by writing special code for the case where the top left 3x3 submatrix is a rotation. In viewing transforms, this special case is very common. – bubba Aug 26 '13 at 00:33

2 Answers2

5

We can invert a 2x2 matrix with 6 multiplies and one inversion. By Strassen's algorithm, we can multiply 2x2 matrices with 7 multiplies.

Partition your 4x4 matrix as

$$\left( \begin{matrix} A & B \\ C & D \end{matrix} \right) $$

We will reduce this to the identity.

Compute the inverse of $A$.

Perform Gaussian elimination

$$\left( \begin{matrix} I & A^{-1} B \\ 0 & D - CA^{-1}B \end{matrix} \right) $$

Let $E = D - C A^{-1} B$. Now compute $E^{-1}$. We can now finish Gaussian elimination.

Repeating these operations on an identity matrix tells us that the inverse is

$$ \left( \begin{matrix} A^{-1} + A^{-1} B E^{-1} C A^{-1} & -A^{-1} B E^{-1} \\ -E^{-1} C A^{-1} & E^{-1} \end{matrix} \right) $$

The calculation of this inverse requires two matrix inversions (12 multiplies and 2 real inversions), and six 2x2 multiplies:

  • $C A^{-1}$
  • $(C A^{-1}) B$
  • $E^{-1} (C A^{-1})$
  • $A^{-1} B$
  • $(A^{-1} B) E^{-1}$
  • $(A^{-1} B)(E^{-1} C A^{-1})$

for 54 multiplies and 2 real inversions in all.

Of course, using Strassen's algorithm for 2x2 matrices is a terrible idea. I don't know if the above organization of the calculation is reasonable or not even if you don't use Strasses algorithm; I suspect it is unlikely.

If $A$ turns out to be singular, you have to partition the matrix differently to do the calculation.

  • I'd rather say that we can invert a $2\times 2$ matrix with $2$ multiplications and $4$ divisions, so your Gauss-tensor-elimination algorithm only requires $46$ multiplications and $8$ divisions. – Jack D'Aurizio Feb 03 '14 at 22:01
  • @Hurkyl One question: If the both the original matrix and A are invertible, Can we guarantee that the matrix E is always invertible? I mean, do we only have to worry about A not being invertible? Thanks. – D1X May 10 '17 at 11:31
  • 1
    @D1X: Yes. It shouldn't be hard to show that a block triangular matrix whose diagonal has square blocks is invertible iff all of the diagonal blocks are invertible. Also, its determinant is the product of the determinant of the diagonal blocks. If nothing else, you can just consider the entries of the product of the matrix with its unknown inverse: $$\begin{pmatrix} P & Q \ 0 & R \end{pmatrix} \begin{pmatrix} W & X \ Y & Z \end{pmatrix} = \begin{pmatrix} I & 0 \ 0 & I \end{pmatrix}$$ –  May 10 '17 at 11:47
1

In his paper Gaussian Elimination is not optimal, Volker Strassen stated the following:

enter image description here

A more recent treatment of matrix inversion is provided by Intel. However, they are aiming at speed rather than at the minimum number of multiplications.

Axel Kemper
  • 4,943