Bug in inverse of 4x4 matrix returning NaN

[TsurHerman]

This came up when working with coordinate transformation matrices in openGL

MWE

julia> using StaticArrays

julia> A =@SMatrix [1 0 0 0;0 0 1 0;0 1 0 0;0 0 0 1]
4×4 StaticArrays.SArray{Tuple{4,4},Int64,2,16}:
 1  0  0  0
 0  0  1  0
 0  1  0  0
 0  0  0  1

julia> inv(A)
4×4 StaticArrays.SArray{Tuple{4,4},Float64,2,16}:
 NaN  NaN  NaN  NaN
 NaN  NaN  NaN  NaN
 NaN  NaN  NaN  NaN
 NaN  NaN  NaN  NaN

[tkoolen]

Looks like the fact that the method doesn't work for the case when then upper left 2x2 block is singular (as it is in the @TsurHerman's example) is acknowledged in http://www.freevec.org/function/inverse_matrix_4x4_using_partitioning, referenced in the code:

Well, in that particular case, we have to call a generic -but slower- function

[TsurHerman]

The following code borrowed from @SimonDanisch FixedSizeArrays,and slightly modified to shave off another 15% benchmarks at 50ns on my laptop whereas the current implementation benchmarks at 30ns.

@inline det4x4(A) = @inbounds return (
    A[13] * A[10]  * A[7]  * A[4]  - A[9] * A[14] * A[7]  * A[4]   -
    A[13] * A[6]   * A[11] * A[4]  + A[5] * A[14] * A[11] * A[4]   +
    A[9]  * A[6]   * A[15] * A[4]  - A[5] * A[10] * A[15] * A[4]   -
    A[13] * A[10]  * A[3]  * A[8]  + A[9] * A[14] * A[3]  * A[8]   +
    A[13] * A[2]   * A[11] * A[8]  - A[1] * A[14] * A[11] * A[8]   -
    A[9]  * A[2]   * A[15] * A[8]  + A[1] * A[10] * A[15] * A[8]   +
    A[13] * A[6]   * A[3]  * A[12] - A[5] * A[14] * A[3]  * A[12]  -
    A[13] * A[2]   * A[7]  * A[12] + A[1] * A[14] * A[7]  * A[12]  +
    A[5]  * A[2]   * A[15] * A[12] - A[1] * A[6]  * A[15] * A[12]  -
    A[9]  * A[6]   * A[3]  * A[16] + A[5] * A[10] * A[3]  * A[16]  +
    A[9]  * A[2]   * A[7]  * A[16] - A[1] * A[10] * A[7]  * A[16]  -
    A[5]  * A[2]   * A[11] * A[16] + A[1] * A[6]  * A[11] * A[16]
)
@inline function inv4x4(A)
    ideterminant = 1/det4x4(A)
    @inbounds B =  @SMatrix [
        (A[2,3]*A[3,4]*A[4,2] - A[2,4]*A[3,3]*A[4,2] + A[2,4]*A[3,2]*A[4,3] - A[2,2]*A[3,4]*A[4,3] - A[2,3]*A[3,2]*A[4,4] + A[2,2]*A[3,3]*A[4,4]) * ideterminant
        (A[2,4]*A[3,3]*A[4,1] - A[2,3]*A[3,4]*A[4,1] - A[2,4]*A[3,1]*A[4,3] + A[2,1]*A[3,4]*A[4,3] + A[2,3]*A[3,1]*A[4,4] - A[2,1]*A[3,3]*A[4,4]) * ideterminant
        (A[2,2]*A[3,4]*A[4,1] - A[2,4]*A[3,2]*A[4,1] + A[2,4]*A[3,1]*A[4,2] - A[2,1]*A[3,4]*A[4,2] - A[2,2]*A[3,1]*A[4,4] + A[2,1]*A[3,2]*A[4,4]) * ideterminant
        (A[2,3]*A[3,2]*A[4,1] - A[2,2]*A[3,3]*A[4,1] - A[2,3]*A[3,1]*A[4,2] + A[2,1]*A[3,3]*A[4,2] + A[2,2]*A[3,1]*A[4,3] - A[2,1]*A[3,2]*A[4,3]) * ideterminant

        (A[1,4]*A[3,3]*A[4,2] - A[1,3]*A[3,4]*A[4,2] - A[1,4]*A[3,2]*A[4,3] + A[1,2]*A[3,4]*A[4,3] + A[1,3]*A[3,2]*A[4,4] - A[1,2]*A[3,3]*A[4,4]) * ideterminant
        (A[1,3]*A[3,4]*A[4,1] - A[1,4]*A[3,3]*A[4,1] + A[1,4]*A[3,1]*A[4,3] - A[1,1]*A[3,4]*A[4,3] - A[1,3]*A[3,1]*A[4,4] + A[1,1]*A[3,3]*A[4,4]) * ideterminant
        (A[1,4]*A[3,2]*A[4,1] - A[1,2]*A[3,4]*A[4,1] - A[1,4]*A[3,1]*A[4,2] + A[1,1]*A[3,4]*A[4,2] + A[1,2]*A[3,1]*A[4,4] - A[1,1]*A[3,2]*A[4,4]) * ideterminant
        (A[1,2]*A[3,3]*A[4,1] - A[1,3]*A[3,2]*A[4,1] + A[1,3]*A[3,1]*A[4,2] - A[1,1]*A[3,3]*A[4,2] - A[1,2]*A[3,1]*A[4,3] + A[1,1]*A[3,2]*A[4,3]) * ideterminant

        (A[1,3]*A[2,4]*A[4,2] - A[1,4]*A[2,3]*A[4,2] + A[1,4]*A[2,2]*A[4,3] - A[1,2]*A[2,4]*A[4,3] - A[1,3]*A[2,2]*A[4,4] + A[1,2]*A[2,3]*A[4,4]) * ideterminant
        (A[1,4]*A[2,3]*A[4,1] - A[1,3]*A[2,4]*A[4,1] - A[1,4]*A[2,1]*A[4,3] + A[1,1]*A[2,4]*A[4,3] + A[1,3]*A[2,1]*A[4,4] - A[1,1]*A[2,3]*A[4,4]) * ideterminant
        (A[1,2]*A[2,4]*A[4,1] - A[1,4]*A[2,2]*A[4,1] + A[1,4]*A[2,1]*A[4,2] - A[1,1]*A[2,4]*A[4,2] - A[1,2]*A[2,1]*A[4,4] + A[1,1]*A[2,2]*A[4,4]) * ideterminant
        (A[1,3]*A[2,2]*A[4,1] - A[1,2]*A[2,3]*A[4,1] - A[1,3]*A[2,1]*A[4,2] + A[1,1]*A[2,3]*A[4,2] + A[1,2]*A[2,1]*A[4,3] - A[1,1]*A[2,2]*A[4,3]) * ideterminant

        (A[1,4]*A[2,3]*A[3,2] - A[1,3]*A[2,4]*A[3,2] - A[1,4]*A[2,2]*A[3,3] + A[1,2]*A[2,4]*A[3,3] + A[1,3]*A[2,2]*A[3,4] - A[1,2]*A[2,3]*A[3,4]) * ideterminant
        (A[1,3]*A[2,4]*A[3,1] - A[1,4]*A[2,3]*A[3,1] + A[1,4]*A[2,1]*A[3,3] - A[1,1]*A[2,4]*A[3,3] - A[1,3]*A[2,1]*A[3,4] + A[1,1]*A[2,3]*A[3,4]) * ideterminant
        (A[1,4]*A[2,2]*A[3,1] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] + A[1,1]*A[2,4]*A[3,2] + A[1,2]*A[2,1]*A[3,4] - A[1,1]*A[2,2]*A[3,4]) * ideterminant
        (A[1,2]*A[2,3]*A[3,1] - A[1,3]*A[2,2]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] + A[1,1]*A[2,2]*A[3,3]) * ideterminant]
        SMatrix{4,4,eltype(B),16}(B)
end
 

[c42f]

Right, those darn sets of measure zero, they get you every time :-/ This is definitely an embarrassing oversight, thanks for the bug report!

Another option to fix this (rather than the determinant based method) is to add a kind of pivoting to the block based partition - choose a permutation of the rows to make the matrices in the partition well conditioned.

But I do wonder whether a static LU would just be a better option overall.

[wsshin]

I like the static LU idea. That way, we will be able to mimic inv of Base, which is implemented as follows and uses lufact:

function inv(A::StridedMatrix{T}) where T
    checksquare(A)
    S = typeof((one(T)*zero(T) + one(T)*zero(T))/one(T))
    AA = convert(AbstractArray{S}, A)
    if istriu(AA)
        Ai = inv(UpperTriangular(AA))
    elseif istril(AA)
        Ai = inv(LowerTriangular(AA))
    else
        Ai = inv(lufact(AA))
    end
    return convert(typeof(parent(Ai)), Ai)
end

[c42f]

@tkoolen I'd be fine with using the adjugate based version you've copied across as a stopgap until someone has time for #73. It's got some stability problems of its own, but at least it won't fall apart for perfectly conditioned permutation matrices.

[tkoolen]

@TsurHerman posted the adjugate-based version, but yes, I agree with that approach.

[c42f]

My apologies, thanks @TsurHerman

[TsurHerman]

Its o.k , I am trying to issue a PR but there is an error in the tests due to the fact that det4x4
for integer (Int64) returns Int64

one possible solution is to relax
@test det(@SMatrix [1 2 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) === 400.0

to

@test det(@SMatrix [1 2 3 4; 0 5 6 7; 0 0 8 9; 0 0 0 10]) == 400.0

Another issue is that for Uint64 the result is a 2-complement of the final result.
So special treatment for unsigned integers is needed, but the current solution of converting to Signed
is less then satisfactory especially for 8-bit and 16-bit unsigned integers which are more likely to break due to inexact error
Signed(UInt8(0xff)) # yields Inexact error

one last thing , is that det2x2 det3x3 and det4x4 can probably nicely share code if we use
https://en.wikipedia.org/wiki/Leibniz_formula_for_determinants

and the same also holds for inverse I presume , but really combinatorics is not my strongest side.

So the minimal fix for now, would be to implement it as is, relax the test for determinant of size 4x4
And later issue another PR for promoting everything unsigned to Int64 and not just Signed

[c42f]

In my opinion the special code for Unsigned in various places is a bit of an anomaly - I'd be more than happy to replace it with a more systematic solution. Also I doubt the return type of det is critical.

[TsurHerman]

#260

[c42f]

Fixed in #260