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

   1 // Copyright ©2013 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         "gonum.org/v1/gonum/lapack/lapack64"
  13 )
  14
  15 const (
  16         badTriangle = "mat: invalid triangle"
  17         badCholesky = "mat: invalid Cholesky factorization"
  18 )
  19
  20 // Cholesky is a type for creating and using the Cholesky factorization of a
  21 // symmetric positive definite matrix.
  22 //
  23 // Cholesky methods may only be called on a value that has been successfully
  24 // initialized by a call to Factorize that has returned true. Calls to methods
  25 // of an unsuccessful Cholesky factorization will panic.
  26 type Cholesky struct {
  27         // The chol pointer must never be retained as a pointer outside the Cholesky
  28         // struct, either by returning chol outside the struct or by setting it to
  29         // a pointer coming from outside. The same prohibition applies to the data
  30         // slice within chol.
  31         chol *TriDense
  32         cond float64
  33 }
  34
  35 // updateCond updates the condition number of the Cholesky decomposition. If
  36 // norm > 0, then that norm is used as the norm of the original matrix A, otherwise
  37 // the norm is estimated from the decomposition.
  38 func (c *Cholesky) updateCond(norm float64) {
  39         n := c.chol.mat.N
  40         work := getFloats(3*n, false)
  41         defer putFloats(work)
  42         if norm < 0 {
  43                 // This is an approximation. By the definition of a norm,
  44                 //  |AB| <= |A| |B|.
  45                 // Since A = U^T*U, we get for the condition number κ that
  46                 //  κ(A) := |A| |A^-1| = |U^T*U| |A^-1| <= |U^T| |U| |A^-1|,
  47                 // so this will overestimate the condition number somewhat.
  48                 // The norm of the original factorized matrix cannot be stored
  49                 // because of update possibilities.
  50                 unorm := lapack64.Lantr(CondNorm, c.chol.mat, work)
  51                 lnorm := lapack64.Lantr(CondNormTrans, c.chol.mat, work)
  52                 norm = unorm * lnorm
  53         }
  54         sym := c.chol.asSymBlas()
  55         iwork := getInts(n, false)
  56         v := lapack64.Pocon(sym, norm, work, iwork)
  57         putInts(iwork)
  58         c.cond = 1 / v
  59 }
  60
  61 // Cond returns the condition number of the factorized matrix.
  62 func (c *Cholesky) Cond() float64 {
  63         return c.cond
  64 }
  65
  66 // Factorize calculates the Cholesky decomposition of the matrix A and returns
  67 // whether the matrix is positive definite. If Factorize returns false, the
  68 // factorization must not be used.
  69 func (c *Cholesky) Factorize(a Symmetric) (ok bool) {
  70         n := a.Symmetric()
  71         if c.chol == nil {
  72                 c.chol = NewTriDense(n, Upper, nil)
  73         } else {
  74                 c.chol = NewTriDense(n, Upper, use(c.chol.mat.Data, n*n))
  75         }
  76         copySymIntoTriangle(c.chol, a)
  77
  78         sym := c.chol.asSymBlas()
  79         work := getFloats(c.chol.mat.N, false)
  80         norm := lapack64.Lansy(CondNorm, sym, work)
  81         putFloats(work)
  82         _, ok = lapack64.Potrf(sym)
  83         if ok {
  84                 c.updateCond(norm)
  85         } else {
  86                 c.Reset()
  87         }
  88         return ok
  89 }
  90
  91 // Reset resets the factorization so that it can be reused as the receiver of a
  92 // dimensionally restricted operation.
  93 func (c *Cholesky) Reset() {
  94         if c.chol != nil {
  95                 c.chol.Reset()
  96         }
  97         c.cond = math.Inf(1)
  98 }
  99
 100 // SetFromU sets the Cholesky decomposition from the given triangular matrix.
 101 // SetFromU panics if t is not upper triangular. Note that t is copied into,
 102 // not stored inside, the receiver.
 103 func (c *Cholesky) SetFromU(t *TriDense) {
 104         n, kind := t.Triangle()
 105         if kind != Upper {
 106                 panic("cholesky: matrix must be upper triangular")
 107         }
 108         if c.chol == nil {
 109                 c.chol = NewTriDense(n, Upper, nil)
 110         } else {
 111                 c.chol = NewTriDense(n, Upper, use(c.chol.mat.Data, n*n))
 112         }
 113         c.chol.Copy(t)
 114         c.updateCond(-1)
 115 }
 116
 117 // Clone makes a copy of the input Cholesky into the receiver, overwriting the
 118 // previous value of the receiver. Clone does not place any restrictions on receiver
 119 // shape. Clone panics if the input Cholesky is not the result of a valid decomposition.
 120 func (c *Cholesky) Clone(chol *Cholesky) {
 121         if !chol.valid() {
 122                 panic(badCholesky)
 123         }
 124         n := chol.Size()
 125         if c.chol == nil {
 126                 c.chol = NewTriDense(n, Upper, nil)
 127         } else {
 128                 c.chol = NewTriDense(n, Upper, use(c.chol.mat.Data, n*n))
 129         }
 130         c.chol.Copy(chol.chol)
 131         c.cond = chol.cond
 132 }
 133
 134 // Size returns the dimension of the factorized matrix.
 135 func (c *Cholesky) Size() int {
 136         if !c.valid() {
 137                 panic(badCholesky)
 138         }
 139         return c.chol.mat.N
 140 }
 141
 142 // Det returns the determinant of the matrix that has been factorized.
 143 func (c *Cholesky) Det() float64 {
 144         if !c.valid() {
 145                 panic(badCholesky)
 146         }
 147         return math.Exp(c.LogDet())
 148 }
 149
 150 // LogDet returns the log of the determinant of the matrix that has been factorized.
 151 func (c *Cholesky) LogDet() float64 {
 152         if !c.valid() {
 153                 panic(badCholesky)
 154         }
 155         var det float64
 156         for i := 0; i < c.chol.mat.N; i++ {
 157                 det += 2 * math.Log(c.chol.mat.Data[i*c.chol.mat.Stride+i])
 158         }
 159         return det
 160 }
 161
 162 // Solve finds the matrix m that solves A * m = b where A is represented
 163 // by the Cholesky decomposition, placing the result in m.
 164 func (c *Cholesky) Solve(m *Dense, b Matrix) error {
 165         if !c.valid() {
 166                 panic(badCholesky)
 167         }
 168         n := c.chol.mat.N
 169         bm, bn := b.Dims()
 170         if n != bm {
 171                 panic(ErrShape)
 172         }
 173
 174         m.reuseAs(bm, bn)
 175         if b != m {
 176                 m.Copy(b)
 177         }
 178         blas64.Trsm(blas.Left, blas.Trans, 1, c.chol.mat, m.mat)
 179         blas64.Trsm(blas.Left, blas.NoTrans, 1, c.chol.mat, m.mat)
 180         if c.cond > ConditionTolerance {
 181                 return Condition(c.cond)
 182         }
 183         return nil
 184 }
 185
 186 // SolveChol finds the matrix m that solves A * m = B where A and B are represented
 187 // by their Cholesky decompositions a and b, placing the result in the receiver.
 188 func (a *Cholesky) SolveChol(m *Dense, b *Cholesky) error {
 189         if !a.valid() || !b.valid() {
 190                 panic(badCholesky)
 191         }
 192         bn := b.chol.mat.N
 193         if a.chol.mat.N != bn {
 194                 panic(ErrShape)
 195         }
 196
 197         m.reuseAsZeroed(bn, bn)
 198         m.Copy(b.chol.T())
 199         blas64.Trsm(blas.Left, blas.Trans, 1, a.chol.mat, m.mat)
 200         blas64.Trsm(blas.Left, blas.NoTrans, 1, a.chol.mat, m.mat)
 201         blas64.Trmm(blas.Right, blas.NoTrans, 1, b.chol.mat, m.mat)
 202         if a.cond > ConditionTolerance {
 203                 return Condition(a.cond)
 204         }
 205         return nil
 206 }
 207
 208 // SolveVec finds the vector v that solves A * v = b where A is represented
 209 // by the Cholesky decomposition, placing the result in v.
 210 func (c *Cholesky) SolveVec(v *VecDense, b Vector) error {
 211         if !c.valid() {
 212                 panic(badCholesky)
 213         }
 214         n := c.chol.mat.N
 215         if br, bc := b.Dims(); br != n || bc != 1 {
 216                 panic(ErrShape)
 217         }
 218         switch rv := b.(type) {
 219         default:
 220                 v.reuseAs(n)
 221                 return c.Solve(v.asDense(), b)
 222         case RawVectorer:
 223                 bmat := rv.RawVector()
 224                 if v != b {
 225                         v.checkOverlap(bmat)
 226                 }
 227                 v.reuseAs(n)
 228                 if v != b {
 229                         v.CopyVec(b)
 230                 }
 231                 blas64.Trsv(blas.Trans, c.chol.mat, v.mat)
 232                 blas64.Trsv(blas.NoTrans, c.chol.mat, v.mat)
 233                 if c.cond > ConditionTolerance {
 234                         return Condition(c.cond)
 235                 }
 236                 return nil
 237         }
 238 }
 239
 240 // RawU returns the Triangular matrix used to store the Cholesky decomposition of
 241 // the original matrix A. The returned matrix should not be modified. If it is
 242 // modified, the decomposition is invalid and should not be used.
 243 func (c *Cholesky) RawU() Triangular {
 244         return c.chol
 245 }
 246
 247 // UTo extracts the n×n upper triangular matrix U from a Cholesky
 248 // decomposition into dst and returns the result. If dst is nil a new
 249 // TriDense is allocated.
 250 //  A = U^T * U.
 251 func (c *Cholesky) UTo(dst *TriDense) *TriDense {
 252         if !c.valid() {
 253                 panic(badCholesky)
 254         }
 255         n := c.chol.mat.N
 256         if dst == nil {
 257                 dst = NewTriDense(n, Upper, make([]float64, n*n))
 258         } else {
 259                 dst.reuseAs(n, Upper)
 260         }
 261         dst.Copy(c.chol)
 262         return dst
 263 }
 264
 265 // LTo extracts the n×n lower triangular matrix L from a Cholesky
 266 // decomposition into dst and returns the result. If dst is nil a new
 267 // TriDense is allocated.
 268 //  A = L * L^T.
 269 func (c *Cholesky) LTo(dst *TriDense) *TriDense {
 270         if !c.valid() {
 271                 panic(badCholesky)
 272         }
 273         n := c.chol.mat.N
 274         if dst == nil {
 275                 dst = NewTriDense(n, Lower, make([]float64, n*n))
 276         } else {
 277                 dst.reuseAs(n, Lower)
 278         }
 279         dst.Copy(c.chol.TTri())
 280         return dst
 281 }
 282
 283 // ToSym reconstructs the original positive definite matrix given its
 284 // Cholesky decomposition into dst and returns the result. If dst is nil
 285 // a new SymDense is allocated.
 286 func (c *Cholesky) ToSym(dst *SymDense) *SymDense {
 287         if !c.valid() {
 288                 panic(badCholesky)
 289         }
 290         n := c.chol.mat.N
 291         if dst == nil {
 292                 dst = NewSymDense(n, make([]float64, n*n))
 293         } else {
 294                 dst.reuseAs(n)
 295         }
 296         dst.SymOuterK(1, c.chol.T())
 297         return dst
 298 }
 299
 300 // InverseTo computes the inverse of the matrix represented by its Cholesky
 301 // factorization and stores the result into s. If the factorized
 302 // matrix is ill-conditioned, a Condition error will be returned.
 303 // Note that matrix inversion is numerically unstable, and should generally be
 304 // avoided where possible, for example by using the Solve routines.
 305 func (c *Cholesky) InverseTo(s *SymDense) error {
 306         if !c.valid() {
 307                 panic(badCholesky)
 308         }
 309         // TODO(btracey): Replace this code with a direct call to Dpotri when it
 310         // is available.
 311         s.reuseAs(c.chol.mat.N)
 312         // If:
 313         //  chol(A) = U^T * U
 314         // Then:
 315         //  chol(A^-1) = S * S^T
 316         // where S = U^-1
 317         var t TriDense
 318         err := t.InverseTri(c.chol)
 319         s.SymOuterK(1, &t)
 320         return err
 321 }
 322
 323 // Scale multiplies the original matrix A by a positive constant using
 324 // its Cholesky decomposition, storing the result in-place into the receiver.
 325 // That is, if the original Cholesky factorization is
 326 //  U^T * U = A
 327 // the updated factorization is
 328 //  U'^T * U' = f A = A'
 329 // Scale panics if the constant is non-positive, or if the receiver is non-zero
 330 // and is of a different Size from the input.
 331 func (c *Cholesky) Scale(f float64, orig *Cholesky) {
 332         if !orig.valid() {
 333                 panic(badCholesky)
 334         }
 335         if f <= 0 {
 336                 panic("cholesky: scaling by a non-positive constant")
 337         }
 338         n := orig.Size()
 339         if c.chol == nil {
 340                 c.chol = NewTriDense(n, Upper, nil)
 341         } else if c.chol.mat.N != n {
 342                 panic(ErrShape)
 343         }
 344         c.chol.ScaleTri(math.Sqrt(f), orig.chol)
 345         c.cond = orig.cond // Scaling by a positive constant does not change the condition number.
 346 }
 347
 348 // ExtendVecSym computes the Cholesky decomposition of the original matrix A,
 349 // whose Cholesky decomposition is in a, extended by a the n×1 vector v according to
 350 //  [A  w]
 351 //  [w' k]
 352 // where k = v[n-1] and w = v[:n-1]. The result is stored into the receiver.
 353 // In order for the updated matrix to be positive definite, it must be the case
 354 // that k > w' A^-1 w. If this condition does not hold then ExtendVecSym will
 355 // return false and the receiver will not be updated.
 356 //
 357 // ExtendVecSym will panic if v.Len() != a.Size()+1 or if a does not contain
 358 // a valid decomposition.
 359 func (chol *Cholesky) ExtendVecSym(a *Cholesky, v Vector) (ok bool) {
 360         n := a.Size()
 361         if v.Len() != n+1 {
 362                 panic(badSliceLength)
 363         }
 364         if !a.valid() {
 365                 panic(badCholesky)
 366         }
 367
 368         // The algorithm is commented here, but see also
 369         //  https://math.stackexchange.com/questions/955874/cholesky-factor-when-adding-a-row-and-column-to-already-factorized-matrix
 370         // We have A and want to compute the Cholesky of
 371         //  [A  w]
 372         //  [w' k]
 373         // We want
 374         //  [U c]
 375         //  [0 d]
 376         // to be the updated Cholesky, and so it must be that
 377         //  [A  w] = [U' 0] [U c]
 378         //  [w' k]   [c' d] [0 d]
 379         // Thus, we need
 380         //  1) A = U'U (true by the original decomposition being valid),
 381         //  2) U' * c = w  =>  c = U'^-1 w
 382         //  3) c'*c + d'*d = k  =>  d = sqrt(k-c'*c)
 383
 384         // First, compute c = U'^-1 a
 385         // TODO(btracey): Replace this with CopyVec when issue 167 is fixed.
 386         w := NewVecDense(n, nil)
 387         for i := 0; i < n; i++ {
 388                 w.SetVec(i, v.At(i, 0))
 389         }
 390         k := v.At(n, 0)
 391
 392         c := NewVecDense(n, nil)
 393         c.SolveVec(a.chol.T(), w)
 394
 395         dot := Dot(c, c)
 396         if dot >= k {
 397                 return false
 398         }
 399         d := math.Sqrt(k - dot)
 400
 401         newU := NewTriDense(n+1, Upper, nil)
 402         newU.Copy(a.chol)
 403         for i := 0; i < n; i++ {
 404                 newU.SetTri(i, n, c.At(i, 0))
 405         }
 406         newU.SetTri(n, n, d)
 407         chol.chol = newU
 408         chol.updateCond(-1)
 409         return true
 410 }
 411
 412 // SymRankOne performs a rank-1 update of the original matrix A and refactorizes
 413 // its Cholesky factorization, storing the result into the receiver. That is, if
 414 // in the original Cholesky factorization
 415 //  U^T * U = A,
 416 // in the updated factorization
 417 //  U'^T * U' = A + alpha * x * x^T = A'.
 418 //
 419 // Note that when alpha is negative, the updating problem may be ill-conditioned
 420 // and the results may be inaccurate, or the updated matrix A' may not be
 421 // positive definite and not have a Cholesky factorization. SymRankOne returns
 422 // whether the updated matrix A' is positive definite.
 423 //
 424 // SymRankOne updates a Cholesky factorization in O(n²) time. The Cholesky
 425 // factorization computation from scratch is O(n³).
 426 func (c *Cholesky) SymRankOne(orig *Cholesky, alpha float64, x Vector) (ok bool) {
 427         if !orig.valid() {
 428                 panic(badCholesky)
 429         }
 430         n := orig.Size()
 431         if r, c := x.Dims(); r != n || c != 1 {
 432                 panic(ErrShape)
 433         }
 434         if orig != c {
 435                 if c.chol == nil {
 436                         c.chol = NewTriDense(n, Upper, nil)
 437                 } else if c.chol.mat.N != n {
 438                         panic(ErrShape)
 439                 }
 440                 c.chol.Copy(orig.chol)
 441         }
 442
 443         if alpha == 0 {
 444                 return true
 445         }
 446
 447         // Algorithms for updating and downdating the Cholesky factorization are
 448         // described, for example, in
 449         // - J. J. Dongarra, J. R. Bunch, C. B. Moler, G. W. Stewart: LINPACK
 450         //   Users' Guide. SIAM (1979), pages 10.10--10.14
 451         // or
 452         // - P. E. Gill, G. H. Golub, W. Murray, and M. A. Saunders: Methods for
 453         //   modifying matrix factorizations. Mathematics of Computation 28(126)
 454         //   (1974), Method C3 on page 521
 455         //
 456         // The implementation is based on LINPACK code
 457         // http://www.netlib.org/linpack/dchud.f
 458         // http://www.netlib.org/linpack/dchdd.f
 459         // and
 460         // https://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=2646
 461         //
 462         // According to http://icl.cs.utk.edu/lapack-forum/archives/lapack/msg00301.html
 463         // LINPACK is released under BSD license.
 464         //
 465         // See also:
 466         // - M. A. Saunders: Large-scale Linear Programming Using the Cholesky
 467         //   Factorization. Technical Report Stanford University (1972)
 468         //   http://i.stanford.edu/pub/cstr/reports/cs/tr/72/252/CS-TR-72-252.pdf
 469         // - Matthias Seeger: Low rank updates for the Cholesky decomposition.
 470         //   EPFL Technical Report 161468 (2004)
 471         //   http://infoscience.epfl.ch/record/161468
 472
 473         work := getFloats(n, false)
 474         defer putFloats(work)
 475         var xmat blas64.Vector
 476         if rv, ok := x.(RawVectorer); ok {
 477                 xmat = rv.RawVector()
 478         } else {
 479                 var tmp *VecDense
 480                 tmp.CopyVec(x)
 481                 xmat = tmp.RawVector()
 482         }
 483         blas64.Copy(n, xmat, blas64.Vector{1, work})
 484
 485         if alpha > 0 {
 486                 // Compute rank-1 update.
 487                 if alpha != 1 {
 488                         blas64.Scal(n, math.Sqrt(alpha), blas64.Vector{1, work})
 489                 }
 490                 umat := c.chol.mat
 491                 stride := umat.Stride
 492                 for i := 0; i < n; i++ {
 493                         // Compute parameters of the Givens matrix that zeroes
 494                         // the i-th element of x.
 495                         c, s, r, _ := blas64.Rotg(umat.Data[i*stride+i], work[i])
 496                         if r < 0 {
 497                                 // Multiply by -1 to have positive diagonal
 498                                 // elemnts.
 499                                 r *= -1
 500                                 c *= -1
 501                                 s *= -1
 502                         }
 503                         umat.Data[i*stride+i] = r
 504                         if i < n-1 {
 505                                 // Multiply the extended factorization matrix by
 506                                 // the Givens matrix from the left. Only
 507                                 // the i-th row and x are modified.
 508                                 blas64.Rot(n-i-1,
 509                                         blas64.Vector{1, umat.Data[i*stride+i+1 : i*stride+n]},
 510                                         blas64.Vector{1, work[i+1 : n]},
 511                                         c, s)
 512                         }
 513                 }
 514                 c.updateCond(-1)
 515                 return true
 516         }
 517
 518         // Compute rank-1 downdate.
 519         alpha = math.Sqrt(-alpha)
 520         if alpha != 1 {
 521                 blas64.Scal(n, alpha, blas64.Vector{1, work})
 522         }
 523         // Solve U^T * p = x storing the result into work.
 524         ok = lapack64.Trtrs(blas.Trans, c.chol.RawTriangular(), blas64.General{
 525                 Rows:   n,
 526                 Cols:   1,
 527                 Stride: 1,
 528                 Data:   work,
 529         })
 530         if !ok {
 531                 // The original matrix is singular. Should not happen, because
 532                 // the factorization is valid.
 533                 panic(badCholesky)
 534         }
 535         norm := blas64.Nrm2(n, blas64.Vector{1, work})
 536         if norm >= 1 {
 537                 // The updated matrix is not positive definite.
 538                 return false
 539         }
 540         norm = math.Sqrt((1 + norm) * (1 - norm))
 541         cos := getFloats(n, false)
 542         defer putFloats(cos)
 543         sin := getFloats(n, false)
 544         defer putFloats(sin)
 545         for i := n - 1; i >= 0; i-- {
 546                 // Compute parameters of Givens matrices that zero elements of p
 547                 // backwards.
 548                 cos[i], sin[i], norm, _ = blas64.Rotg(norm, work[i])
 549                 if norm < 0 {
 550                         norm *= -1
 551                         cos[i] *= -1
 552                         sin[i] *= -1
 553                 }
 554         }
 555         umat := c.chol.mat
 556         stride := umat.Stride
 557         for i := n - 1; i >= 0; i-- {
 558                 // Apply Givens matrices to U.
 559                 // TODO(vladimir-ch): Use workspace to avoid modifying the
 560                 // receiver in case an invalid factorization is created.
 561                 blas64.Rot(n-i, blas64.Vector{1, work[i:n]}, blas64.Vector{1, umat.Data[i*stride+i : i*stride+n]}, cos[i], sin[i])
 562                 if umat.Data[i*stride+i] == 0 {
 563                         // The matrix is singular (may rarely happen due to
 564                         // floating-point effects?).
 565                         ok = false
 566                 } else if umat.Data[i*stride+i] < 0 {
 567                         // Diagonal elements should be positive. If it happens
 568                         // that on the i-th row the diagonal is negative,
 569                         // multiply U from the left by an identity matrix that
 570                         // has -1 on the i-th row.
 571                         blas64.Scal(n-i, -1, blas64.Vector{1, umat.Data[i*stride+i : i*stride+n]})
 572                 }
 573         }
 574         if ok {
 575                 c.updateCond(-1)
 576         } else {
 577                 c.Reset()
 578         }
 579         return ok
 580 }
 581
 582 func (c *Cholesky) valid() bool {
 583         return c.chol != nil && !c.chol.IsZero()
 584 }