matrix/mat64: Dense.Inverse does not fail on singular matrix

[ChristopherRabotin]

One should expect Dense.Inverse to return an error for singular matrices. Instead, it returns garbage: for a 2x2, it returns the input matrix, for a 4x4 it returns the input matrix where some values are changed. I have not checked for other sizes.

Code:

package main

import (
	"fmt"

	"github.com/gonum/matrix/mat64"
)

func main() {
	A := mat64.NewDense(2, 2, []float64{0, 0, 0, 1})
	var InvA mat64.Dense
	err := InvA.Inverse(A)
	if err != nil {
		panic(err)
	}
	fmt.Printf("%v\n", mat64.Formatted(&InvA, mat64.Prefix("")))

	B := mat64.NewDense(4, 4, []float64{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 20, 20, 0, 0, 20, 20})
	var InvB mat64.Dense
	err = InvB.Inverse(B)
	if err != nil {
		panic(err)
	}
	fmt.Printf("%v\n", mat64.Formatted(&InvB, mat64.Prefix("")))
}

Result:

⎡0  0⎤
⎣0  1⎦
⎡ 0   0   0   0⎤
⎢ 0   0   0   0⎥
⎢ 0   0  20  20⎥
⎣ 0   0   1   0⎦

[vladimir-ch]

Yes, that should return an error. I've submitted a possible fix in #411 .
At the moment the receiver is always modified even if an error is returned.

[vladimir-ch]

Thanks for spotting this, @ChristopherRabotin

[ChristopherRabotin]

[Issue needs to be re-opened]
Sorry to be the bearer of bad news, but now, the Inverse now fails on matrices with small-ish numbers, although it really should not. The matrix used in this example is invertible on Matlab up until the exponent -15. I also wonder whether this isn't an issue for lapack instead.

package main

import (
	"fmt"
	"math"

	"github.com/gonum/matrix/mat64"
)

func main() {
	//Q := mat64.NewSymDense(3, []float64{25, 0, 0, 0, 25, 0, 0, 0, 25})
	Q := mat64.NewSymDense(3, []float64{2.5e-9, 6.25e-7, (25e-5) / 3, 6.25e-7, (5e-1) / 3, 2.5e-2, (25e-5) / 3, 2.5e-2, 5})
	for exp := 0; exp < 6; exp++ {
		var QPrime, Qinv mat64.Dense
		QPrime.Scale(math.Pow(10, -1.0*float64(exp)), Q)
		if err := Qinv.Inverse(&QPrime); err != nil {
			fmt.Printf("Q'=%v\n", mat64.Formatted(&QPrime, mat64.Prefix("   ")))
			panic(fmt.Errorf("Q not invertible (exp=-%d): %s", exp, err))
		}
	}
	fmt.Println("ok")
}

:-( I'll just put a sad face here because all my tests and code on my Kalman fitler is now broken...

[vladimir-ch]

Can you check what condition numbers Inverse returns for your matrices? Can you check whether your matrix times the computed inverse gives something similar to the identity?

…

On Dec 16, 2016 7:03 PM, "Chris" ***@***.***> wrote: Sorry to be the bearer of bad news, but now, the Inverse now fails on matrices with small-ish numbers, although it really should not. The matrix used in this example is invertible on Matlab up until the exponent -15. I also wonder whether this isn't an issue for lapack instead. package main import ( "fmt" "math" "github.com/gonum/matrix/mat64" ) func main() { //Q := mat64.NewSymDense(3, []float64{25, 0, 0, 0, 25, 0, 0, 0, 25}) Q := mat64.NewSymDense(3, []float64{2.5e-9, 6.25e-7, (25e-5) / 3, 6.25e-7, (5e-1) / 3, 2.5e-2, (25e-5) / 3, 2.5e-2, 5}) for exp := 0; exp < 6; exp++ { var QPrime, Qinv mat64.Dense QPrime.Scale(math.Pow(10, -1.0*float64(exp)), Q) if err := Qinv.Inverse(&QPrime); err != nil { fmt.Printf("Q'=%v\n", mat64.Formatted(&QPrime, mat64.Prefix(" "))) panic(fmt.Errorf("Q not invertible (exp=-%d): %s", exp, err)) } } fmt.Println("ok") } :-( I'll just put a sad face here because all my tests and code on my Kalman fitler is now broken... — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub <#410 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/ABwJvkAK8c7ApC_WY0iE4E52hhmg_pbWks5rItJjgaJpZM4LJ3Wt> .

[ChristopherRabotin]

Sure.

• The condition number is 5.5370e+16. Here is the full error from the above code:
  panic: Q not invertible (exp=-3): matrix singular or near-singular with condition number 5.5370e+16
• I'm surprised that the inverse is actually mostly computed (although I didn't delve into the code, I expected that an error during the inverse computation would not populate the matrix). Regardless, yes, the result is somewhat close to identity:

QiQ=⎡                      1  -3.1763735522036263e-22                        0⎤
    ⎢                      0       0.9999999999999999                        0⎥
    ⎣ 1.4551915228366852e-11  -2.0816681711721685e-17                        1⎦

For reference, here is the updated test code (it's not very pretty):

package main

import (
	"fmt"
	"math"

	"github.com/gonum/matrix/mat64"
)

func main() {
	//Q := mat64.NewSymDense(3, []float64{25, 0, 0, 0, 25, 0, 0, 0, 25})
	Q := mat64.NewSymDense(3, []float64{2.5e-9, 6.25e-7, (25e-5) / 3, 6.25e-7, (5e-1) / 3, 2.5e-2, (25e-5) / 3, 2.5e-2, 5})
	for exp := 0; exp < 6; exp++ {
		var QPrime, Qinv, ID mat64.Dense
		QPrime.Scale(math.Pow(10, -1.0*float64(exp)), Q)
		if err := Qinv.Inverse(&QPrime); err != nil {
			fmt.Printf("Q'=%v\n", mat64.Formatted(&QPrime, mat64.Prefix("   ")))
			fmt.Printf("Qi=%v\n", mat64.Formatted(&Qinv, mat64.Prefix("   ")))
			ID.Mul(&QPrime, &Qinv)
			fmt.Printf("QiQ=%v\n", mat64.Formatted(&ID, mat64.Prefix("    ")))
			panic(fmt.Errorf("Q not invertible (exp=-%d): %s", exp, err))
		}
	}
	fmt.Println("ok")
}

[btracey]

Sorry, I don't understand what you think the bug is (what you think the behavior should be).

The condition number is a measure of the expected precision of the operation. A condition number of 10^16 means that you could lose up to 16 digits of precision (i.e. everything). The matrix inverse algorithm works unless your matrix is exactly singular; the error warns you that your result may be garbage. This error can be handled, for example by checking if the resulting inverse times the resulting matrix is close to identity (for whatever definition of close you need).

[ChristopherRabotin]

Oh okay, sorry, I misunderstood. I thought that an error would be returned (ie not nil) only if the matrix was singular to float point precision.

…

On Fri, Dec 16, 2016, 11:54 Brendan Tracey ***@***.***> wrote: Sorry, I don't understand what you think the bug is (what you think the behavior should be). The condition number is a measure of the expected precision of the operation. A condition number of 10^16 means that you could lose up to 16 digits of precision (i.e. everything). The matrix inverse algorithm works unless your matrix is exactly singular, the error warns you that your result may be garbage. This error can be handled, for example by checking if the resulting inverse times the resulting matrix is close to identity (for whatever definition of close you need). — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#410 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AEma6CKNmgvWYxg-Kaoh9WoH45x4WpC6ks5rIt2dgaJpZM4LJ3Wt> .

[btracey]

No, it's like matlab in that it returns an error if it's singular or near singular

[ChristopherRabotin]

Okay. So how is it near singular with a condition at e-16 if it's a float64? It intuitively sounds to me like nearly half of the bits aren't used.

…

On Fri, Dec 16, 2016, 12:19 Brendan Tracey ***@***.***> wrote: No, it's like matlab in that it returns an error if it's singular or near singular — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#410 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AEma6OMYKjHPW1MgXDx1XWVVl2mga555ks5rIuPfgaJpZM4LJ3Wt> .

[btracey]

float64 numbers have a mantissa of 53 bits. 2^53 ~= 10^16