Add condition number check to LU solve

mat64/lu.go

@@ -17,6 +17,29 @@ import (
 type LU struct {
 	lu    *Dense
 	pivot []int
+	cond  float64
+}
+
+// updateCond updates the stored condition number of the matrix. Norm is the
+// norm of the original matrix. If norm is negative it will be estimated.
+func (lu *LU) updateCond(norm float64) {
+	n := lu.lu.mat.Cols
+	work := make([]float64, 4*n)
+	iwork := make([]int, n)
+	if norm < 0 {
+		// This is an approximation. By the defintion of a norm, ||AB|| <= ||A|| ||B||.
+		// The condition number is ||A|| || A^-1||, so this will underestimate
+		// the condition number somewhat.
+		// The norm of the original factorized matrix cannot be stored because of
+		// update possibilites, e.g. RankOne.
+		u := lu.lu.asTriDense(n, blas.NonUnit, blas.Upper)
+		l := lu.lu.asTriDense(n, blas.Unit, blas.Lower)
+		unorm := lapack64.Lantr(condNorm, u.mat, work)
+		lnorm := lapack64.Lantr(condNorm, l.mat, work)
+		norm = unorm * lnorm
+	}
+	v := lapack64.Gecon(condNorm, lu.lu.mat, norm, work, iwork)
+	lu.cond = 1 / v
 }
 
 // Factorize computes the LU factorization of the square matrix a and stores the
@@ -39,7 +62,10 @@ func (lu *LU) Factorize(a Matrix) {
 		lu.pivot = make([]int, r)
 	}
 	lu.pivot = lu.pivot[:r]
+	work := make([]float64, r)
+	anorm := lapack64.Lange(condNorm, lu.lu.mat, work)
 	lapack64.Getrf(lu.lu.mat, lu.pivot)
+	lu.updateCond(anorm)
 }
 
 // Det returns the determinant of the matrix that has been factorized.
@@ -135,6 +161,7 @@ func (lu *LU) RankOne(orig *LU, alpha float64, x, y *Vector) {
 			lum.Data[j*lum.Stride+i] += gamma * tmp
 		}
 	}
+	lu.updateCond(-1)
 }
 
 // LFromLU extracts the lower triangular matrix from an LU factorization.
@@ -215,6 +242,9 @@ func (m *Dense) SolveLU(lu *LU, trans bool, b Matrix) error {
 		t = blas.Trans
 	}
 	lapack64.Getrs(t, lu.lu.mat, m.mat, lu.pivot)
+	if lu.cond > condTol {
+		return Condition(lu.cond)
+	}
 	return nil
 }
 
@@ -257,5 +287,8 @@ func (v *Vector) SolveLUVec(lu *LU, trans bool, b *Vector) error {
 		t = blas.Trans
 	}
 	lapack64.Getrs(t, lu.lu.mat, vMat, lu.pivot)
+	if lu.cond > condTol {
+		return Condition(lu.cond)
+	}
 	return nil
 }

mat64/lu_test.go

@@ -143,6 +143,27 @@ func (s *S) TestSolveLU(c *check.C) {
 	// TODO(btracey): Add testOneInput test when such a function exists.
 }
 
+func (s *S) TestSolveLUCond(c *check.C) {
+	for _, test := range []*Dense{
+		NewDense(2, 2, []float64{1, 0, 0, 1e-20}),
+	} {
+		m, _ := test.Dims()
+		var lu LU
+		lu.Factorize(test)
+		b := NewDense(m, 2, nil)
+		var x Dense
+		if err := x.SolveLU(&lu, false, b); err == nil {
+			c.Error("No error for near-singular matrix in matrix solve.")
+		}
+
+		bvec := NewVector(m, nil)
+		var xvec Vector
+		if err := xvec.SolveLUVec(&lu, false, bvec); err == nil {
+			c.Error("No error for near-singular matrix in matrix solve.")
+		}
+	}
+}
+
 func (s *S) TestSolveLUVec(c *check.C) {
 	for _, n := range []int{5, 10} {
 		a := NewDense(n, n, nil)

mat64/matrix.go

@@ -8,6 +8,7 @@ import (
 	"fmt"
 
 	"github.com/gonum/blas/blas64"
+	"github.com/gonum/lapack"
 )
 
 // Matrix is the basic matrix interface type.
@@ -390,12 +391,12 @@ func MaybeFloat(fn func() float64) (f float64, err error) {
 	return fn(), nil
 }
 
-// Condition is the reciprocal of the condition number of a matrix. The condition
+// Condition is the condition number of a matrix. The condition
 // number is defined as ||A|| * ||A^-1||.
 //
 // One important use of Condition is during linear solve routines (finding x such
 // that A * x = b). The condition number of A indicates the accuracy of
-// the computed solution. A Condition error will be returned if the inverse condition
+// the computed solution. A Condition error will be returned if the condition
 // number of A is sufficiently large. If A is exactly singular to working precision,
 // Condition == ∞, and the solve algorithm may have completed early. If Condition
 // is large and finite the solve algorithm will be performed, but the computed
@@ -407,6 +408,13 @@ func (c Condition) Error() string {
 	return fmt.Sprintf("matrix singular or near-singular with inverse condition number %.4e", c)
 }
 
+// condTol describes the limit of the condition number. If the inverse of the
+// condition number is above this value, the matrix is considered singular.
+var condTol float64 = 1e16
+
+// condNorm describes the matrix norm to use for computing the condition number.
+var condNorm = lapack.MaxRowSum
+
 // Type Error represents matrix handling errors. These errors can be recovered by Maybe wrappers.
 type Error struct{ string }