mat64: add higher order generalised singular value decomposition
MODIFY
mat64/cholesky.go
mat64/cholesky.go
@@ -176,6 +176,29 @@ func (m *Dense) SolveCholesky(chol *Cholesky, b Matrix) error {
return nil
}
+// TODO(kortschak): Export this as SolveTwoChol.
+// solveTwoChol finds the matrix m that solves A * m = B where A and B are represented
+// by their Cholesky decompositions a and b, placing the result in the receiver.
+func (m *Dense) solveTwoChol(a, b *Cholesky) error {
+ if !a.valid() || !b.valid() {
+ panic(badCholesky)
+ }
+ bn := b.chol.mat.N
+ if a.chol.mat.N != bn {
+ panic(matrix.ErrShape)
+ }
+
+ m.reuseAsZeroed(bn, bn)
+ m.Copy(b.chol.T())
+ blas64.Trsm(blas.Left, blas.Trans, 1, a.chol.mat, m.mat)
+ blas64.Trsm(blas.Left, blas.NoTrans, 1, a.chol.mat, m.mat)
+ blas64.Trmm(blas.Right, blas.NoTrans, 1, b.chol.mat, m.mat)
+ if a.cond > matrix.ConditionTolerance {
+ return matrix.Condition(a.cond)
+ }
+ 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 {

MODIFY
mat64/cholesky_test.go
mat64/cholesky_test.go
@@ -117,6 +117,67 @@ func TestCholeskySolve(t *testing.T) {
}
}
+func TestSolveTwoChol(t *testing.T) {
+ for _, test := range []struct {
+ a, b *SymDense
+ }{
+ {
+ a: NewSymDense(2, []float64{
+ 1, 0,
+ 0, 1,
+ }),
+ b: NewSymDense(2, []float64{
+ 1, 0,
+ 0, 1,
+ }),
+ },
+ {
+ a: NewSymDense(2, []float64{
+ 1, 0,
+ 0, 1,
+ }),
+ b: NewSymDense(2, []float64{
+ 2, 0,
+ 0, 2,
+ }),
+ },
+ {
+ a: NewSymDense(3, []float64{
+ 53, 59, 37,
+ 59, 83, 71,
+ 37, 71, 101,
+ }),
+ b: NewSymDense(3, []float64{
+ 2, -1, 0,
+ -1, 2, -1,
+ 0, -1, 2,
+ }),
+ },
+ } {
+ var chola, cholb Cholesky
+ ok := chola.Factorize(test.a)
+ if !ok {
+ t.Fatal("unexpected Cholesky factorization failure for a: not positive definite")
+ }
+ ok = cholb.Factorize(test.b)
+ if !ok {
+ t.Fatal("unexpected Cholesky factorization failure for b: not positive definite")
+ }
+
+ var x Dense
+ x.solveTwoChol(&chola, &cholb)
+
+ var ans Dense
+ ans.Mul(test.a, &x)
+ if !EqualApprox(&ans, test.b, 1e-12) {
+ var y Dense
+ y.Solve(test.a, test.b)
+ t.Errorf("incorrect Cholesky solve solution product\ngot solution:\n%.4v\nwant solution\n%.4v",
+ Formatted(&x), Formatted(&y))
+ }
+ }
+}
+
func TestCholeskySolveVec(t *testing.T) {
for _, test := range []struct {
a *SymDense

mat64: add higher order generalised singular value decomposition