0

I need to invert an Eigen matrix (9x9 in my particular case) as a part of code that I want to automatically differentiate using CppAD. For this to succeed the code executing the inversion can not contain any branching like for example if or switch statements. Unfortunately, the inverse function of Eigen contains branching with makes the algorithmic differentiation of CppAD fail.

Mathematically it should be possible to come up with a formulation that does not need branching for a fixed matrix size that is guaranteed to be invertible. Is that correct?

Do you know of any library that implements such an inverse without branching?

Manumerous
  • 455
  • 6
  • 21
  • 1
    I have used block matrix inversion and Eigen's fixed-size vectorized inversion routines (4 x 4 and smaller) to compute 10 x 10 inverses. You would need to know if a particular subblock on the diagonal (say the top left corner) is invertible at compile time. In my case, it was actually much faster and just as accurate as LU decomposition. https://en.wikipedia.org/wiki/Block_matrix#Block_matrix_inversion – Charlie S Nov 22 '21 at 22:17
  • Thanks that is a great input! I have a particular case where the bottom right matrix D is actually zero for which it makes sense that you can get a significant speed up. Unfortunately, my top left matrix A is already of dimensionality 7x7. But I guess I could just use the same block inverse technique to obtain the inverse of A. Thanks a lot! – Manumerous Nov 22 '21 at 22:23
  • 1
    Exactly -- you can mix and match block sizes to meet your specific criteria (e.g. 4 x 4 and 3 x 3 blocks on the diagonal for your 7 x 7). – Charlie S Nov 22 '21 at 22:56

2 Answers2

1

In general, it is possible to invert arbitrary large (invertible) matrices without branching, but this gets inefficient for bigger matrices. Eigen only does this for matrices up to size 4x4.

If you want the derivation of the inverse, just use the identity

deriv (inv(A)) = - inv(A) * deriv(A) * inv(A)

i.e., compute the inverse of the plain matrix, then compute the expression above.

chtz
  • 17,329
  • 4
  • 26
  • 56
  • Thanks, I have already been thinking about using that identity. It kind of circumnavigates the problem of evaluating the matrix inverse without branching instead of solving it. In my case the A matrix is a nonlinear function of the vector x: `A(x)`. where x is the vector to which I want to differentiate A. But I could just auto differentiate `dA(x)/dx` using CppAD and use that in my computation. Thanks for the proposal. – Manumerous Nov 12 '21 at 23:22
1

There is a mechanical conversion from branch to no-branch for arithmetic functions.

Duplicate all the variables you use in each branch, and calculate both halves. At the end of the block, multiply the if branch by condition, and the else branch by !condition, then sum them.

Similarly for a switch, calculate all the cases, and multiply by value == case.

E.g.

Mat frob_branch(Mat a, Mat b) {
    if (a.baz()) {
        return a * b;
    } else {
        return b * a;
    }
}

becomes

Mat frob_no_branch(Mat a, Mat b) {
    auto if_true = a * b;
    auto if_false = b * a;
    bool condition = a.baz();
    return (if_true * condition) + (if_false * !condition);
}
Caleth
  • 52,200
  • 2
  • 44
  • 75