search

gonum / matrix

Matrix packages for the Go language [DEPRECATED]
446 stars 53 forks source link

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

Closed ChristopherRabotin closed 7 years ago

ChristopherRabotin commented 7 years ago

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 commented 7 years ago

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 commented 7 years ago

Thanks for spotting this, @ChristopherRabotin

ChristopherRabotin commented 7 years ago

[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 commented 7 years ago

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" notifications@github.com 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 https://github.com/gonum/matrix/issues/410#issuecomment-267657188, or mute the thread https://github.com/notifications/unsubscribe-auth/ABwJvkAK8c7ApC_WY0iE4E52hhmg_pbWks5rItJjgaJpZM4LJ3Wt .

ChristopherRabotin commented 7 years ago

Sure.

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 commented 7 years ago

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 commented 7 years ago

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 notifications@github.com 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 https://github.com/gonum/matrix/issues/410#issuecomment-267667830, or mute the thread https://github.com/notifications/unsubscribe-auth/AEma6CKNmgvWYxg-Kaoh9WoH45x4WpC6ks5rIt2dgaJpZM4LJ3Wt .

btracey commented 7 years ago

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

ChristopherRabotin commented 7 years ago

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 notifications@github.com 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 https://github.com/gonum/matrix/issues/410#issuecomment-267673740, or mute the thread https://github.com/notifications/unsubscribe-auth/AEma6OMYKjHPW1MgXDx1XWVVl2mga555ks5rIuPfgaJpZM4LJ3Wt .

btracey commented 7 years ago

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