blob: eec0c07dd07edcc14770d2017ee2ccda3a40c54c [file] [log] [blame]
// Copyright ©2015 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package mat
import (
"math"
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
)
var (
symDense *SymDense
_ Matrix = symDense
_ allMatrix = symDense
_ denseMatrix = symDense
_ Symmetric = symDense
_ RawSymmetricer = symDense
_ MutableSymmetric = symDense
)
const badSymTriangle = "mat: blas64.Symmetric not upper"
// SymDense is a symmetric matrix that uses dense storage. SymDense
// matrices are stored in the upper triangle.
type SymDense struct {
mat blas64.Symmetric
cap int
}
// Symmetric represents a symmetric matrix (where the element at {i, j} equals
// the element at {j, i}). Symmetric matrices are always square.
type Symmetric interface {
Matrix
// Symmetric returns the number of rows/columns in the matrix.
Symmetric() int
}
// A RawSymmetricer can return a view of itself as a BLAS Symmetric matrix.
type RawSymmetricer interface {
RawSymmetric() blas64.Symmetric
}
// A MutableSymmetric can set elements of a symmetric matrix.
type MutableSymmetric interface {
Symmetric
SetSym(i, j int, v float64)
}
// NewSymDense creates a new Symmetric matrix with n rows and columns. If data == nil,
// a new slice is allocated for the backing slice. If len(data) == n*n, data is
// used as the backing slice, and changes to the elements of the returned SymDense
// will be reflected in data. If neither of these is true, NewSymDense will panic.
// NewSymDense will panic if n is zero.
//
// The data must be arranged in row-major order, i.e. the (i*c + j)-th
// element in the data slice is the {i, j}-th element in the matrix.
// Only the values in the upper triangular portion of the matrix are used.
func NewSymDense(n int, data []float64) *SymDense {
if n <= 0 {
if n == 0 {
panic(ErrZeroLength)
}
panic("mat: negative dimension")
}
if data != nil && n*n != len(data) {
panic(ErrShape)
}
if data == nil {
data = make([]float64, n*n)
}
return &SymDense{
mat: blas64.Symmetric{
N: n,
Stride: n,
Data: data,
Uplo: blas.Upper,
},
cap: n,
}
}
// Dims returns the number of rows and columns in the matrix.
func (s *SymDense) Dims() (r, c int) {
return s.mat.N, s.mat.N
}
// Caps returns the number of rows and columns in the backing matrix.
func (s *SymDense) Caps() (r, c int) {
return s.cap, s.cap
}
// 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
}
// RawSymmetric returns the matrix as a blas64.Symmetric. The returned
// value must be stored in upper triangular format.
func (s *SymDense) RawSymmetric() blas64.Symmetric {
return s.mat
}
// SetRawSymmetric sets the underlying blas64.Symmetric used by the receiver.
// Changes to elements in the receiver following the call will be reflected
// in the input.
//
// The supplied Symmetric must use blas.Upper storage format.
func (s *SymDense) SetRawSymmetric(mat blas64.Symmetric) {
if mat.Uplo != blas.Upper {
panic(badSymTriangle)
}
s.cap = mat.N
s.mat = mat
}
// Reset empties the matrix so that it can be reused as the
// receiver of a dimensionally restricted operation.
//
// Reset should not be used when the matrix shares backing data.
// See the Reseter interface for more information.
func (s *SymDense) Reset() {
// N and Stride must be zeroed in unison.
s.mat.N, s.mat.Stride = 0, 0
s.mat.Data = s.mat.Data[:0]
}
// ReuseAsSym changes the receiver if it IsEmpty() to be of size n×n.
//
// ReuseAsSym re-uses the backing data slice if it has sufficient capacity,
// otherwise a new slice is allocated. The backing data is zero on return.
//
// ReuseAsSym panics if the receiver is not empty, and panics if
// the input size is less than one. To empty the receiver for re-use,
// Reset should be used.
func (s *SymDense) ReuseAsSym(n int) {
if n <= 0 {
if n == 0 {
panic(ErrZeroLength)
}
panic(ErrNegativeDimension)
}
if !s.IsEmpty() {
panic(ErrReuseNonEmpty)
}
s.reuseAsZeroed(n)
}
// Zero sets all of the matrix elements to zero.
func (s *SymDense) Zero() {
for i := 0; i < s.mat.N; i++ {
zero(s.mat.Data[i*s.mat.Stride+i : i*s.mat.Stride+s.mat.N])
}
}
// IsEmpty returns whether the receiver is empty. Empty matrices can be the
// receiver for size-restricted operations. The receiver can be emptied using
// Reset.
func (s *SymDense) IsEmpty() bool {
// It must be the case that m.Dims() returns
// zeros in this case. See comment in Reset().
return s.mat.N == 0
}
// reuseAsNonZeroed resizes an empty matrix to a n×n matrix,
// or checks that a non-empty matrix is n×n.
func (s *SymDense) reuseAsNonZeroed(n int) {
// reuseAsNonZeroed must be kept in sync with reuseAsZeroed.
if n == 0 {
panic(ErrZeroLength)
}
if s.mat.N > s.cap {
// Panic as a string, not a mat.Error.
panic(badCap)
}
if s.IsEmpty() {
s.mat = blas64.Symmetric{
N: n,
Stride: n,
Data: use(s.mat.Data, n*n),
Uplo: blas.Upper,
}
s.cap = n
return
}
if s.mat.Uplo != blas.Upper {
panic(badSymTriangle)
}
if s.mat.N != n {
panic(ErrShape)
}
}
// reuseAsNonZeroed resizes an empty matrix to a n×n matrix,
// or checks that a non-empty matrix is n×n. It then zeros the
// elements of the matrix.
func (s *SymDense) reuseAsZeroed(n int) {
// reuseAsZeroed must be kept in sync with reuseAsNonZeroed.
if n == 0 {
panic(ErrZeroLength)
}
if s.mat.N > s.cap {
// Panic as a string, not a mat.Error.
panic(badCap)
}
if s.IsEmpty() {
s.mat = blas64.Symmetric{
N: n,
Stride: n,
Data: useZeroed(s.mat.Data, n*n),
Uplo: blas.Upper,
}
s.cap = n
return
}
if s.mat.Uplo != blas.Upper {
panic(badSymTriangle)
}
if s.mat.N != n {
panic(ErrShape)
}
s.Zero()
}
func (s *SymDense) isolatedWorkspace(a Symmetric) (w *SymDense, restore func()) {
n := a.Symmetric()
if n == 0 {
panic(ErrZeroLength)
}
w = getWorkspaceSym(n, false)
return w, func() {
s.CopySym(w)
putWorkspaceSym(w)
}
}
// DiagView returns the diagonal as a matrix backed by the original data.
func (s *SymDense) DiagView() Diagonal {
n := s.mat.N
return &DiagDense{
mat: blas64.Vector{
N: n,
Inc: s.mat.Stride + 1,
Data: s.mat.Data[:(n-1)*s.mat.Stride+n],
},
}
}
func (s *SymDense) AddSym(a, b Symmetric) {
n := a.Symmetric()
if n != b.Symmetric() {
panic(ErrShape)
}
s.reuseAsNonZeroed(n)
if a, ok := a.(RawSymmetricer); ok {
if b, ok := b.(RawSymmetricer); ok {
amat, bmat := a.RawSymmetric(), b.RawSymmetric()
if s != a {
s.checkOverlap(generalFromSymmetric(amat))
}
if s != b {
s.checkOverlap(generalFromSymmetric(bmat))
}
for i := 0; i < n; i++ {
btmp := bmat.Data[i*bmat.Stride+i : i*bmat.Stride+n]
stmp := s.mat.Data[i*s.mat.Stride+i : i*s.mat.Stride+n]
for j, v := range amat.Data[i*amat.Stride+i : i*amat.Stride+n] {
stmp[j] = v + btmp[j]
}
}
return
}
}
s.checkOverlapMatrix(a)
s.checkOverlapMatrix(b)
for i := 0; i < n; i++ {
stmp := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
for j := i; j < n; j++ {
stmp[j] = a.At(i, j) + b.At(i, j)
}
}
}
func (s *SymDense) CopySym(a Symmetric) int {
n := a.Symmetric()
n = min(n, s.mat.N)
if n == 0 {
return 0
}
switch a := a.(type) {
case RawSymmetricer:
amat := a.RawSymmetric()
if amat.Uplo != blas.Upper {
panic(badSymTriangle)
}
for i := 0; i < n; i++ {
copy(s.mat.Data[i*s.mat.Stride+i:i*s.mat.Stride+n], amat.Data[i*amat.Stride+i:i*amat.Stride+n])
}
default:
for i := 0; i < n; i++ {
stmp := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
for j := i; j < n; j++ {
stmp[j] = a.At(i, j)
}
}
}
return n
}
// SymRankOne performs a symmetric rank-one update to the matrix a with x,
// which is treated as a column vector, and stores the result in the receiver
// s = a + alpha * x * xᵀ
func (s *SymDense) SymRankOne(a Symmetric, alpha float64, x Vector) {
n := x.Len()
if a.Symmetric() != n {
panic(ErrShape)
}
s.reuseAsNonZeroed(n)
if s != a {
if rs, ok := a.(RawSymmetricer); ok {
s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
}
s.CopySym(a)
}
xU, _ := untransposeExtract(x)
if rv, ok := xU.(*VecDense); ok {
r, c := xU.Dims()
xmat := rv.mat
s.checkOverlap(generalFromVector(xmat, r, c))
blas64.Syr(alpha, xmat, s.mat)
return
}
for i := 0; i < n; i++ {
for j := i; j < n; j++ {
s.set(i, j, s.at(i, j)+alpha*x.AtVec(i)*x.AtVec(j))
}
}
}
// SymRankK performs a symmetric rank-k update to the matrix a and stores the
// result into the receiver. If a is zero, see SymOuterK.
// s = a + alpha * x * x'
func (s *SymDense) SymRankK(a Symmetric, alpha float64, x Matrix) {
n := a.Symmetric()
r, _ := x.Dims()
if r != n {
panic(ErrShape)
}
xMat, aTrans := untransposeExtract(x)
var g blas64.General
if rm, ok := xMat.(*Dense); ok {
g = rm.mat
} else {
g = DenseCopyOf(x).mat
aTrans = false
}
if a != s {
if rs, ok := a.(RawSymmetricer); ok {
s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
}
s.reuseAsNonZeroed(n)
s.CopySym(a)
}
t := blas.NoTrans
if aTrans {
t = blas.Trans
}
blas64.Syrk(t, alpha, g, 1, s.mat)
}
// SymOuterK calculates the outer product of x with itself and stores
// the result into the receiver. It is equivalent to the matrix
// multiplication
// s = alpha * x * x'.
// In order to update an existing matrix, see SymRankOne.
func (s *SymDense) SymOuterK(alpha float64, x Matrix) {
n, _ := x.Dims()
switch {
case s.IsEmpty():
s.mat = blas64.Symmetric{
N: n,
Stride: n,
Data: useZeroed(s.mat.Data, n*n),
Uplo: blas.Upper,
}
s.cap = n
s.SymRankK(s, alpha, x)
case s.mat.Uplo != blas.Upper:
panic(badSymTriangle)
case s.mat.N == n:
if s == x {
w := getWorkspaceSym(n, true)
w.SymRankK(w, alpha, x)
s.CopySym(w)
putWorkspaceSym(w)
} else {
switch r := x.(type) {
case RawMatrixer:
s.checkOverlap(r.RawMatrix())
case RawSymmetricer:
s.checkOverlap(generalFromSymmetric(r.RawSymmetric()))
case RawTriangular:
s.checkOverlap(generalFromTriangular(r.RawTriangular()))
}
// Only zero the upper triangle.
for i := 0; i < n; i++ {
ri := i * s.mat.Stride
zero(s.mat.Data[ri+i : ri+n])
}
s.SymRankK(s, alpha, x)
}
default:
panic(ErrShape)
}
}
// RankTwo performs a symmetric rank-two update to the matrix a with the
// vectors x and y, which are treated as column vectors, and stores the
// result in the receiver
// m = a + alpha * (x * yᵀ + y * xᵀ)
func (s *SymDense) RankTwo(a Symmetric, alpha float64, x, y Vector) {
n := s.mat.N
if x.Len() != n {
panic(ErrShape)
}
if y.Len() != n {
panic(ErrShape)
}
if s != a {
if rs, ok := a.(RawSymmetricer); ok {
s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
}
}
var xmat, ymat blas64.Vector
fast := true
xU, _ := untransposeExtract(x)
if rv, ok := xU.(*VecDense); ok {
r, c := xU.Dims()
xmat = rv.mat
s.checkOverlap(generalFromVector(xmat, r, c))
} else {
fast = false
}
yU, _ := untransposeExtract(y)
if rv, ok := yU.(*VecDense); ok {
r, c := yU.Dims()
ymat = rv.mat
s.checkOverlap(generalFromVector(ymat, r, c))
} else {
fast = false
}
if s != a {
if rs, ok := a.(RawSymmetricer); ok {
s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
}
s.reuseAsNonZeroed(n)
s.CopySym(a)
}
if fast {
if s != a {
s.reuseAsNonZeroed(n)
s.CopySym(a)
}
blas64.Syr2(alpha, xmat, ymat, s.mat)
return
}
for i := 0; i < n; i++ {
s.reuseAsNonZeroed(n)
for j := i; j < n; j++ {
s.set(i, j, a.At(i, j)+alpha*(x.AtVec(i)*y.AtVec(j)+y.AtVec(i)*x.AtVec(j)))
}
}
}
// ScaleSym multiplies the elements of a by f, placing the result in the receiver.
func (s *SymDense) ScaleSym(f float64, a Symmetric) {
n := a.Symmetric()
s.reuseAsNonZeroed(n)
if a, ok := a.(RawSymmetricer); ok {
amat := a.RawSymmetric()
if s != a {
s.checkOverlap(generalFromSymmetric(amat))
}
for i := 0; i < n; i++ {
for j := i; j < n; j++ {
s.mat.Data[i*s.mat.Stride+j] = f * amat.Data[i*amat.Stride+j]
}
}
return
}
for i := 0; i < n; i++ {
for j := i; j < n; j++ {
s.mat.Data[i*s.mat.Stride+j] = f * a.At(i, j)
}
}
}
// SubsetSym extracts a subset of the rows and columns of the matrix a and stores
// the result in-place into the receiver. The resulting matrix size is
// len(set)×len(set). Specifically, at the conclusion of SubsetSym,
// s.At(i, j) equals a.At(set[i], set[j]). Note that the supplied set does not
// have to be a strict subset, dimension repeats are allowed.
func (s *SymDense) SubsetSym(a Symmetric, set []int) {
n := len(set)
na := a.Symmetric()
s.reuseAsNonZeroed(n)
var restore func()
if a == s {
s, restore = s.isolatedWorkspace(a)
defer restore()
}
if a, ok := a.(RawSymmetricer); ok {
raw := a.RawSymmetric()
if s != a {
s.checkOverlap(generalFromSymmetric(raw))
}
for i := 0; i < n; i++ {
ssub := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
r := set[i]
rsub := raw.Data[r*raw.Stride : r*raw.Stride+na]
for j := i; j < n; j++ {
c := set[j]
if r <= c {
ssub[j] = rsub[c]
} else {
ssub[j] = raw.Data[c*raw.Stride+r]
}
}
}
return
}
for i := 0; i < n; i++ {
for j := i; j < n; j++ {
s.mat.Data[i*s.mat.Stride+j] = a.At(set[i], set[j])
}
}
}
// SliceSym returns a new Matrix that shares backing data with the receiver.
// The returned matrix starts at {i,i} of the receiver and extends k-i rows
// and columns. The final row and column in the resulting matrix is k-1.
// SliceSym panics with ErrIndexOutOfRange if the slice is outside the
// capacity of the receiver.
func (s *SymDense) SliceSym(i, k int) Symmetric {
return s.sliceSym(i, k)
}
func (s *SymDense) sliceSym(i, k int) *SymDense {
sz := s.cap
if i < 0 || sz < i || k < i || sz < k {
panic(ErrIndexOutOfRange)
}
v := *s
v.mat.Data = s.mat.Data[i*s.mat.Stride+i : (k-1)*s.mat.Stride+k]
v.mat.N = k - i
v.cap = s.cap - i
return &v
}
// Trace returns the trace of the matrix.
func (s *SymDense) Trace() float64 {
// TODO(btracey): could use internal asm sum routine.
var v float64
for i := 0; i < s.mat.N; i++ {
v += s.mat.Data[i*s.mat.Stride+i]
}
return v
}
// GrowSym returns the receiver expanded by n rows and n columns. If the
// dimensions of the expanded matrix are outside the capacity of the receiver
// a new allocation is made, otherwise not. Note that the receiver itself is
// not modified during the call to GrowSquare.
func (s *SymDense) GrowSym(n int) Symmetric {
if n < 0 {
panic(ErrIndexOutOfRange)
}
if n == 0 {
return s
}
var v SymDense
n += s.mat.N
if s.IsEmpty() || n > s.cap {
v.mat = blas64.Symmetric{
N: n,
Stride: n,
Uplo: blas.Upper,
Data: make([]float64, n*n),
}
v.cap = n
// Copy elements, including those not currently visible. Use a temporary
// structure to avoid modifying the receiver.
var tmp SymDense
tmp.mat = blas64.Symmetric{
N: s.cap,
Stride: s.mat.Stride,
Data: s.mat.Data,
Uplo: s.mat.Uplo,
}
tmp.cap = s.cap
v.CopySym(&tmp)
return &v
}
v.mat = blas64.Symmetric{
N: n,
Stride: s.mat.Stride,
Uplo: blas.Upper,
Data: s.mat.Data[:(n-1)*s.mat.Stride+n],
}
v.cap = s.cap
return &v
}
// PowPSD computes a^pow where a is a positive symmetric definite matrix.
//
// PowPSD returns an error if the matrix is not not positive symmetric definite
// or the Eigen decomposition is not successful.
func (s *SymDense) PowPSD(a Symmetric, pow float64) error {
dim := a.Symmetric()
s.reuseAsNonZeroed(dim)
var eigen EigenSym
ok := eigen.Factorize(a, true)
if !ok {
return ErrFailedEigen
}
values := eigen.Values(nil)
for i, v := range values {
if v <= 0 {
return ErrNotPSD
}
values[i] = math.Pow(v, pow)
}
var u Dense
eigen.VectorsTo(&u)
s.SymOuterK(values[0], u.ColView(0))
var v VecDense
for i := 1; i < dim; i++ {
v.ColViewOf(&u, i)
s.SymRankOne(s, values[i], &v)
}
return nil
}