Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

mat: rewrite extractions and solves as methods on factorising types #113

Merged
merged 7 commits into from
Jul 6, 2017
Merged

mat: rewrite extractions and solves as methods on factorising types #113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
b568539
mat: rewrite cholesky extractions as To rather than From
kortschak Jun 30, 2017
3af0e0a
mat: move SolveCholesky* onto Cholesky
kortschak Jun 30, 2017
a55bda3
mat: move SolveLQ* onto LQ
kortschak Jun 30, 2017
2e3fcb1
mat: move SolveLU* onto LU
kortschak Jun 30, 2017
f6ec77b
mat: move SolveQR* onto QR
kortschak Jun 30, 2017
0387bc9
mat: move InverseCholesky onto Cholesky
kortschak Jun 30, 2017
fa212df
mat: allow nil destination for Cholesky extraction
kortschak Jul 1, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
108 changes: 63 additions & 45 deletions mat/cholesky.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -158,13 +158,13 @@ func (c *Cholesky) LogDet() float64 {
return det
}

// SolveCholesky finds the matrix m that solves A * m = b where A is represented
// by the Cholesky decomposition, placing the result in the receiver.
func (m *Dense) SolveCholesky(chol *Cholesky, b Matrix) error {
if !chol.valid() {
// Solve finds the matrix m that solves A * m = b where A is represented
// by the Cholesky decomposition, placing the result in m.
func (c *Cholesky) Solve(m *Dense, b Matrix) error {
if !c.valid() {
panic(badCholesky)
}
n := chol.chol.mat.N
n := c.chol.mat.N
bm, bn := b.Dims()
if n != bm {
panic(ErrShape)
Expand All @@ -174,10 +174,10 @@ func (m *Dense) SolveCholesky(chol *Cholesky, b Matrix) error {
if b != m {
m.Copy(b)
}
blas64.Trsm(blas.Left, blas.Trans, 1, chol.chol.mat, m.mat)
blas64.Trsm(blas.Left, blas.NoTrans, 1, chol.chol.mat, m.mat)
if chol.cond > ConditionTolerance {
return Condition(chol.cond)
blas64.Trsm(blas.Left, blas.Trans, 1, c.chol.mat, m.mat)
blas64.Trsm(blas.Left, blas.NoTrans, 1, c.chol.mat, m.mat)
if c.cond > ConditionTolerance {
return Condition(c.cond)
}
return nil
}
Expand Down Expand Up @@ -205,13 +205,13 @@ func (m *Dense) solveTwoChol(a, b *Cholesky) error {
return nil
}

// SolveCholeskyVec finds the vector v that solves A * v = b where A is represented
// by the Cholesky decomposition, placing the result in the receiver.
func (v *Vector) SolveCholeskyVec(chol *Cholesky, b *Vector) error {
if !chol.valid() {
// SolveVec finds the vector v that solves A * v = b where A is represented
// by the Cholesky decomposition, placing the result in v.
func (c *Cholesky) SolveVec(v, b *Vector) error {
if !c.valid() {
panic(badCholesky)
}
n := chol.chol.mat.N
n := c.chol.mat.N
vn := b.Len()
if vn != n {
panic(ErrShape)
Expand All @@ -223,69 +223,87 @@ func (v *Vector) SolveCholeskyVec(chol *Cholesky, b *Vector) error {
if v != b {
v.CopyVec(b)
}
blas64.Trsv(blas.Trans, chol.chol.mat, v.mat)
blas64.Trsv(blas.NoTrans, chol.chol.mat, v.mat)
if chol.cond > ConditionTolerance {
return Condition(chol.cond)
blas64.Trsv(blas.Trans, c.chol.mat, v.mat)
blas64.Trsv(blas.NoTrans, c.chol.mat, v.mat)
if c.cond > ConditionTolerance {
return Condition(c.cond)
}
return nil

}

// UFromCholesky extracts the n×n upper triangular matrix U from a Cholesky
// decomposition
// UTo extracts the n×n upper triangular matrix U from a Cholesky
// decomposition into dst and returns the result. If dst is nil a new
// TriDense is allocated.
// A = U^T * U.
func (t *TriDense) UFromCholesky(chol *Cholesky) {
if !chol.valid() {
func (c *Cholesky) UTo(dst *TriDense) *TriDense {
if !c.valid() {
panic(badCholesky)
}
n := chol.chol.mat.N
t.reuseAs(n, Upper)
t.Copy(chol.chol)
n := c.chol.mat.N
if dst == nil {
dst = NewTriDense(n, Upper, make([]float64, n*n))
} else {
dst.reuseAs(n, Upper)
}
dst.Copy(c.chol)
return dst
}

// LFromCholesky extracts the n×n lower triangular matrix L from a Cholesky
// decomposition
// LTo extracts the n×n lower triangular matrix L from a Cholesky
// decomposition into dst and returns the result. If dst is nil a new
// TriDense is allocated.
// A = L * L^T.
func (t *TriDense) LFromCholesky(chol *Cholesky) {
if !chol.valid() {
func (c *Cholesky) LTo(dst *TriDense) *TriDense {
if !c.valid() {
panic(badCholesky)
}
n := chol.chol.mat.N
t.reuseAs(n, Lower)
t.Copy(chol.chol.TTri())
n := c.chol.mat.N
if dst == nil {
dst = NewTriDense(n, Lower, make([]float64, n*n))
} else {
dst.reuseAs(n, Lower)
}
dst.Copy(c.chol.TTri())
return dst
}

// FromCholesky reconstructs the original positive definite matrix given its
// Cholesky decomposition.
func (s *SymDense) FromCholesky(chol *Cholesky) {
if !chol.valid() {
// To reconstructs the original positive definite matrix given its
// Cholesky decomposition into dst and returns the result. If dst is nil
// a new SymDense is allocated.
func (c *Cholesky) To(dst *SymDense) *SymDense {
if !c.valid() {
panic(badCholesky)
}
n := chol.chol.mat.N
s.reuseAs(n)
s.SymOuterK(1, chol.chol.T())
n := c.chol.mat.N
if dst == nil {
dst = NewSymDense(n, make([]float64, n*n))
} else {
dst.reuseAs(n)
}
dst.SymOuterK(1, c.chol.T())
return dst
}

// InverseCholesky computes the inverse of the matrix represented by its Cholesky
// factorization and stores the result into the receiver. If the factorized
// InverseTo computes the inverse of the matrix represented by its Cholesky
// factorization and stores the result into s. If the factorized
// matrix is ill-conditioned, a Condition error will be returned.
// Note that matrix inversion is numerically unstable, and should generally be
// avoided where possible, for example by using the Solve routines.
func (s *SymDense) InverseCholesky(chol *Cholesky) error {
if !chol.valid() {
func (c *Cholesky) InverseTo(s *SymDense) error {
if !c.valid() {
panic(badCholesky)
}
// TODO(btracey): Replace this code with a direct call to Dpotri when it
// is available.
s.reuseAs(chol.chol.mat.N)
s.reuseAs(c.chol.mat.N)
// If:
// chol(A) = U^T * U
// Then:
// chol(A^-1) = S * S^T
// where S = U^-1
var t TriDense
err := t.InverseTri(chol.chol)
err := t.InverseTri(c.chol)
s.SymOuterK(1, &t)
return err
}
Expand Down
12 changes: 5 additions & 7 deletions mat/cholesky_example_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,17 +35,16 @@ func ExampleCholesky() {
// Use the factorization to solve the system of equations a * x = b.
b := mat.NewVector(4, []float64{1, 2, 3, 4})
var x mat.Vector
if err := x.SolveCholeskyVec(&chol, b); err != nil {
if err := chol.SolveVec(&x, b); err != nil {
fmt.Println("Matrix is near singular: ", err)
}
fmt.Println("Solve a * x = b")
fmt.Printf("x = %0.4v\n", mat.Formatted(&x, mat.Prefix(" ")))

// Extract the factorization and check that it equals the original matrix.
var t mat.TriDense
t.LFromCholesky(&chol)
t := chol.LTo(nil)
var test mat.Dense
test.Mul(&t, t.T())
test.Mul(t, t.T())
fmt.Println()
fmt.Printf("L * L^T = %0.4v\n", mat.Formatted(&a, mat.Prefix(" ")))

Expand Down Expand Up @@ -92,13 +91,12 @@ func ExampleCholesky_SymRankOne() {
// Rank-1 update the matrix a.
a.SymRankOne(a, 1, x)

var au mat.SymDense
au.FromCholesky(&chol)
au := chol.To(nil)

// Print the matrix that was updated directly.
fmt.Printf("\nA' = %0.4v\n", mat.Formatted(a, mat.Prefix(" ")))
// Print the matrix recovered from the factorization.
fmt.Printf("\nU'^T * U' = %0.4v\n", mat.Formatted(&au, mat.Prefix(" ")))
fmt.Printf("\nU'^T * U' = %0.4v\n", mat.Formatted(au, mat.Prefix(" ")))

// Output:
// A = ⎡ 1 1 1 1⎤
Expand Down
31 changes: 14 additions & 17 deletions mat/cholesky_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,18 +53,16 @@ func TestCholesky(t *testing.T) {
if math.Abs(test.cond-chol.cond) > 1e-13 {
t.Errorf("Condition number mismatch: Want %v, got %v", test.cond, chol.cond)
}
var U TriDense
U.UFromCholesky(chol)
U := chol.UTo(nil)
aCopy := DenseCopyOf(test.a)
var a Dense
a.Mul(U.TTri(), &U)
a.Mul(U.TTri(), U)
if !EqualApprox(&a, aCopy, 1e-14) {
t.Error("unexpected Cholesky factor product")
}

var L TriDense
L.LFromCholesky(chol)
a.Mul(&L, L.TTri())
L := chol.LTo(nil)
a.Mul(L, L.TTri())
if !EqualApprox(&a, aCopy, 1e-14) {
t.Error("unexpected Cholesky factor product")
}
Expand Down Expand Up @@ -103,7 +101,7 @@ func TestCholeskySolve(t *testing.T) {
}

var x Dense
x.SolveCholesky(&chol, test.b)
chol.Solve(&x, test.b)
if !EqualApprox(&x, test.ans, 1e-12) {
t.Error("incorrect Cholesky solve solution")
}
Expand Down Expand Up @@ -208,7 +206,7 @@ func TestCholeskySolveVec(t *testing.T) {
}

var x Vector
x.SolveCholeskyVec(&chol, test.b)
chol.SolveVec(&x, test.b)
if !EqualApprox(&x, test.ans, 1e-12) {
t.Error("incorrect Cholesky solve solution")
}
Expand All @@ -221,7 +219,7 @@ func TestCholeskySolveVec(t *testing.T) {
}
}

func TestFromCholesky(t *testing.T) {
func TestCholeskyTo(t *testing.T) {
for _, test := range []*SymDense{
NewSymDense(3, []float64{
53, 59, 37,
Expand All @@ -234,11 +232,10 @@ func TestFromCholesky(t *testing.T) {
if !ok {
t.Fatal("unexpected Cholesky factorization failure: not positive definite")
}
var s SymDense
s.FromCholesky(&chol)
s := chol.To(nil)

if !EqualApprox(&s, test, 1e-12) {
t.Errorf("Cholesky reconstruction not equal to original matrix.\nWant:\n% v\nGot:\n% v\n", Formatted(test), Formatted(&s))
if !EqualApprox(s, test, 1e-12) {
t.Errorf("Cholesky reconstruction not equal to original matrix.\nWant:\n% v\nGot:\n% v\n", Formatted(test), Formatted(s))
}
}
}
Expand Down Expand Up @@ -279,7 +276,7 @@ func TestCloneCholesky(t *testing.T) {
}
}

func TestInverseCholesky(t *testing.T) {
func TestCholeskyInverseTo(t *testing.T) {
for _, n := range []int{1, 3, 5, 9} {
data := make([]float64, n*n)
for i := range data {
Expand All @@ -295,7 +292,7 @@ func TestInverseCholesky(t *testing.T) {
}

var sInv SymDense
sInv.InverseCholesky(&chol)
chol.InverseTo(&sInv)

var ans Dense
ans.Mul(&sInv, &s)
Expand Down Expand Up @@ -341,7 +338,7 @@ func TestCholeskySymRankOne(t *testing.T) {
a.SymRankOne(&a, alpha, x)

var achol SymDense
achol.FromCholesky(&chol)
chol.To(&achol)
if !EqualApprox(&achol, &a, 1e-13) {
t.Errorf("n=%v, alpha=%v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v",
n, alpha, Formatted(&a), Formatted(&achol))
Expand Down Expand Up @@ -416,7 +413,7 @@ func TestCholeskySymRankOne(t *testing.T) {
a.SymRankOne(a, test.alpha, x)

var achol SymDense
achol.FromCholesky(&chol)
chol.To(&achol)
if !EqualApprox(&achol, a, 1e-13) {
t.Errorf("Case %v: mismatch between updated matrix and from Cholesky:\nupdated:\n%v\nfrom Cholesky:\n%v",
i, Formatted(a), Formatted(&achol))
Expand Down
14 changes: 7 additions & 7 deletions mat/lq.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,16 +121,16 @@ func (lq *LQ) QTo(dst *Dense) *Dense {
return dst
}

// SolveLQ finds a minimum-norm solution to a system of linear equations defined
// Solve finds a minimum-norm solution to a system of linear equations defined
// by the matrices A and b, where A is an m×n matrix represented in its LQ factorized
// form. If A is singular or near-singular a Condition error is returned. Please
// see the documentation for Condition for more information.
//
// The minimization problem solved depends on the input parameters.
// If trans == false, find the minimum norm solution of A * X = b.
// If trans == true, find X such that ||A*X - b||_2 is minimized.
// The solution matrix, X, is stored in place into the receiver.
func (m *Dense) SolveLQ(lq *LQ, trans bool, b Matrix) error {
// The solution matrix, X, is stored in place into m.
func (lq *LQ) Solve(m *Dense, trans bool, b Matrix) error {
r, c := lq.lq.Dims()
br, bc := b.Dims()

Expand Down Expand Up @@ -188,9 +188,9 @@ func (m *Dense) SolveLQ(lq *LQ, trans bool, b Matrix) error {
return nil
}

// SolveLQVec finds a minimum-norm solution to a system of linear equations.
// Please see Dense.SolveLQ for the full documentation.
func (v *Vector) SolveLQVec(lq *LQ, trans bool, b *Vector) error {
// SolveVec finds a minimum-norm solution to a system of linear equations.
// Please see LQ.Solve for the full documentation.
func (lq *LQ) SolveVec(v *Vector, trans bool, b *Vector) error {
if v != b {
v.checkOverlap(b.mat)
}
Expand All @@ -202,5 +202,5 @@ func (v *Vector) SolveLQVec(lq *LQ, trans bool, b *Vector) error {
} else {
v.reuseAs(c)
}
return v.asDense().SolveLQ(lq, trans, b.asDense())
return lq.Solve(v.asDense(), trans, b.asDense())
}
8 changes: 4 additions & 4 deletions mat/lq_test.go
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ func TestSolveLQ(t *testing.T) {
var x Dense
lq := &LQ{}
lq.Factorize(a)
x.SolveLQ(lq, trans, b)
lq.Solve(&x, trans, b)

// Test that the normal equations hold.
// A^T * A * x = A^T * b if !trans
Expand Down Expand Up @@ -130,7 +130,7 @@ func TestSolveLQVec(t *testing.T) {
var x Vector
lq := &LQ{}
lq.Factorize(a)
x.SolveLQVec(lq, trans, b)
lq.SolveVec(&x, trans, b)

// Test that the normal equations hold.
// A^T * A * x = A^T * b if !trans
Expand Down Expand Up @@ -166,13 +166,13 @@ func TestSolveLQCond(t *testing.T) {
lq.Factorize(test)
b := NewDense(m, 2, nil)
var x Dense
if err := x.SolveLQ(&lq, false, b); err == nil {
if err := lq.Solve(&x, false, b); err == nil {
t.Error("No error for near-singular matrix in matrix solve.")
}

bvec := NewVector(m, nil)
var xvec Vector
if err := xvec.SolveLQVec(&lq, false, bvec); err == nil {
if err := lq.SolveVec(&xvec, false, bvec); err == nil {
t.Error("No error for near-singular matrix in matrix solve.")
}
}
Expand Down