vendor/gonum.org/v1/gonum/mat/symmetric.go

   1 // Copyright ©2015 The Gonum Authors. All rights reserved.
   2 // Use of this source code is governed by a BSD-style
   3 // license that can be found in the LICENSE file.
   4
   5 package mat
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas"
  11         "gonum.org/v1/gonum/blas/blas64"
  12 )
  13
  14 var (
  15         symDense *SymDense
  16
  17         _ Matrix           = symDense
  18         _ Symmetric        = symDense
  19         _ RawSymmetricer   = symDense
  20         _ MutableSymmetric = symDense
  21 )
  22
  23 const (
  24         badSymTriangle = "mat: blas64.Symmetric not upper"
  25         badSymCap      = "mat: bad capacity for SymDense"
  26 )
  27
  28 // SymDense is a symmetric matrix that uses dense storage. SymDense
  29 // matrices are stored in the upper triangle.
  30 type SymDense struct {
  31         mat blas64.Symmetric
  32         cap int
  33 }
  34
  35 // Symmetric represents a symmetric matrix (where the element at {i, j} equals
  36 // the element at {j, i}). Symmetric matrices are always square.
  37 type Symmetric interface {
  38         Matrix
  39         // Symmetric returns the number of rows/columns in the matrix.
  40         Symmetric() int
  41 }
  42
  43 // A RawSymmetricer can return a view of itself as a BLAS Symmetric matrix.
  44 type RawSymmetricer interface {
  45         RawSymmetric() blas64.Symmetric
  46 }
  47
  48 // A MutableSymmetric can set elements of a symmetric matrix.
  49 type MutableSymmetric interface {
  50         Symmetric
  51         SetSym(i, j int, v float64)
  52 }
  53
  54 // NewSymDense creates a new Symmetric matrix with n rows and columns. If data == nil,
  55 // a new slice is allocated for the backing slice. If len(data) == n*n, data is
  56 // used as the backing slice, and changes to the elements of the returned SymDense
  57 // will be reflected in data. If neither of these is true, NewSymDense will panic.
  58 //
  59 // The data must be arranged in row-major order, i.e. the (i*c + j)-th
  60 // element in the data slice is the {i, j}-th element in the matrix.
  61 // Only the values in the upper triangular portion of the matrix are used.
  62 func NewSymDense(n int, data []float64) *SymDense {
  63         if n < 0 {
  64                 panic("mat: negative dimension")
  65         }
  66         if data != nil && n*n != len(data) {
  67                 panic(ErrShape)
  68         }
  69         if data == nil {
  70                 data = make([]float64, n*n)
  71         }
  72         return &SymDense{
  73                 mat: blas64.Symmetric{
  74                         N:      n,
  75                         Stride: n,
  76                         Data:   data,
  77                         Uplo:   blas.Upper,
  78                 },
  79                 cap: n,
  80         }
  81 }
  82
  83 // Dims returns the number of rows and columns in the matrix.
  84 func (s *SymDense) Dims() (r, c int) {
  85         return s.mat.N, s.mat.N
  86 }
  87
  88 // Caps returns the number of rows and columns in the backing matrix.
  89 func (s *SymDense) Caps() (r, c int) {
  90         return s.cap, s.cap
  91 }
  92
  93 // T implements the Matrix interface. Symmetric matrices, by definition, are
  94 // equal to their transpose, and this is a no-op.
  95 func (s *SymDense) T() Matrix {
  96         return s
  97 }
  98
  99 func (s *SymDense) Symmetric() int {
 100         return s.mat.N
 101 }
 102
 103 // RawSymmetric returns the matrix as a blas64.Symmetric. The returned
 104 // value must be stored in upper triangular format.
 105 func (s *SymDense) RawSymmetric() blas64.Symmetric {
 106         return s.mat
 107 }
 108
 109 // SetRawSymmetric sets the underlying blas64.Symmetric used by the receiver.
 110 // Changes to elements in the receiver following the call will be reflected
 111 // in b. SetRawSymmetric will panic if b is not an upper-encoded symmetric
 112 // matrix.
 113 func (s *SymDense) SetRawSymmetric(b blas64.Symmetric) {
 114         if b.Uplo != blas.Upper {
 115                 panic(badSymTriangle)
 116         }
 117         s.mat = b
 118 }
 119
 120 // Reset zeros the dimensions of the matrix so that it can be reused as the
 121 // receiver of a dimensionally restricted operation.
 122 //
 123 // See the Reseter interface for more information.
 124 func (s *SymDense) Reset() {
 125         // N and Stride must be zeroed in unison.
 126         s.mat.N, s.mat.Stride = 0, 0
 127         s.mat.Data = s.mat.Data[:0]
 128 }
 129
 130 // IsZero returns whether the receiver is zero-sized. Zero-sized matrices can be the
 131 // receiver for size-restricted operations. SymDense matrices can be zeroed using Reset.
 132 func (s *SymDense) IsZero() bool {
 133         // It must be the case that m.Dims() returns
 134         // zeros in this case. See comment in Reset().
 135         return s.mat.N == 0
 136 }
 137
 138 // reuseAs resizes an empty matrix to a n×n matrix,
 139 // or checks that a non-empty matrix is n×n.
 140 func (s *SymDense) reuseAs(n int) {
 141         if s.mat.N > s.cap {
 142                 panic(badSymCap)
 143         }
 144         if s.IsZero() {
 145                 s.mat = blas64.Symmetric{
 146                         N:      n,
 147                         Stride: n,
 148                         Data:   use(s.mat.Data, n*n),
 149                         Uplo:   blas.Upper,
 150                 }
 151                 s.cap = n
 152                 return
 153         }
 154         if s.mat.Uplo != blas.Upper {
 155                 panic(badSymTriangle)
 156         }
 157         if s.mat.N != n {
 158                 panic(ErrShape)
 159         }
 160 }
 161
 162 func (s *SymDense) isolatedWorkspace(a Symmetric) (w *SymDense, restore func()) {
 163         n := a.Symmetric()
 164         w = getWorkspaceSym(n, false)
 165         return w, func() {
 166                 s.CopySym(w)
 167                 putWorkspaceSym(w)
 168         }
 169 }
 170
 171 func (s *SymDense) AddSym(a, b Symmetric) {
 172         n := a.Symmetric()
 173         if n != b.Symmetric() {
 174                 panic(ErrShape)
 175         }
 176         s.reuseAs(n)
 177
 178         if a, ok := a.(RawSymmetricer); ok {
 179                 if b, ok := b.(RawSymmetricer); ok {
 180                         amat, bmat := a.RawSymmetric(), b.RawSymmetric()
 181                         if s != a {
 182                                 s.checkOverlap(generalFromSymmetric(amat))
 183                         }
 184                         if s != b {
 185                                 s.checkOverlap(generalFromSymmetric(bmat))
 186                         }
 187                         for i := 0; i < n; i++ {
 188                                 btmp := bmat.Data[i*bmat.Stride+i : i*bmat.Stride+n]
 189                                 stmp := s.mat.Data[i*s.mat.Stride+i : i*s.mat.Stride+n]
 190                                 for j, v := range amat.Data[i*amat.Stride+i : i*amat.Stride+n] {
 191                                         stmp[j] = v + btmp[j]
 192                                 }
 193                         }
 194                         return
 195                 }
 196         }
 197
 198         for i := 0; i < n; i++ {
 199                 stmp := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
 200                 for j := i; j < n; j++ {
 201                         stmp[j] = a.At(i, j) + b.At(i, j)
 202                 }
 203         }
 204 }
 205
 206 func (s *SymDense) CopySym(a Symmetric) int {
 207         n := a.Symmetric()
 208         n = min(n, s.mat.N)
 209         if n == 0 {
 210                 return 0
 211         }
 212         switch a := a.(type) {
 213         case RawSymmetricer:
 214                 amat := a.RawSymmetric()
 215                 if amat.Uplo != blas.Upper {
 216                         panic(badSymTriangle)
 217                 }
 218                 for i := 0; i < n; i++ {
 219                         copy(s.mat.Data[i*s.mat.Stride+i:i*s.mat.Stride+n], amat.Data[i*amat.Stride+i:i*amat.Stride+n])
 220                 }
 221         default:
 222                 for i := 0; i < n; i++ {
 223                         stmp := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
 224                         for j := i; j < n; j++ {
 225                                 stmp[j] = a.At(i, j)
 226                         }
 227                 }
 228         }
 229         return n
 230 }
 231
 232 // SymRankOne performs a symetric rank-one update to the matrix a and stores
 233 // the result in the receiver
 234 //  s = a + alpha * x * x'
 235 func (s *SymDense) SymRankOne(a Symmetric, alpha float64, x Vector) {
 236         n, c := x.Dims()
 237         if a.Symmetric() != n || c != 1 {
 238                 panic(ErrShape)
 239         }
 240         s.reuseAs(n)
 241
 242         if s != a {
 243                 if rs, ok := a.(RawSymmetricer); ok {
 244                         s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
 245                 }
 246                 s.CopySym(a)
 247         }
 248
 249         xU, _ := untranspose(x)
 250         if rv, ok := xU.(RawVectorer); ok {
 251                 xmat := rv.RawVector()
 252                 s.checkOverlap((&VecDense{mat: xmat, n: n}).asGeneral())
 253                 blas64.Syr(alpha, xmat, s.mat)
 254                 return
 255         }
 256
 257         for i := 0; i < n; i++ {
 258                 for j := i; j < n; j++ {
 259                         s.set(i, j, s.at(i, j)+alpha*x.AtVec(i)*x.AtVec(j))
 260                 }
 261         }
 262 }
 263
 264 // SymRankK performs a symmetric rank-k update to the matrix a and stores the
 265 // result into the receiver. If a is zero, see SymOuterK.
 266 //  s = a + alpha * x * x'
 267 func (s *SymDense) SymRankK(a Symmetric, alpha float64, x Matrix) {
 268         n := a.Symmetric()
 269         r, _ := x.Dims()
 270         if r != n {
 271                 panic(ErrShape)
 272         }
 273         xMat, aTrans := untranspose(x)
 274         var g blas64.General
 275         if rm, ok := xMat.(RawMatrixer); ok {
 276                 g = rm.RawMatrix()
 277         } else {
 278                 g = DenseCopyOf(x).mat
 279                 aTrans = false
 280         }
 281         if a != s {
 282                 if rs, ok := a.(RawSymmetricer); ok {
 283                         s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
 284                 }
 285                 s.reuseAs(n)
 286                 s.CopySym(a)
 287         }
 288         t := blas.NoTrans
 289         if aTrans {
 290                 t = blas.Trans
 291         }
 292         blas64.Syrk(t, alpha, g, 1, s.mat)
 293 }
 294
 295 // SymOuterK calculates the outer product of x with itself and stores
 296 // the result into the receiver. It is equivalent to the matrix
 297 // multiplication
 298 //  s = alpha * x * x'.
 299 // In order to update an existing matrix, see SymRankOne.
 300 func (s *SymDense) SymOuterK(alpha float64, x Matrix) {
 301         n, _ := x.Dims()
 302         switch {
 303         case s.IsZero():
 304                 s.mat = blas64.Symmetric{
 305                         N:      n,
 306                         Stride: n,
 307                         Data:   useZeroed(s.mat.Data, n*n),
 308                         Uplo:   blas.Upper,
 309                 }
 310                 s.cap = n
 311                 s.SymRankK(s, alpha, x)
 312         case s.mat.Uplo != blas.Upper:
 313                 panic(badSymTriangle)
 314         case s.mat.N == n:
 315                 if s == x {
 316                         w := getWorkspaceSym(n, true)
 317                         w.SymRankK(w, alpha, x)
 318                         s.CopySym(w)
 319                         putWorkspaceSym(w)
 320                 } else {
 321                         switch r := x.(type) {
 322                         case RawMatrixer:
 323                                 s.checkOverlap(r.RawMatrix())
 324                         case RawSymmetricer:
 325                                 s.checkOverlap(generalFromSymmetric(r.RawSymmetric()))
 326                         case RawTriangular:
 327                                 s.checkOverlap(generalFromTriangular(r.RawTriangular()))
 328                         }
 329                         // Only zero the upper triangle.
 330                         for i := 0; i < n; i++ {
 331                                 ri := i * s.mat.Stride
 332                                 zero(s.mat.Data[ri+i : ri+n])
 333                         }
 334                         s.SymRankK(s, alpha, x)
 335                 }
 336         default:
 337                 panic(ErrShape)
 338         }
 339 }
 340
 341 // RankTwo performs a symmmetric rank-two update to the matrix a and stores
 342 // the result in the receiver
 343 //  m = a + alpha * (x * y' + y * x')
 344 func (s *SymDense) RankTwo(a Symmetric, alpha float64, x, y Vector) {
 345         n := s.mat.N
 346         xr, xc := x.Dims()
 347         if xr != n || xc != 1 {
 348                 panic(ErrShape)
 349         }
 350         yr, yc := y.Dims()
 351         if yr != n || yc != 1 {
 352                 panic(ErrShape)
 353         }
 354
 355         if s != a {
 356                 if rs, ok := a.(RawSymmetricer); ok {
 357                         s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
 358                 }
 359         }
 360
 361         var xmat, ymat blas64.Vector
 362         fast := true
 363         xU, _ := untranspose(x)
 364         if rv, ok := xU.(RawVectorer); ok {
 365                 xmat = rv.RawVector()
 366                 s.checkOverlap((&VecDense{mat: xmat, n: x.Len()}).asGeneral())
 367         } else {
 368                 fast = false
 369         }
 370         yU, _ := untranspose(y)
 371         if rv, ok := yU.(RawVectorer); ok {
 372                 ymat = rv.RawVector()
 373                 s.checkOverlap((&VecDense{mat: ymat, n: y.Len()}).asGeneral())
 374         } else {
 375                 fast = false
 376         }
 377
 378         if s != a {
 379                 if rs, ok := a.(RawSymmetricer); ok {
 380                         s.checkOverlap(generalFromSymmetric(rs.RawSymmetric()))
 381                 }
 382                 s.reuseAs(n)
 383                 s.CopySym(a)
 384         }
 385
 386         if fast {
 387                 if s != a {
 388                         s.reuseAs(n)
 389                         s.CopySym(a)
 390                 }
 391                 blas64.Syr2(alpha, xmat, ymat, s.mat)
 392                 return
 393         }
 394
 395         for i := 0; i < n; i++ {
 396                 s.reuseAs(n)
 397                 for j := i; j < n; j++ {
 398                         s.set(i, j, a.At(i, j)+alpha*(x.AtVec(i)*y.AtVec(j)+y.AtVec(i)*x.AtVec(j)))
 399                 }
 400         }
 401 }
 402
 403 // ScaleSym multiplies the elements of a by f, placing the result in the receiver.
 404 func (s *SymDense) ScaleSym(f float64, a Symmetric) {
 405         n := a.Symmetric()
 406         s.reuseAs(n)
 407         if a, ok := a.(RawSymmetricer); ok {
 408                 amat := a.RawSymmetric()
 409                 if s != a {
 410                         s.checkOverlap(generalFromSymmetric(amat))
 411                 }
 412                 for i := 0; i < n; i++ {
 413                         for j := i; j < n; j++ {
 414                                 s.mat.Data[i*s.mat.Stride+j] = f * amat.Data[i*amat.Stride+j]
 415                         }
 416                 }
 417                 return
 418         }
 419         for i := 0; i < n; i++ {
 420                 for j := i; j < n; j++ {
 421                         s.mat.Data[i*s.mat.Stride+j] = f * a.At(i, j)
 422                 }
 423         }
 424 }
 425
 426 // SubsetSym extracts a subset of the rows and columns of the matrix a and stores
 427 // the result in-place into the receiver. The resulting matrix size is
 428 // len(set)×len(set). Specifically, at the conclusion of SubsetSym,
 429 // s.At(i, j) equals a.At(set[i], set[j]). Note that the supplied set does not
 430 // have to be a strict subset, dimension repeats are allowed.
 431 func (s *SymDense) SubsetSym(a Symmetric, set []int) {
 432         n := len(set)
 433         na := a.Symmetric()
 434         s.reuseAs(n)
 435         var restore func()
 436         if a == s {
 437                 s, restore = s.isolatedWorkspace(a)
 438                 defer restore()
 439         }
 440
 441         if a, ok := a.(RawSymmetricer); ok {
 442                 raw := a.RawSymmetric()
 443                 if s != a {
 444                         s.checkOverlap(generalFromSymmetric(raw))
 445                 }
 446                 for i := 0; i < n; i++ {
 447                         ssub := s.mat.Data[i*s.mat.Stride : i*s.mat.Stride+n]
 448                         r := set[i]
 449                         rsub := raw.Data[r*raw.Stride : r*raw.Stride+na]
 450                         for j := i; j < n; j++ {
 451                                 c := set[j]
 452                                 if r <= c {
 453                                         ssub[j] = rsub[c]
 454                                 } else {
 455                                         ssub[j] = raw.Data[c*raw.Stride+r]
 456                                 }
 457                         }
 458                 }
 459                 return
 460         }
 461         for i := 0; i < n; i++ {
 462                 for j := i; j < n; j++ {
 463                         s.mat.Data[i*s.mat.Stride+j] = a.At(set[i], set[j])
 464                 }
 465         }
 466 }
 467
 468 // SliceSquare returns a new Matrix that shares backing data with the receiver.
 469 // The returned matrix starts at {i,i} of the receiver and extends k-i rows
 470 // and columns. The final row and column in the resulting matrix is k-1.
 471 // SliceSquare panics with ErrIndexOutOfRange if the slice is outside the capacity
 472 // of the receiver.
 473 func (s *SymDense) SliceSquare(i, k int) Matrix {
 474         sz := s.cap
 475         if i < 0 || sz < i || k < i || sz < k {
 476                 panic(ErrIndexOutOfRange)
 477         }
 478         v := *s
 479         v.mat.Data = s.mat.Data[i*s.mat.Stride+i : (k-1)*s.mat.Stride+k]
 480         v.mat.N = k - i
 481         v.cap = s.cap - i
 482         return &v
 483 }
 484
 485 // GrowSquare returns the receiver expanded by n rows and n columns. If the
 486 // dimensions of the expanded matrix are outside the capacity of the receiver
 487 // a new allocation is made, otherwise not. Note that the receiver itself is
 488 // not modified during the call to GrowSquare.
 489 func (s *SymDense) GrowSquare(n int) Matrix {
 490         if n < 0 {
 491                 panic(ErrIndexOutOfRange)
 492         }
 493         if n == 0 {
 494                 return s
 495         }
 496         var v SymDense
 497         n += s.mat.N
 498         if n > s.cap {
 499                 v.mat = blas64.Symmetric{
 500                         N:      n,
 501                         Stride: n,
 502                         Uplo:   blas.Upper,
 503                         Data:   make([]float64, n*n),
 504                 }
 505                 v.cap = n
 506                 // Copy elements, including those not currently visible. Use a temporary
 507                 // structure to avoid modifying the receiver.
 508                 var tmp SymDense
 509                 tmp.mat = blas64.Symmetric{
 510                         N:      s.cap,
 511                         Stride: s.mat.Stride,
 512                         Data:   s.mat.Data,
 513                         Uplo:   s.mat.Uplo,
 514                 }
 515                 tmp.cap = s.cap
 516                 v.CopySym(&tmp)
 517                 return &v
 518         }
 519         v.mat = blas64.Symmetric{
 520                 N:      n,
 521                 Stride: s.mat.Stride,
 522                 Uplo:   blas.Upper,
 523                 Data:   s.mat.Data[:(n-1)*s.mat.Stride+n],
 524         }
 525         v.cap = s.cap
 526         return &v
 527 }
 528
 529 // PowPSD computes a^pow where a is a positive symmetric definite matrix.
 530 //
 531 // PowPSD returns an error if the matrix is not  not positive symmetric definite
 532 // or the Eigendecomposition is not successful.
 533 func (s *SymDense) PowPSD(a Symmetric, pow float64) error {
 534         dim := a.Symmetric()
 535         s.reuseAs(dim)
 536
 537         var eigen EigenSym
 538         ok := eigen.Factorize(a, true)
 539         if !ok {
 540                 return ErrFailedEigen
 541         }
 542         values := eigen.Values(nil)
 543         for i, v := range values {
 544                 if v <= 0 {
 545                         return ErrNotPSD
 546                 }
 547                 values[i] = math.Pow(v, pow)
 548         }
 549         var u Dense
 550         u.EigenvectorsSym(&eigen)
 551
 552         s.SymOuterK(values[0], u.ColView(0))
 553
 554         var v VecDense
 555         for i := 1; i < dim; i++ {
 556                 v.ColViewOf(&u, i)
 557                 s.SymRankOne(s, values[i], &v)
 558         }
 559         return nil
 560 }