mat64: add outer product method

This is subtly different enough from RankOne to justify its own method.
The equivalent code using RankOne would require a zero of m (which we
don't provide in an efficient way) or the use of a temporary.

mat64/dense_arithmetic.go

@@ -896,7 +896,7 @@ func (m *Dense) EqualsApprox(b Matrix, epsilon float64) bool {
 }
 
 // RankOne performs a rank-one update to the matrix a and stores the result
-// in the receiver
+// in the receiver. If a is zero, see Outer.
 //  m = a + alpha * x * y'
 func (m *Dense) RankOne(a Matrix, alpha float64, x, y *Vector) {
 	ar, ac := a.Dims()
@@ -920,3 +920,39 @@ func (m *Dense) RankOne(a Matrix, alpha float64, x, y *Vector) {
 	blas64.Ger(alpha, x.mat, y.mat, w.mat)
 	*m = w
 }
+
+// Outer calculates the outer product of x and y, and stores the result
+// in the receiver. In order to update to an existing matrix, see RankOne.
+//  m = x * y'
+func (m *Dense) Outer(x, y *Vector) {
+	r := x.Len()
+	c := y.Len()
+
+	// Copied from reuseAs with use replaced by useZeroed
+	// and a final zero of the matrix elements if we pass
+	// the shape checks.
+	// TODO(kortschak): Factor out into reuseZeroedAs if
+	// we find another case that needs it.
+	if m.mat.Rows > m.capRows || m.mat.Cols > m.capCols {
+		// Panic as a string, not a mat64.Error.
+		panic("mat64: caps not correctly set")
+	}
+	if m.isZero() {
+		m.mat = blas64.General{
+			Rows:   r,
+			Cols:   c,
+			Stride: c,
+			Data:   useZeroed(m.mat.Data, r*c),
+		}
+		m.capRows = r
+		m.capCols = c
+	} else if r != m.mat.Rows || c != m.mat.Cols {
+		panic(ErrShape)
+	} else {
+		for i := 0; i < r; i++ {
+			zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+c])
+		}
+	}
+
+	blas64.Ger(1, x.mat, y.mat, m.mat)
+}

mat64/dense_test.go

@@ -1342,6 +1342,53 @@ func (s *S) TestRankOne(c *check.C) {
 	}
 }
 
+func (s *S) TestOuter(c *check.C) {
+	for i, test := range []struct {
+		x []float64
+		y []float64
+	}{
+		{
+			x: []float64{5},
+			y: []float64{10},
+		},
+		{
+			x: []float64{5, 6, 1},
+			y: []float64{10},
+		},
+
+		{
+			x: []float64{5},
+			y: []float64{10, 15, 8},
+		},
+		{
+			x: []float64{1, 5},
+			y: []float64{10, 15},
+		},
+		{
+			x: []float64{2, 3, 9},
+			y: []float64{8, 9},
+		},
+		{
+			x: []float64{2, 3},
+			y: []float64{8, 9, 9},
+		},
+	} {
+		want := &Dense{}
+		xm := NewDense(len(test.x), 1, test.x)
+		ym := NewDense(1, len(test.y), test.y)
+
+		want.Mul(xm, ym)
+
+		var m Dense
+		// Check with a new matrix
+		m.Outer(NewVector(len(test.x), test.x), NewVector(len(test.y), test.y))
+		c.Check(m.Equals(want), check.Equals, true, check.Commentf("Test %v. Want %v, Got %v", i, want, m))
+		// Check with the same matrix
+		m.Outer(NewVector(len(test.x), test.x), NewVector(len(test.y), test.y))
+		c.Check(m.Equals(want), check.Equals, true, check.Commentf("Test %v. Want %v, Got %v", i, want, m))
+	}
+}
+
 var (
 	wd *Dense
 )