diff --git a/matrix/matrix.go b/matrix/matrix.go index c94768b..5fa857b 100644 --- a/matrix/matrix.go +++ b/matrix/matrix.go @@ -3,6 +3,7 @@ package matrix import ( "bytes" "errors" + "fmt" "math" ) @@ -14,6 +15,9 @@ const ( var ( // ErrDimensionMismatch is a typical error. ErrDimensionMismatch = errors.New("dimension mismatch") + + // ErrSingularValue is a typical error. + ErrSingularValue = errors.New("singular value") ) // New returns a new matrix. @@ -220,6 +224,64 @@ func (m *Matrix) Row(row int) Vector { return Vector(values) } +// SubMatrix returns a sub matrix from a given outer matrix. +func (m *Matrix) SubMatrix(i, j, rows, cols int) *Matrix { + return &Matrix{ + elements: m.elements[i*m.stride+j : i*m.stride+j+(rows-1)*m.stride+cols], + stride: m.stride, + epsilon: m.epsilon, + } +} + +// ScaleRow applies a scale to an entire row. +func (m *Matrix) ScaleRow(row int, scale float64) { + startIndex := row * m.stride + for i := startIndex; i < m.stride; i++ { + m.elements[i] = m.elements[i] * scale + } +} + +func (m *Matrix) scaleAddRow(rd int, rs int, f float64) { + indexd := rd * m.stride + indexs := rs * m.stride + for col := 0; col < m.stride; col++ { + m.elements[indexd] += f * m.elements[indexs] + indexd++ + indexs++ + } +} + +// SwapRows swaps a row in the matrix in place. +func (m *Matrix) SwapRows(i, j int) { + var vi, vj float64 + for col := 0; col < m.stride; col++ { + vi = m.Get(i, col) + vj = m.Get(j, col) + m.Set(i, col, vj) + m.Set(j, col, vi) + } +} + +// Augment concatenates two matrices about the horizontal. +func (m *Matrix) Augment(m2 *Matrix) (*Matrix, error) { + mr, mc := m.Size() + m2r, m2c := m2.Size() + if mr != m2r { + return nil, ErrDimensionMismatch + } + + m3 := Zero(mr, mc+m2c) + for row := 0; row < mr; row++ { + for col := 0; col < mc; col++ { + m3.Set(row, col, m.Get(row, col)) + } + for col := 0; col < m2c; col++ { + m3.Set(row, mc+col, m2.Get(row, col)) + } + } + return m3, nil +} + // Copy returns a duplicate of a given matrix. func (m *Matrix) Copy() *Matrix { m2 := &Matrix{stride: m.stride, epsilon: m.epsilon, elements: make([]float64, len(m.elements))} @@ -354,6 +416,30 @@ func (m *Matrix) Pivotize() *Matrix { return p } +// Times returns the product of a matrix and another. +func (m *Matrix) Times(m2 *Matrix) (*Matrix, error) { + mr, mc := m.Size() + m2r, m2c := m2.Size() + + if mc != m2r { + return nil, fmt.Errorf("cannot multiply (%dx%d) and (%dx%d)", mr, mc, m2r, m2c) + //return nil, ErrDimensionMismatch + } + + c := Zero(mr, m2c) + + for i := 0; i < mr; i++ { + sums := c.elements[i*c.stride : (i+1)*c.stride] + for k, a := range m.elements[i*m.stride : i*m.stride+m.stride] { + for j, b := range m2.elements[k*m2.stride : k*m2.stride+m2.stride] { + sums[j] += a * b + } + } + } + + return c, nil +} + // Decompositions // LU performs the LU decomposition. @@ -459,3 +545,48 @@ func (m *Matrix) QR() (q, r *Matrix) { return } + +// Transpose flips a matrix about its diagonal, returning a new copy. +func (m *Matrix) Transpose() *Matrix { + rows, cols := m.Size() + m2 := Zero(cols, rows) + for i := 0; i < rows; i++ { + for j := 0; j < cols; j++ { + m2.Set(j, i, m.Get(i, j)) + } + } + return m2 +} + +// Inverse returns a matrix such that M*I==1. +func (m *Matrix) Inverse() (*Matrix, error) { + if !m.IsSymmetric() { + return nil, ErrDimensionMismatch + } + + rows, cols := m.Size() + + aug, _ := m.Augment(Eye(rows)) + for i := 0; i < rows; i++ { + j := i + for k := i; k < rows; k++ { + if math.Abs(aug.Get(k, i)) > math.Abs(aug.Get(j, i)) { + j = k + } + } + if j != i { + aug.SwapRows(i, j) + } + if aug.Get(i, i) == 0 { + return nil, ErrSingularValue + } + aug.ScaleRow(i, 1.0/aug.Get(i, i)) + for k := 0; k < rows; k++ { + if k == i { + continue + } + aug.scaleAddRow(k, i, -aug.Get(k, i)) + } + } + return aug.SubMatrix(0, cols, rows, cols), nil +} diff --git a/matrix/regression.go b/matrix/regression.go new file mode 100644 index 0000000..7aaea7c --- /dev/null +++ b/matrix/regression.go @@ -0,0 +1,45 @@ +package matrix + +import "errors" + +var ( + // ErrPolyRegArraysSameLength is a common error. + ErrPolyRegArraysSameLength = errors.New("polynomial array inputs must be the same length") +) + +// Poly returns the polynomial regress of a given degree over the given values. +func Poly(xvalues, yvalues []float64, degree int) ([]float64, error) { + if len(xvalues) != len(yvalues) { + return nil, ErrPolyRegArraysSameLength + } + + m := len(yvalues) + n := degree + 1 + y := New(m, 1, yvalues...) + x := Zero(m, n) + + for i := 0; i < m; i++ { + ip := float64(1) + for j := 0; j < n; j++ { + x.Set(i, j, ip) + ip *= xvalues[i] + } + } + + q, r := x.QR() + qty, err := q.Transpose().Times(y) + if err != nil { + return nil, err + } + + c := make([]float64, n) + for i := n - 1; i >= 0; i-- { + c[i] = qty.Get(i, 0) + for j := i + 1; j < n; j++ { + c[i] -= c[j] * r.Get(i, j) + } + c[i] /= r.Get(i, i) + } + + return c, nil +} diff --git a/matrix/regression_test.go b/matrix/regression_test.go new file mode 100644 index 0000000..c55a480 --- /dev/null +++ b/matrix/regression_test.go @@ -0,0 +1,22 @@ +package matrix + +import ( + "testing" + + assert "github.com/blendlabs/go-assert" +) + +func TestPoly(t *testing.T) { + assert := assert.New(t) + var xGiven = []float64{0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10} + var yGiven = []float64{1, 6, 17, 34, 57, 86, 121, 162, 209, 262, 321} + var degree = 2 + + c, err := Poly(xGiven, yGiven, degree) + assert.Nil(err) + assert.Len(c, 3) + + assert.InDelta(c[0], 0.999999999, DefaultEpsilon) + assert.InDelta(c[1], 2, DefaultEpsilon) + assert.InDelta(c[2], 3, DefaultEpsilon) +}