mat: Add methods for Cholesky to implement Matrix and Symmetric
* mat: Add methods for Cholesky to implement Matrix and Symmetric
Fixes
MODIFY
mat/cholesky.go
mat/cholesky.go
@@ -17,8 +17,21 @@ const (
badCholesky = "mat: invalid Cholesky factorization"
)
-// Cholesky is a type for creating and using the Cholesky factorization of a
-// symmetric positive definite matrix.
+var (
+ _ Matrix = (*Cholesky)(nil)
+ _ Symmetric = (*Cholesky)(nil)
+)
+
+// Cholesky is a symmetric positive definite matrix represented by its
+// Cholesky decomposition.
+//
+// The decomposition can be constructed using the Factorize method. The
+// factorization itself can be extracted using the UTo or LTo methods, and the
+// original symmetric matrix can be recovered with ToSym.
+//
+// Note that this matrix representation is useful for certain operations, in
+// particular finding solutions to linear equations. It is very inefficient
+// at other operations, in particular At is slow.
//
// Cholesky methods may only be called on a value that has been successfully
// initialized by a call to Factorize that has returned true. Calls to methods
@@ -58,6 +71,47 @@ func (c *Cholesky) updateCond(norm float64) {
c.cond = 1 / v
}
+// Dims returns the dimensions of the matrix.
+func (ch *Cholesky) Dims() (r, c int) {
+ if !ch.valid() {
+ panic(badCholesky)
+ }
+ r, c = ch.chol.Dims()
+ return r, c
+}
+
+// At returns the element at row i, column j.
+func (c *Cholesky) At(i, j int) float64 {
+ if !c.valid() {
+ panic(badCholesky)
+ }
+ n := c.Symmetric()
+ if uint(i) >= uint(n) {
+ panic(ErrRowAccess)
+ }
+ if uint(j) >= uint(n) {
+ panic(ErrColAccess)
+ }
+
+ var val float64
+ for k := 0; k <= min(i, j); k++ {
+ val += c.chol.at(k, i) * c.chol.at(k, j)
+ }
+ return val
+}
+
+// T returns the the receiver, the transpose of a symmetric matrix.
+func (c *Cholesky) T() Matrix {
+ return c
+}
+
+// Symmetric implements the Symmetric interface and returns the number of rows
+// in the matrix (this is also the number of columns).
+func (c *Cholesky) Symmetric() int {
+ r, _ := c.chol.Dims()
+ return r
+}
+
// Cond returns the condition number of the factorized matrix.
func (c *Cholesky) Cond() float64 {
if !c.valid() {
@@ -124,7 +178,7 @@ func (c *Cholesky) Clone(chol *Cholesky) {
if !chol.valid() {
panic(badCholesky)
}
- n := chol.Size()
+ n := chol.Symmetric()
if c.chol == nil {
c.chol = NewTriDense(n, Upper, nil)
} else {
@@ -134,14 +188,6 @@ func (c *Cholesky) Clone(chol *Cholesky) {
c.cond = chol.cond
}
-// Size returns the dimension of the factorized matrix.
-func (c *Cholesky) Size() int {
- if !c.valid() {
- panic(badCholesky)
- }
- return c.chol.mat.N
-}
-
// Det returns the determinant of the matrix that has been factorized.
func (c *Cholesky) Det() float64 {
if !c.valid() {
@@ -365,7 +411,7 @@ func (c *Cholesky) InverseTo(s *SymDense) error {
// the updated factorization is
// U'^T * U' = f A = A'
// Scale panics if the constant is non-positive, or if the receiver is non-zero
-// and is of a different Size from the input.
+// and is of a different size from the input.
func (c *Cholesky) Scale(f float64, orig *Cholesky) {
if !orig.valid() {
panic(badCholesky)
@@ -373,7 +419,7 @@ func (c *Cholesky) Scale(f float64, orig *Cholesky) {
if f <= 0 {
panic("cholesky: scaling by a non-positive constant")
}
- n := orig.Size()
+ n := orig.Symmetric()
if c.chol == nil {
c.chol = NewTriDense(n, Upper, nil)
} else if c.chol.mat.N != n {
@@ -392,10 +438,11 @@ func (c *Cholesky) Scale(f float64, orig *Cholesky) {
// that k > w' A^-1 w. If this condition does not hold then ExtendVecSym will
// return false and the receiver will not be updated.
//
-// ExtendVecSym will panic if v.Len() != a.Size()+1 or if a does not contain
+// ExtendVecSym will panic if v.Len() != a.Symmetric()+1 or if a does not contain
// a valid decomposition.
func (c *Cholesky) ExtendVecSym(a *Cholesky, v Vector) (ok bool) {
- n := a.Size()
+ n := a.Symmetric()
+
if v.Len() != n+1 {
panic(badSliceLength)
}
@@ -465,7 +512,7 @@ func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool)
if !orig.valid() {
panic(badCholesky)
}
- n := orig.Size()
+ n := orig.Symmetric()
if r, c := x.Dims(); r != n || c != 1 {
panic(ErrShape)
}

MODIFY
mat/cholesky_test.go
mat/cholesky_test.go
@@ -72,6 +72,36 @@ func TestCholesky(t *testing.T) {
}
}
+func TestCholeskyAt(t *testing.T) {
+ for _, test := range []*SymDense{
+ NewSymDense(3, []float64{
+ 53, 59, 37,
+ 59, 83, 71,
+ 37, 71, 101,
+ }),
+ } {
+ var chol Cholesky
+ ok := chol.Factorize(test)
+ if !ok {
+ t.Fatalf("Matrix not positive definite")
+ }
+ n := test.Symmetric()
+ cn := chol.Symmetric()
+ if cn != n {
+ t.Errorf("Cholesky size does not match. Got %d, want %d", cn, n)
+ }
+ for i := 0; i < n; i++ {
+ for j := 0; j < n; j++ {
+ got := chol.At(i, j)
+ want := test.At(i, j)
+ if math.Abs(got-want) > 1e-12 {
+ t.Errorf("Cholesky at does not match at %d, %d. Got %v, want %v", i, j, got, want)
+ }
+ }
+ }
+ }
+}
+
func TestCholeskySolveTo(t *testing.T) {
for _, test := range []struct {
a *SymDense

MODIFY
mat/symmetric.go
mat/symmetric.go
@@ -94,12 +94,13 @@ func (s *SymDense) Caps() (r, c int) {
return s.cap, s.cap
}
-// T implements the Matrix interface. Symmetric matrices, by definition, are
-// equal to their transpose, and this is a no-op.
+// T returns the receiver, the transpose of a symmetric matrix.
func (s *SymDense) T() Matrix {
return s
}
+// Symmetric implements the Symmetric interface and returns the number of rows
+// and columns in the matrix.
func (s *SymDense) Symmetric() int {
return s.mat.N
}

MODIFY
optimize/cmaes.go
optimize/cmaes.go
@@ -204,7 +204,7 @@ func (cma *CmaEsChol) Init(dim, tasks int) int {
cma.mean = resize(cma.mean, dim) // mean location initialized at the start of Run
if cma.InitCholesky != nil {
- if cma.InitCholesky.Size() != dim {
+ if cma.InitCholesky.Symmetric() != dim {
panic("cma-es-chol: incorrect InitCholesky size")
}
cma.chol.Clone(cma.InitCholesky)

MODIFY
stat/distmat/wishart.go
stat/distmat/wishart.go
@@ -113,7 +113,7 @@ func (w *Wishart) LogProbSym(x mat.Symmetric) float64 {
// LogProbSymChol returns the log of the probability of the input symmetric matrix
// given its Cholesky decomposition.
func (w *Wishart) LogProbSymChol(cholX *mat.Cholesky) float64 {
- dim := cholX.Size()
+ dim := cholX.Symmetric()
if dim != w.dim {
panic(badDim)
}

MODIFY
stat/distmv/normal.go
stat/distmv/normal.go
@@ -71,7 +71,7 @@ func NewNormal(mu []float64, sigma mat.Symmetric, src rand.Source) (*Normal, boo
// panics if len(mu) is not equal to chol.Size().
func NewNormalChol(mu []float64, chol *mat.Cholesky, src rand.Source) *Normal {
dim := len(mu)
- if dim != chol.Size() {
+ if dim != chol.Symmetric() {
panic(badSizeMismatch)
}
n := &Normal{
@@ -198,7 +198,7 @@ func NormalLogProb(x, mu []float64, chol *mat.Cholesky) float64 {
if len(x) != dim {
panic(badSizeMismatch)
}
- if chol.Size() != dim {
+ if chol.Symmetric() != dim {
panic(badSizeMismatch)
}
logSqrtDet := 0.5 * chol.LogDet()
@@ -302,7 +302,7 @@ func (n *Normal) Rand(x []float64) []float64 {
// available. Otherwise, the NewNormal function should be used.
func NormalRand(x, mean []float64, chol *mat.Cholesky, src rand.Source) []float64 {
x = reuseAs(x, len(mean))
- if len(mean) != chol.Size() {
+ if len(mean) != chol.Symmetric() {
panic(badInputLength)
}
if src == nil {

mat: Add methods for Cholesky to implement Matrix and Symmetric (#928)

* mat: Add methods for Cholesky to implement Matrix and Symmetric

Fixes #919.