stat,stat/distmv: don't allocate matrix on nil input
MODIFY
stat/distmv/dirichlet.go
stat/distmv/dirichlet.go
@@ -56,32 +56,30 @@ func NewDirichlet(alpha []float64, src rand.Source) *Dirichlet {
return d
}
-// CovarianceMatrix returns the covariance matrix of the distribution. Upon
-// return, the value at element {i, j} of the covariance matrix is equal to
-// the covariance of the i^th and j^th variables.
+// CovarianceMatrix calculates the covariance matrix of the distribution,
+// storing the result in dst. Upon return, the value at element {i, j} of the
+// covariance matrix is equal to the covariance of the i^th and j^th variables.
// covariance(i, j) = E[(x_i - E[x_i])(x_j - E[x_j])]
-// If the input matrix is nil a new matrix is allocated, otherwise the result
-// is stored in-place into the input.
-func (d *Dirichlet) CovarianceMatrix(cov *mat.SymDense) *mat.SymDense {
- if cov == nil {
- cov = mat.NewSymDense(d.Dim(), nil)
- } else if cov.Symmetric() == 0 {
- *cov = *(cov.GrowSym(d.dim).(*mat.SymDense))
- } else if cov.Symmetric() != d.dim {
- panic("normal: input matrix size mismatch")
+// If the dst matrix is zero-sized it will be resized to the correct dimensions,
+// otherwise dst must match the dimension of the receiver or CovarianceMatrix
+// will panic.
+func (d *Dirichlet) CovarianceMatrix(dst *mat.SymDense) {
+ if dst.IsZero() {
+ *dst = *(dst.GrowSym(d.dim).(*mat.SymDense))
+ } else if dst.Symmetric() != d.dim {
+ panic("dirichelet: input matrix size mismatch")
}
scale := 1 / (d.sumAlpha * d.sumAlpha * (d.sumAlpha + 1))
for i := 0; i < d.dim; i++ {
ai := d.alpha[i]
v := ai * (d.sumAlpha - ai) * scale
- cov.SetSym(i, i, v)
+ dst.SetSym(i, i, v)
for j := i + 1; j < d.dim; j++ {
aj := d.alpha[j]
v := -ai * aj * scale
- cov.SetSym(i, j, v)
+ dst.SetSym(i, j, v)
}
}
- return cov
}
// genLBeta computes the generalized LBeta function.

MODIFY
stat/distmv/general_test.go
stat/distmv/general_test.go
@@ -71,26 +71,28 @@ func checkMean(t *testing.T, cas int, x *mat.Dense, m Meaner, tol float64) {
}
type Cover interface {
- CovarianceMatrix(*mat.SymDense) *mat.SymDense
+ CovarianceMatrix(*mat.SymDense)
}
func checkCov(t *testing.T, cas int, x *mat.Dense, c Cover, tol float64) {
- cov := c.CovarianceMatrix(nil)
+ var cov mat.SymDense
+ c.CovarianceMatrix(&cov)
n := cov.Symmetric()
cov2 := mat.NewSymDense(n, nil)
c.CovarianceMatrix(cov2)
- if !mat.Equal(cov, cov2) {
+ if !mat.Equal(&cov, cov2) {
t.Errorf("Cov mismatch when providing nil and matrix. Case %v", cas)
}
var cov3 mat.SymDense
c.CovarianceMatrix(&cov3)
- if !mat.Equal(cov, &cov3) {
+ if !mat.Equal(&cov, &cov3) {
t.Errorf("Cov mismatch when providing zero matrix. Case %v", cas)
}
// Check that the covariance matrix matches the samples
- covEst := stat.CovarianceMatrix(nil, x, nil)
- if !mat.EqualApprox(covEst, cov, tol) {
- t.Errorf("Return cov and sample cov mismatch. Cas %v.\nGot:\n%0.4v\nWant:\n%0.4v", cas, mat.Formatted(cov), mat.Formatted(covEst))
+ var covEst mat.SymDense
+ stat.CovarianceMatrix(&covEst, x, nil)
+ if !mat.EqualApprox(&covEst, &cov, tol) {
+ t.Errorf("Return cov and sample cov mismatch. Cas %v.\nGot:\n%0.4v\nWant:\n%0.4v", cas, mat.Formatted(&cov), mat.Formatted(&covEst))
}
}

MODIFY
stat/distmv/normal.go
stat/distmv/normal.go
@@ -150,22 +150,20 @@ func (n *Normal) ConditionNormal(observed []int, values []float64, src rand.Sour
return NewNormal(mu1, sigma11, src)
}
-// CovarianceMatrix returns the covariance matrix of the distribution. Upon
-// return, the value at element {i, j} of the covariance matrix is equal to
-// the covariance of the i^th and j^th variables.
+// CovarianceMatrix stores the covariance matrix of the distribution in dst.
+// Upon return, the value at element {i, j} of the covariance matrix is equal
+// to the covariance of the i^th and j^th variables.
// covariance(i, j) = E[(x_i - E[x_i])(x_j - E[x_j])]
-// If the input matrix is nil a new matrix is allocated, otherwise the result
-// is stored in-place into the input.
-func (n *Normal) CovarianceMatrix(s *mat.SymDense) *mat.SymDense {
- if s == nil {
- s = mat.NewSymDense(n.Dim(), nil)
- }
- sn := s.Symmetric()
- if sn != n.Dim() {
+// If the dst matrix is zero-sized it will be resized to the correct dimensions,
+// otherwise dst must match the dimension of the receiver or CovarianceMatrix
+// will panic.
+func (n *Normal) CovarianceMatrix(dst *mat.SymDense) {
+ if dst.IsZero() {
+ *dst = *(dst.GrowSym(n.dim).(*mat.SymDense))
+ } else if dst.Symmetric() != n.dim {
panic("normal: input matrix size mismatch")
}
- s.CopySym(&n.sigma)
- return s
+ dst.CopySym(&n.sigma)
}
// Dim returns the dimension of the distribution.

MODIFY
stat/distmv/normal_test.go
stat/distmv/normal_test.go
@@ -121,9 +121,10 @@ func TestNormRand(t *testing.T) {
if !floats.EqualApprox(estMean, test.mean, 1e-2) {
t.Errorf("Mean mismatch: want: %v, got %v", test.mean, estMean)
}
- estCov := stat.CovarianceMatrix(nil, samps, nil)
- if !mat.EqualApprox(estCov, cov, 1e-2) {
- t.Errorf("Cov mismatch: want: %v, got %v", cov, estCov)
+ var estCov mat.SymDense
+ stat.CovarianceMatrix(&estCov, samps, nil)
+ if !mat.EqualApprox(&estCov, cov, 1e-2) {
+ t.Errorf("Cov mismatch: want: %v, got %v", cov, &estCov)
}
}
}
@@ -165,9 +166,10 @@ func TestNormalQuantile(t *testing.T) {
if !floats.EqualApprox(estMean, test.mean, 1e-2) {
t.Errorf("Mean mismatch: want: %v, got %v", test.mean, estMean)
}
- estCov := stat.CovarianceMatrix(nil, samps, nil)
- if !mat.EqualApprox(estCov, cov, 1e-2) {
- t.Errorf("Cov mismatch: want: %v, got %v", cov, estCov)
+ var estCov mat.SymDense
+ stat.CovarianceMatrix(&estCov, samps, nil)
+ if !mat.EqualApprox(&estCov, cov, 1e-2) {
+ t.Errorf("Cov mismatch: want: %v, got %v", cov, &estCov)
}
}
}
@@ -352,7 +354,8 @@ func TestConditionNormal(t *testing.T) {
for i := range estMean {
estMean[i] = stat.Mean(mat.Col(nil, i, samples), nil)
}
- estCov := stat.CovarianceMatrix(nil, samples, nil)
+ var estCov mat.SymDense
+ stat.CovarianceMatrix(&estCov, samples, nil)
// Compute update rule.
newNormal, ok := normal.ConditionNormal(test.observed, test.value, nil)
@@ -399,14 +402,15 @@ func TestCovarianceMatrix(t *testing.T) {
if !ok {
t.Fatalf("Bad test, covariance matrix not positive definite")
}
- cov := normal.CovarianceMatrix(nil)
- if !mat.EqualApprox(cov, test.sigma, 1e-14) {
+ var cov mat.SymDense
+ normal.CovarianceMatrix(&cov)
+ if !mat.EqualApprox(&cov, test.sigma, 1e-14) {
t.Errorf("Covariance mismatch with nil input")
}
dim := test.sigma.Symmetric()
- cov = mat.NewSymDense(dim, nil)
- normal.CovarianceMatrix(cov)
- if !mat.EqualApprox(cov, test.sigma, 1e-14) {
+ cov = *mat.NewSymDense(dim, nil)
+ normal.CovarianceMatrix(&cov)
+ if !mat.EqualApprox(&cov, test.sigma, 1e-14) {
t.Errorf("Covariance mismatch with supplied input")
}
}
@@ -459,8 +463,10 @@ func TestMarginal(t *testing.T) {
}
}
- marginalCov := marginal.CovarianceMatrix(nil)
- estCov := stat.CovarianceMatrix(nil, samps, nil)
+ var marginalCov mat.SymDense
+ marginal.CovarianceMatrix(&marginalCov)
+ var estCov mat.SymDense
+ stat.CovarianceMatrix(&estCov, samps, nil)
for i, v1 := range test.marginal {
for j, v2 := range test.marginal {
c := marginalCov.At(i, j)
@@ -532,7 +538,9 @@ func TestMarginalSingle(t *testing.T) {
if math.Abs(single.Mean()-mult.Mean(nil)[0]) > 1e-14 {
t.Errorf("Mean mismatch")
}
- if math.Abs(single.Variance()-mult.CovarianceMatrix(nil).At(0, 0)) > 1e-14 {
+ var cov mat.SymDense
+ mult.CovarianceMatrix(&cov)
+ if math.Abs(single.Variance()-cov.At(0, 0)) > 1e-14 {
t.Errorf("Variance mismatch")
}
}

MODIFY
stat/distmv/normalbench_test.go
stat/distmv/normalbench_test.go
@@ -17,7 +17,7 @@ func BenchmarkMarginalNormal10(b *testing.B) {
sz := 10
rnd := rand.New(rand.NewSource(1))
normal := randomNormal(sz, rnd)
- _ = normal.CovarianceMatrix(nil) // pre-compute sigma
+ normal.CovarianceMatrix(&mat.SymDense{}) // pre-compute sigma
b.ResetTimer()
for i := 0; i < b.N; i++ {
marg, ok := normal.MarginalNormal([]int{1}, nil)

MODIFY
stat/distmv/studentst.go
stat/distmv/studentst.go
@@ -221,23 +221,21 @@ func findUnob(observed []int, dim int) (unobserved []int) {
return unobserved
}
-// CovarianceMatrix returns the covariance matrix of the distribution. Upon
-// return, the value at element {i, j} of the covariance matrix is equal to
-// the covariance of the i^th and j^th variables.
+// CovarianceMatrix calculates the covariance matrix of the distribution,
+// storing the result in dst. Upon return, the value at element {i, j} of the
+// covariance matrix is equal to the covariance of the i^th and j^th variables.
// covariance(i, j) = E[(x_i - E[x_i])(x_j - E[x_j])]
-// If the input matrix is nil a new matrix is allocated, otherwise the result
-// is stored in-place into the input.
-func (st *StudentsT) CovarianceMatrix(s *mat.SymDense) *mat.SymDense {
- if s == nil {
- s = mat.NewSymDense(st.dim, nil)
- }
- sn := s.Symmetric()
- if sn != st.dim {
- panic("normal: input matrix size mismatch")
- }
- s.CopySym(&st.sigma)
- s.ScaleSym(st.nu/(st.nu-2), s)
- return s
+// If the dst matrix is zero-sized it will be resized to the correct dimensions,
+// otherwise dst must match the dimension of the receiver or CovarianceMatrix
+// will panic.
+func (st *StudentsT) CovarianceMatrix(dst *mat.SymDense) {
+ if dst.IsZero() {
+ *dst = *(dst.GrowSym(st.dim).(*mat.SymDense))
+ } else if dst.Symmetric() != st.dim {
+ panic("studentst: input matrix size mismatch")
+ }
+ dst.CopySym(&st.sigma)
+ dst.ScaleSym(st.nu/(st.nu-2), dst)
}
// Dim returns the dimension of the distribution.

MODIFY
stat/distmv/studentst_test.go
stat/distmv/studentst_test.go
@@ -129,10 +129,11 @@ func TestStudentsTRand(t *testing.T) {
if !floats.EqualApprox(estMean, mean, 1e-2) {
t.Errorf("Mean mismatch: want: %v, got %v", test.mean, estMean)
}
- cov := s.CovarianceMatrix(nil)
- estCov := stat.CovarianceMatrix(nil, samps, nil)
- if !mat.EqualApprox(estCov, cov, test.tolcov) {
- t.Errorf("Case %d: Cov mismatch: want: %v, got %v", cas, cov, estCov)
+ var cov, estCov mat.SymDense
+ s.CovarianceMatrix(&cov)
+ stat.CovarianceMatrix(&estCov, samps, nil)
+ if !mat.EqualApprox(&estCov, &cov, test.tolcov) {
+ t.Errorf("Case %d: Cov mismatch: want: %v, got %v", cas, &cov, &estCov)
}
}
}

MODIFY
stat/samplemv/sample_test.go
stat/samplemv/sample_test.go
@@ -182,7 +182,8 @@ func randomNormal(dim int, src *rand.Rand) (*distmv.Normal, bool) {
func compareNormal(t *testing.T, want *distmv.Normal, batch *mat.Dense, weights []float64, meanTol, covTol float64) {
dim := want.Dim()
mu := want.Mean(nil)
- sigma := want.CovarianceMatrix(nil)
+ var sigma mat.SymDense
+ want.CovarianceMatrix(&sigma)
n, _ := batch.Dims()
if weights == nil {
weights = make([]float64, n)
@@ -198,8 +199,9 @@ func compareNormal(t *testing.T, want *distmv.Normal, batch *mat.Dense, weights
}
}
- cov := stat.CovarianceMatrix(nil, batch, weights)
- if !mat.EqualApprox(cov, sigma, covTol) {
+ var cov mat.SymDense
+ stat.CovarianceMatrix(&cov, batch, weights)
+ if !mat.EqualApprox(&cov, &sigma, covTol) {
t.Errorf("Covariance matrix mismatch")
}
}

MODIFY
stat/statmat.go
stat/statmat.go
@@ -11,27 +11,25 @@ import (
"gonum.org/v1/gonum/mat"
)
-// CovarianceMatrix returns the covariance matrix (also known as the
+// CovarianceMatrix calculates the covariance matrix (also known as the
// variance-covariance matrix) calculated from a matrix of data, x, using
-// a two-pass algorithm.
+// a two-pass algorithm. The result is stored in dst.
//
// If weights is not nil the weighted covariance of x is calculated. weights
// must have length equal to the number of rows in input data matrix and
// must not contain negative elements.
-// If cov is not nil it must either be zero-sized or have the same number of
-// columns as the input data matrix. cov will be used as the destination for
-// the covariance data. If cov is nil, a new mat.SymDense is allocated for
-// the destination.
-func CovarianceMatrix(cov *mat.SymDense, x mat.Matrix, weights []float64) *mat.SymDense {
+// The dst matrix must either be zero-sized or have the same number of
+// columns as the input data matrix.
+func CovarianceMatrix(dst *mat.SymDense, x mat.Matrix, weights []float64) {
// This is the matrix version of the two-pass algorithm. It doesn't use the
// additional floating point error correction that the Covariance function uses
// to reduce the impact of rounding during centering.
r, c := x.Dims()
- if cov == nil {
- cov = mat.NewSymDense(c, nil)
- } else if n := cov.Symmetric(); n != c && n != 0 {
+ if dst.IsZero() {
+ *dst = *(dst.GrowSym(c).(*mat.SymDense))
+ } else if n := dst.Symmetric(); n != c {
panic(mat.ErrShape)
}
@@ -49,8 +47,8 @@ func CovarianceMatrix(cov *mat.SymDense, x mat.Matrix, weights []float64) *mat.S
if weights == nil {
// Calculate the normalization factor
// scaled by the sample size.
- cov.SymOuterK(1/(float64(r)-1), &xt)
- return cov
+ dst.SymOuterK(1/(float64(r)-1), &xt)
+ return
}
// Multiply by the sqrt of the weights, so that multiplication is symmetric.
@@ -69,25 +67,21 @@ func CovarianceMatrix(cov *mat.SymDense, x mat.Matrix, weights []float64) *mat.S
// Calculate the normalization factor
// scaled by the weighted sample size.
- cov.SymOuterK(1/(floats.Sum(weights)-1), &xt)
- return cov
+ dst.SymOuterK(1/(floats.Sum(weights)-1), &xt)
}
// CorrelationMatrix returns the correlation matrix calculated from a matrix
-// of data, x, using a two-pass algorithm.
+// of data, x, using a two-pass algorithm. The result is stored in dst.
//
// If weights is not nil the weighted correlation of x is calculated. weights
// must have length equal to the number of rows in input data matrix and
// must not contain negative elements.
-// If corr is not nil it must either be zero-sized or have the same number of
-// columns as the input data matrix. corr will be used as the destination for
-// the correlation data. If corr is nil, a new mat.SymDense is allocated for
-// the destination.
-func CorrelationMatrix(corr *mat.SymDense, x mat.Matrix, weights []float64) *mat.SymDense {
+// The dst matrix must either be zero-sized or have the same number of
+// columns as the input data matrix.
+func CorrelationMatrix(dst *mat.SymDense, x mat.Matrix, weights []float64) {
// This will panic if the sizes don't match, or if weights is the wrong size.
- corr = CovarianceMatrix(corr, x, weights)
- covToCorr(corr)
- return corr
+ CovarianceMatrix(dst, x, weights)
+ covToCorr(dst)
}
// covToCorr converts a covariance matrix to a correlation matrix.

MODIFY
stat/statmat_test.go
stat/statmat_test.go
@@ -61,10 +61,11 @@ func TestCovarianceMatrix(t *testing.T) {
if test.weights != nil {
copy(w, test.weights)
}
- for _, cov := range []*mat.SymDense{nil, {}} {
- c := CovarianceMatrix(cov, test.data, test.weights)
- if !mat.Equal(c, test.ans) {
- t.Errorf("%d: expected cov %v, found %v", i, test.ans, c)
+
+ var cov mat.SymDense
+ CovarianceMatrix(&cov, test.data, test.weights)
+ if !mat.Equal(&cov, test.ans) {
+ t.Errorf("%d: expected cov %v, found %v", i, test.ans, &cov)
}
if !floats.Equal(d, r.Data) {
t.Errorf("%d: data was modified during execution", i)
@@ -74,15 +75,14 @@ func TestCovarianceMatrix(t *testing.T) {
}
// compare with call to Covariance
- _, cols := c.Dims()
+ _, cols := cov.Dims()
for ci := 0; ci < cols; ci++ {
for cj := 0; cj < cols; cj++ {
x := mat.Col(nil, ci, test.data)
y := mat.Col(nil, cj, test.data)
- cov := Covariance(x, y, test.weights)
- if math.Abs(cov-c.At(ci, cj)) > 1e-14 {
- t.Errorf("CovMat does not match at (%v, %v). Want %v, got %v.", ci, cj, cov, c.At(ci, cj))
- }
+ wc := Covariance(x, y, test.weights)
+ if math.Abs(wc-cov.At(ci, cj)) > 1e-14 {
+ t.Errorf("CovMat does not match at (%v, %v). Want %v, got %v.", ci, cj, wc, cov.At(ci, cj))
}
}
}
@@ -157,10 +157,11 @@ func TestCorrelationMatrix(t *testing.T) {
if test.weights != nil {
copy(w, test.weights)
}
- for _, corr := range []*mat.SymDense{nil, {}} {
- c := CorrelationMatrix(corr, test.data, test.weights)
- if !mat.Equal(c, test.ans) {
- t.Errorf("%d: expected corr %v, found %v", i, test.ans, c)
+
+ var corr mat.SymDense
+ CorrelationMatrix(&corr, test.data, test.weights)
+ if !mat.Equal(&corr, test.ans) {
+ t.Errorf("%d: expected corr %v, found %v", i, test.ans, &corr)
}
if !floats.Equal(d, r.Data) {
t.Errorf("%d: data was modified during execution", i)
@@ -170,15 +171,14 @@ func TestCorrelationMatrix(t *testing.T) {
}
// compare with call to Covariance
- _, cols := c.Dims()
+ _, cols := corr.Dims()
for ci := 0; ci < cols; ci++ {
for cj := 0; cj < cols; cj++ {
x := mat.Col(nil, ci, test.data)
y := mat.Col(nil, cj, test.data)
- corr := Correlation(x, y, test.weights)
- if math.Abs(corr-c.At(ci, cj)) > 1e-14 {
- t.Errorf("CorrMat does not match at (%v, %v). Want %v, got %v.", ci, cj, corr, c.At(ci, cj))
- }
+ wc := Correlation(x, y, test.weights)
+ if math.Abs(wc-corr.At(ci, cj)) > 1e-14 {
+ t.Errorf("CorrMat does not match at (%v, %v). Want %v, got %v.", ci, cj, wc, corr.At(ci, cj))
}
}
}
@@ -231,8 +231,9 @@ func TestCorrCov(t *testing.T) {
},
},
} {
- corr := CorrelationMatrix(nil, test.data, test.weights)
- cov := CovarianceMatrix(nil, test.data, test.weights)
+ var corr, cov mat.SymDense
+ CorrelationMatrix(&corr, test.data, test.weights)
+ CovarianceMatrix(&cov, test.data, test.weights)
r := cov.Symmetric()
@@ -243,17 +244,17 @@ func TestCorrCov(t *testing.T) {
}
covFromCorr := mat.NewSymDense(corr.Symmetric(), nil)
- covFromCorr.CopySym(corr)
+ covFromCorr.CopySym(&corr)
corrToCov(covFromCorr, sigmas)
corrFromCov := mat.NewSymDense(cov.Symmetric(), nil)
- corrFromCov.CopySym(cov)
+ corrFromCov.CopySym(&cov)
covToCorr(corrFromCov)
- if !mat.EqualApprox(corr, corrFromCov, 1e-14) {
+ if !mat.EqualApprox(&corr, corrFromCov, 1e-14) {
t.Errorf("%d: corrToCov did not match direct Correlation calculation. Want: %v, got: %v. ", i, corr, corrFromCov)
}
- if !mat.EqualApprox(cov, covFromCorr, 1e-14) {
+ if !mat.EqualApprox(&cov, covFromCorr, 1e-14) {
t.Errorf("%d: covToCorr did not match direct Covariance calculation. Want: %v, got: %v. ", i, cov, covFromCorr)
}
@@ -434,12 +435,13 @@ func BenchmarkCovarianceMatrixHugexSmallInPlace(b *testing.B) {
func BenchmarkCovToCorr(b *testing.B) {
// generate a 10x10 covariance matrix
m := randMat(small, small)
- c := CovarianceMatrix(nil, m, nil)
- cc := mat.NewSymDense(c.Symmetric(), nil)
+ var cov mat.SymDense
+ CovarianceMatrix(&cov, m, nil)
+ cc := mat.NewSymDense(cov.Symmetric(), nil)
b.ResetTimer()
for i := 0; i < b.N; i++ {
b.StopTimer()
- cc.CopySym(c)
+ cc.CopySym(&cov)
b.StartTimer()
covToCorr(cc)
}
@@ -448,8 +450,9 @@ func BenchmarkCovToCorr(b *testing.B) {
func BenchmarkCorrToCov(b *testing.B) {
// generate a 10x10 correlation matrix
m := randMat(small, small)
- c := CorrelationMatrix(nil, m, nil)
- cc := mat.NewSymDense(c.Symmetric(), nil)
+ var corr mat.SymDense
+ CorrelationMatrix(&corr, m, nil)
+ cc := mat.NewSymDense(corr.Symmetric(), nil)
sigma := make([]float64, small)
for i := range sigma {
sigma[i] = 2
@@ -457,7 +460,7 @@ func BenchmarkCorrToCov(b *testing.B) {
b.ResetTimer()
for i := 0; i < b.N; i++ {
b.StopTimer()
- cc.CopySym(c)
+ cc.CopySym(&corr)
b.StartTimer()
corrToCov(cc, sigma)
}

stat,stat/distmv: don't allocate matrix on nil input