Skip to content
This repository has been archived by the owner on Dec 22, 2018. It is now read-only.

Commit

Permalink
Store sigma in Normal and StudentsT
Browse files Browse the repository at this point in the history 
Currently, we throw sigma away, and recompute it if necessary. This PR keeps sigma. This fixes an issue with concurrent calling of methods. In addition, however, it removes any possible issues with reconstructing a badly-conditioned sigma from its Cholesky decomposition, and avoids an extra n^3 work if sigma does need to be recomputed. The complexity of the implementation and difficulties listed above is not worth the memory savings in some cases, especially since the memory of the type is already O(n^2)
  • Loading branch information
@btracey
btracey committed Mar 7, 2017
1 parent cfbee06 commit 6d596f2
Show file tree
Hide file tree
Showing 5 changed files with 17 additions and 98 deletions.
36 changes: 7 additions & 29 deletions distmv/normal.go
View file
Expand Up @@ -7,7 +7,6 @@ package distmv
import (
"math"
"math/rand"
"sync"

"github.com/gonum/floats"
"github.com/gonum/matrix/mat64"
Expand All @@ -27,8 +26,7 @@ var (
type Normal struct {
mu []float64

once sync.Once
sigma *mat64.SymDense // only stored if needed
sigma mat64.SymDense

chol mat64.Cholesky
lower mat64.TriDense
Expand Down Expand Up @@ -59,6 +57,8 @@ func NewNormal(mu []float64, sigma mat64.Symmetric, src *rand.Rand) (*Normal, bo
if !ok {
return nil, false
}
n.sigma = *mat64.NewSymDense(dim, nil)
n.sigma.CopySym(sigma)
n.lower.LFromCholesky(&n.chol)
n.logSqrtDet = 0.5 * n.chol.LogDet()
return n, true
Expand Down Expand Up @@ -141,9 +141,7 @@ func (n *Normal) ConditionNormal(observed []int, values []float64, src *rand.Ran
}
}

n.setSigma()

_, mu1, sigma11 := studentsTConditional(observed, values, math.Inf(1), n.mu, n.sigma)
_, mu1, sigma11 := studentsTConditional(observed, values, math.Inf(1), n.mu, &n.sigma)
if mu1 == nil {
return nil, false
}
Expand All @@ -164,8 +162,7 @@ func (n *Normal) CovarianceMatrix(s *mat64.SymDense) *mat64.SymDense {
if sn != n.Dim() {
panic("normal: input matrix size mismatch")
}
n.setSigma()
s.CopySym(n.sigma)
s.CopySym(&n.sigma)
return s
}

Expand Down Expand Up @@ -202,9 +199,8 @@ func (n *Normal) MarginalNormal(vars []int, src *rand.Rand) (*Normal, bool) {
for i, v := range vars {
newMean[i] = n.mu[v]
}
n.setSigma()
var s mat64.SymDense
s.SubsetSym(n.sigma, vars)
s.SubsetSym(&n.sigma, vars)
return NewNormal(newMean, &s, src)
}

Expand All @@ -216,19 +212,9 @@ func (n *Normal) MarginalNormal(vars []int, src *rand.Rand) (*Normal, bool) {
//
// The input src is passed to the constructed distuv.Normal.
func (n *Normal) MarginalNormalSingle(i int, src *rand.Rand) distuv.Normal {
var std float64
if n.sigma != nil {
std = n.sigma.At(i, i)
} else {
// Reconstruct the {i,i} diagonal element of the covariance directly.
for j := 0; j <= i; j++ {
v := n.lower.At(i, j)
std += v * v
}
}
return distuv.Normal{
Mu: n.mu[i],
Sigma: math.Sqrt(std),
Sigma: math.Sqrt(n.sigma.At(i, i)),
Source: src,
}
}
Expand Down Expand Up @@ -300,14 +286,6 @@ func (n *Normal) SetMean(mu []float64) {
copy(n.mu, mu)
}

// setSigma computes and stores the covariance matrix of the distribution.
func (n *Normal) setSigma() {
n.once.Do(func() {
n.sigma = mat64.NewSymDense(n.Dim(), nil)
n.sigma.FromCholesky(&n.chol)
})
}

// TransformNormal transforms the vector, normal, generated from a standard
// multidimensional normal into a vector that has been generated under the
// distribution of the receiver.
Expand Down
15 changes: 0 additions & 15 deletions distmv/normal_test.go
View file
Expand Up @@ -489,8 +489,6 @@ func TestMarginalSingle(t *testing.T) {
if !ok {
t.Fatalf("Bad test, covariance matrix not positive definite")
}
// Verify with nil Sigma.
normal.sigma = nil
for i, mean := range test.mu {
norm := normal.MarginalNormalSingle(i, nil)
if norm.Mean() != mean {
Expand All @@ -501,19 +499,6 @@ func TestMarginalSingle(t *testing.T) {
t.Errorf("StdDev mismatch nil Sigma, idx %v: want %v, got %v.", i, std, norm.StdDev())
}
}

// Verify with non-nil Sigma.
normal.setSigma()
for i, mean := range test.mu {
norm := normal.MarginalNormalSingle(i, nil)
if norm.Mean() != mean {
t.Errorf("Mean mismatch non-nil Sigma, idx %v: want %v, got %v.", i, mean, norm.Mean())
}
std := math.Sqrt(test.sigma.At(i, i))
if math.Abs(norm.StdDev()-std) > 1e-14 {
t.Errorf("StdDev mismatch non-nil Sigma, idx %v: want %v, got %v.", i, std, norm.StdDev())
}
}
}

// Test matching with TestMarginal.
Expand Down
3 changes: 0 additions & 3 deletions distmv/normalbench_test.go
View file
Expand Up @@ -7,7 +7,6 @@ package distmv
import (
"log"
"math/rand"
"sync"
"testing"

"github.com/gonum/matrix/mat64"
Expand Down Expand Up @@ -38,8 +37,6 @@ func BenchmarkMarginalNormalReset10(b *testing.B) {
if !ok {
b.Error("bad test")
}
normal.sigma = nil
normal.once = sync.Once{}
_ = marg
}
}
Expand Down
37 changes: 7 additions & 30 deletions distmv/studentst.go
View file
Expand Up @@ -8,7 +8,6 @@ import (
"math"
"math/rand"
"sort"
"sync"

"golang.org/x/tools/container/intsets"

Expand Down Expand Up @@ -36,8 +35,7 @@ type StudentsT struct {
mu []float64
src *rand.Rand

once sync.Once
sigma *mat64.SymDense // only stored if needed
sigma mat64.SymDense // only stored if needed

chol mat64.Cholesky
lower mat64.TriDense
Expand Down Expand Up @@ -71,6 +69,8 @@ func NewStudentsT(mu []float64, sigma mat64.Symmetric, nu float64, src *rand.Ran
if !ok {
return nil, false
}
s.sigma = *mat64.NewSymDense(dim, nil)
s.sigma.CopySym(sigma)
s.lower.LFromCholesky(&s.chol)
s.logSqrtDet = 0.5 * s.chol.LogDet()
return s, true
Expand Down Expand Up @@ -101,9 +101,7 @@ func (s *StudentsT) ConditionStudentsT(observed []int, values []float64, src *ra
}
}

s.setSigma()

newNu, newMean, newSigma := studentsTConditional(observed, values, s.nu, s.mu, s.sigma)
newNu, newMean, newSigma := studentsTConditional(observed, values, s.nu, s.mu, &s.sigma)
if newMean == nil {
return nil, false
}
Expand Down Expand Up @@ -231,8 +229,7 @@ func (st *StudentsT) CovarianceMatrix(s *mat64.SymDense) *mat64.SymDense {
if sn != st.dim {
panic("normal: input matrix size mismatch")
}
st.setSigma()
s.CopySym(st.sigma)
s.CopySym(&st.sigma)
s.ScaleSym(st.nu/(st.nu-2), s)
return s
}
Expand Down Expand Up @@ -286,9 +283,8 @@ func (s *StudentsT) MarginalStudentsT(vars []int, src *rand.Rand) (dist *Student
for i, v := range vars {
newMean[i] = s.mu[v]
}
s.setSigma()
var newSigma mat64.SymDense
newSigma.SubsetSym(s.sigma, vars)
newSigma.SubsetSym(&s.sigma, vars)
return NewStudentsT(newMean, &newSigma, s.nu, src)
}

Expand All @@ -300,20 +296,9 @@ func (s *StudentsT) MarginalStudentsT(vars []int, src *rand.Rand) (dist *Student
//
// The input src is passed to the call to NewStudentsT.
func (s *StudentsT) MarginalStudentsTSingle(i int, src *rand.Rand) distuv.StudentsT {
var std float64
if s.sigma != nil {
std = s.sigma.At(i, i)
} else {
// Reconstruct the {i,i} diagonal element of the covariance directly.
for j := 0; j <= i; j++ {
v := s.lower.At(i, j)
std += v * v
}
}

return distuv.StudentsT{
Mu: s.mu[i],
Sigma: math.Sqrt(std),
Sigma: math.Sqrt(s.sigma.At(i, i)),
Nu: s.nu,
Src: src,
}
Expand Down Expand Up @@ -367,11 +352,3 @@ func (s *StudentsT) Rand(x []float64) []float64 {
floats.Add(x, s.mu)
return x
}

// setSigma computes and stores the covariance matrix of the distribution.
func (s *StudentsT) setSigma() {
s.once.Do(func() {
s.sigma = mat64.NewSymDense(s.dim, nil)
s.sigma.FromCholesky(&s.chol)
})
}
24 changes: 3 additions & 21 deletions distmv/studentst_test.go
View file
Expand Up @@ -182,12 +182,9 @@ func TestStudentsTConditional(t *testing.T) {
muOb[i] = test.mean[v]
}

s.setSigma()
sUp.setSigma()

var sig11, sig22 mat64.SymDense
sig11.SubsetSym(s.sigma, unob)
sig22.SubsetSym(s.sigma, ob)
sig11.SubsetSym(&s.sigma, unob)
sig22.SubsetSym(&s.sigma, ob)

sig12 := mat64.NewDense(len(unob), len(ob), nil)
for i := range unob {
Expand Down Expand Up @@ -221,7 +218,7 @@ func TestStudentsTConditional(t *testing.T) {

dot := mat64.Dot(shiftVec, &tmp)
tmp3.Scale((test.nu+dot)/(test.nu+float64(len(ob))), &tmp3)
if !mat64.EqualApprox(&tmp3, sUp.sigma, 1e-10) {
if !mat64.EqualApprox(&tmp3, &sUp.sigma, 1e-10) {
t.Errorf("Sigma mismatch")
}
}
Expand All @@ -248,8 +245,6 @@ func TestStudentsTMarginalSingle(t *testing.T) {
if !ok {
t.Fatalf("Bad test, covariance matrix not positive definite")
}
// Verify with nil Sigma.
studentst.sigma = nil
for i, mean := range test.mu {
st := studentst.MarginalStudentsTSingle(i, nil)
if st.Mean() != mean {
Expand All @@ -263,18 +258,5 @@ func TestStudentsTMarginalSingle(t *testing.T) {
t.Errorf("Nu mismatch nil Sigma, idx %v: want %v, got %v ", i, test.nu, st.Nu)
}
}

// Verify with non-nil Sigma.
studentst.setSigma()
for i, mean := range test.mu {
st := studentst.MarginalStudentsTSingle(i, nil)
if st.Mean() != mean {
t.Errorf("Mean mismatch non-nil Sigma, idx %v: want %v, got %v.", i, mean, st.Mean())
}
std := math.Sqrt(test.sigma.At(i, i))
if math.Abs(st.Sigma-std) > 1e-14 {
t.Errorf("StdDev mismatch non-nil Sigma, idx %v: want %v, got %v.", i, std, st.StdDev())
}
}
}
}

0 comments on commit 6d596f2

Please sign in to comment.