vendor/gonum.org/v1/gonum/mat/lu.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/floats"
  13         "gonum.org/v1/gonum/lapack"
  14         "gonum.org/v1/gonum/lapack/lapack64"
  15 )
  16
  17 const badSliceLength = "mat: improper slice length"
  18
  19 // LU is a type for creating and using the LU factorization of a matrix.
  20 type LU struct {
  21         lu    *Dense
  22         pivot []int
  23         cond  float64
  24 }
  25
  26 // updateCond updates the stored condition number of the matrix. anorm is the
  27 // norm of the original matrix. If anorm is negative it will be estimated.
  28 func (lu *LU) updateCond(anorm float64, norm lapack.MatrixNorm) {
  29         n := lu.lu.mat.Cols
  30         work := getFloats(4*n, false)
  31         defer putFloats(work)
  32         iwork := getInts(n, false)
  33         defer putInts(iwork)
  34         if anorm < 0 {
  35                 // This is an approximation. By the definition of a norm,
  36                 //  |AB| <= |A| |B|.
  37                 // Since A = L*U, we get for the condition number κ that
  38                 //  κ(A) := |A| |A^-1| = |L*U| |A^-1| <= |L| |U| |A^-1|,
  39                 // so this will overestimate the condition number somewhat.
  40                 // The norm of the original factorized matrix cannot be stored
  41                 // because of update possibilities.
  42                 u := lu.lu.asTriDense(n, blas.NonUnit, blas.Upper)
  43                 l := lu.lu.asTriDense(n, blas.Unit, blas.Lower)
  44                 unorm := lapack64.Lantr(norm, u.mat, work)
  45                 lnorm := lapack64.Lantr(norm, l.mat, work)
  46                 anorm = unorm * lnorm
  47         }
  48         v := lapack64.Gecon(norm, lu.lu.mat, anorm, work, iwork)
  49         lu.cond = 1 / v
  50 }
  51
  52 // Factorize computes the LU factorization of the square matrix a and stores the
  53 // result. The LU decomposition will complete regardless of the singularity of a.
  54 //
  55 // The LU factorization is computed with pivoting, and so really the decomposition
  56 // is a PLU decomposition where P is a permutation matrix. The individual matrix
  57 // factors can be extracted from the factorization using the Permutation method
  58 // on Dense, and the LU LTo and UTo methods.
  59 func (lu *LU) Factorize(a Matrix) {
  60         lu.factorize(a, CondNorm)
  61 }
  62
  63 func (lu *LU) factorize(a Matrix, norm lapack.MatrixNorm) {
  64         r, c := a.Dims()
  65         if r != c {
  66                 panic(ErrSquare)
  67         }
  68         if lu.lu == nil {
  69                 lu.lu = NewDense(r, r, nil)
  70         } else {
  71                 lu.lu.Reset()
  72                 lu.lu.reuseAs(r, r)
  73         }
  74         lu.lu.Copy(a)
  75         if cap(lu.pivot) < r {
  76                 lu.pivot = make([]int, r)
  77         }
  78         lu.pivot = lu.pivot[:r]
  79         work := getFloats(r, false)
  80         anorm := lapack64.Lange(norm, lu.lu.mat, work)
  81         putFloats(work)
  82         lapack64.Getrf(lu.lu.mat, lu.pivot)
  83         lu.updateCond(anorm, norm)
  84 }
  85
  86 // Cond returns the condition number for the factorized matrix.
  87 // Cond will panic if the receiver does not contain a successful factorization.
  88 func (lu *LU) Cond() float64 {
  89         if lu.lu == nil || lu.lu.IsZero() {
  90                 panic("lu: no decomposition computed")
  91         }
  92         return lu.cond
  93 }
  94
  95 // Reset resets the factorization so that it can be reused as the receiver of a
  96 // dimensionally restricted operation.
  97 func (lu *LU) Reset() {
  98         if lu.lu != nil {
  99                 lu.lu.Reset()
 100         }
 101         lu.pivot = lu.pivot[:0]
 102 }
 103
 104 func (lu *LU) isZero() bool {
 105         return len(lu.pivot) == 0
 106 }
 107
 108 // Det returns the determinant of the matrix that has been factorized. In many
 109 // expressions, using LogDet will be more numerically stable.
 110 func (lu *LU) Det() float64 {
 111         det, sign := lu.LogDet()
 112         return math.Exp(det) * sign
 113 }
 114
 115 // LogDet returns the log of the determinant and the sign of the determinant
 116 // for the matrix that has been factorized. Numerical stability in product and
 117 // division expressions is generally improved by working in log space.
 118 func (lu *LU) LogDet() (det float64, sign float64) {
 119         _, n := lu.lu.Dims()
 120         logDiag := getFloats(n, false)
 121         defer putFloats(logDiag)
 122         sign = 1.0
 123         for i := 0; i < n; i++ {
 124                 v := lu.lu.at(i, i)
 125                 if v < 0 {
 126                         sign *= -1
 127                 }
 128                 if lu.pivot[i] != i {
 129                         sign *= -1
 130                 }
 131                 logDiag[i] = math.Log(math.Abs(v))
 132         }
 133         return floats.Sum(logDiag), sign
 134 }
 135
 136 // Pivot returns pivot indices that enable the construction of the permutation
 137 // matrix P (see Dense.Permutation). If swaps == nil, then new memory will be
 138 // allocated, otherwise the length of the input must be equal to the size of the
 139 // factorized matrix.
 140 func (lu *LU) Pivot(swaps []int) []int {
 141         _, n := lu.lu.Dims()
 142         if swaps == nil {
 143                 swaps = make([]int, n)
 144         }
 145         if len(swaps) != n {
 146                 panic(badSliceLength)
 147         }
 148         // Perform the inverse of the row swaps in order to find the final
 149         // row swap position.
 150         for i := range swaps {
 151                 swaps[i] = i
 152         }
 153         for i := n - 1; i >= 0; i-- {
 154                 v := lu.pivot[i]
 155                 swaps[i], swaps[v] = swaps[v], swaps[i]
 156         }
 157         return swaps
 158 }
 159
 160 // RankOne updates an LU factorization as if a rank-one update had been applied to
 161 // the original matrix A, storing the result into the receiver. That is, if in
 162 // the original LU decomposition P * L * U = A, in the updated decomposition
 163 // P * L * U = A + alpha * x * y^T.
 164 func (lu *LU) RankOne(orig *LU, alpha float64, x, y Vector) {
 165         // RankOne uses algorithm a1 on page 28 of "Multiple-Rank Updates to Matrix
 166         // Factorizations for Nonlinear Analysis and Circuit Design" by Linzhong Deng.
 167         // http://web.stanford.edu/group/SOL/dissertations/Linzhong-Deng-thesis.pdf
 168         _, n := orig.lu.Dims()
 169         if r, c := x.Dims(); r != n || c != 1 {
 170                 panic(ErrShape)
 171         }
 172         if r, c := y.Dims(); r != n || c != 1 {
 173                 panic(ErrShape)
 174         }
 175         if orig != lu {
 176                 if lu.isZero() {
 177                         if cap(lu.pivot) < n {
 178                                 lu.pivot = make([]int, n)
 179                         }
 180                         lu.pivot = lu.pivot[:n]
 181                         if lu.lu == nil {
 182                                 lu.lu = NewDense(n, n, nil)
 183                         } else {
 184                                 lu.lu.reuseAs(n, n)
 185                         }
 186                 } else if len(lu.pivot) != n {
 187                         panic(ErrShape)
 188                 }
 189                 copy(lu.pivot, orig.pivot)
 190                 lu.lu.Copy(orig.lu)
 191         }
 192
 193         xs := getFloats(n, false)
 194         defer putFloats(xs)
 195         ys := getFloats(n, false)
 196         defer putFloats(ys)
 197         for i := 0; i < n; i++ {
 198                 xs[i] = x.AtVec(i)
 199                 ys[i] = y.AtVec(i)
 200         }
 201
 202         // Adjust for the pivoting in the LU factorization
 203         for i, v := range lu.pivot {
 204                 xs[i], xs[v] = xs[v], xs[i]
 205         }
 206
 207         lum := lu.lu.mat
 208         omega := alpha
 209         for j := 0; j < n; j++ {
 210                 ujj := lum.Data[j*lum.Stride+j]
 211                 ys[j] /= ujj
 212                 theta := 1 + xs[j]*ys[j]*omega
 213                 beta := omega * ys[j] / theta
 214                 gamma := omega * xs[j]
 215                 omega -= beta * gamma
 216                 lum.Data[j*lum.Stride+j] *= theta
 217                 for i := j + 1; i < n; i++ {
 218                         xs[i] -= lum.Data[i*lum.Stride+j] * xs[j]
 219                         tmp := ys[i]
 220                         ys[i] -= lum.Data[j*lum.Stride+i] * ys[j]
 221                         lum.Data[i*lum.Stride+j] += beta * xs[i]
 222                         lum.Data[j*lum.Stride+i] += gamma * tmp
 223                 }
 224         }
 225         lu.updateCond(-1, CondNorm)
 226 }
 227
 228 // LTo extracts the lower triangular matrix from an LU factorization.
 229 // If dst is nil, a new matrix is allocated. The resulting L matrix is returned.
 230 func (lu *LU) LTo(dst *TriDense) *TriDense {
 231         _, n := lu.lu.Dims()
 232         if dst == nil {
 233                 dst = NewTriDense(n, Lower, nil)
 234         } else {
 235                 dst.reuseAs(n, Lower)
 236         }
 237         // Extract the lower triangular elements.
 238         for i := 0; i < n; i++ {
 239                 for j := 0; j < i; j++ {
 240                         dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j]
 241                 }
 242         }
 243         // Set ones on the diagonal.
 244         for i := 0; i < n; i++ {
 245                 dst.mat.Data[i*dst.mat.Stride+i] = 1
 246         }
 247         return dst
 248 }
 249
 250 // UTo extracts the upper triangular matrix from an LU factorization.
 251 // If dst is nil, a new matrix is allocated. The resulting U matrix is returned.
 252 func (lu *LU) UTo(dst *TriDense) *TriDense {
 253         _, n := lu.lu.Dims()
 254         if dst == nil {
 255                 dst = NewTriDense(n, Upper, nil)
 256         } else {
 257                 dst.reuseAs(n, Upper)
 258         }
 259         // Extract the upper triangular elements.
 260         for i := 0; i < n; i++ {
 261                 for j := i; j < n; j++ {
 262                         dst.mat.Data[i*dst.mat.Stride+j] = lu.lu.mat.Data[i*lu.lu.mat.Stride+j]
 263                 }
 264         }
 265         return dst
 266 }
 267
 268 // Permutation constructs an r×r permutation matrix with the given row swaps.
 269 // A permutation matrix has exactly one element equal to one in each row and column
 270 // and all other elements equal to zero. swaps[i] specifies the row with which
 271 // i will be swapped, which is equivalent to the non-zero column of row i.
 272 func (m *Dense) Permutation(r int, swaps []int) {
 273         m.reuseAs(r, r)
 274         for i := 0; i < r; i++ {
 275                 zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+r])
 276                 v := swaps[i]
 277                 if v < 0 || v >= r {
 278                         panic(ErrRowAccess)
 279                 }
 280                 m.mat.Data[i*m.mat.Stride+v] = 1
 281         }
 282 }
 283
 284 // Solve solves a system of linear equations using the LU decomposition of a matrix.
 285 // It computes
 286 //  A * x = b if trans == false
 287 //  A^T * x = b if trans == true
 288 // In both cases, A is represented in LU factorized form, and the matrix x is
 289 // stored into m.
 290 //
 291 // If A is singular or near-singular a Condition error is returned. See
 292 // the documentation for Condition for more information.
 293 func (lu *LU) Solve(m *Dense, trans bool, b Matrix) error {
 294         _, n := lu.lu.Dims()
 295         br, bc := b.Dims()
 296         if br != n {
 297                 panic(ErrShape)
 298         }
 299         // TODO(btracey): Should test the condition number instead of testing that
 300         // the determinant is exactly zero.
 301         if lu.Det() == 0 {
 302                 return Condition(math.Inf(1))
 303         }
 304
 305         m.reuseAs(n, bc)
 306         bU, _ := untranspose(b)
 307         var restore func()
 308         if m == bU {
 309                 m, restore = m.isolatedWorkspace(bU)
 310                 defer restore()
 311         } else if rm, ok := bU.(RawMatrixer); ok {
 312                 m.checkOverlap(rm.RawMatrix())
 313         }
 314
 315         m.Copy(b)
 316         t := blas.NoTrans
 317         if trans {
 318                 t = blas.Trans
 319         }
 320         lapack64.Getrs(t, lu.lu.mat, m.mat, lu.pivot)
 321         if lu.cond > ConditionTolerance {
 322                 return Condition(lu.cond)
 323         }
 324         return nil
 325 }
 326
 327 // SolveVec solves a system of linear equations using the LU decomposition of a matrix.
 328 // It computes
 329 //  A * x = b if trans == false
 330 //  A^T * x = b if trans == true
 331 // In both cases, A is represented in LU factorized form, and the matrix x is
 332 // stored into v.
 333 //
 334 // If A is singular or near-singular a Condition error is returned. See
 335 // the documentation for Condition for more information.
 336 func (lu *LU) SolveVec(v *VecDense, trans bool, b Vector) error {
 337         _, n := lu.lu.Dims()
 338         if br, bc := b.Dims(); br != n || bc != 1 {
 339                 panic(ErrShape)
 340         }
 341         switch rv := b.(type) {
 342         default:
 343                 v.reuseAs(n)
 344                 return lu.Solve(v.asDense(), trans, b)
 345         case RawVectorer:
 346                 if v != b {
 347                         v.checkOverlap(rv.RawVector())
 348                 }
 349                 // TODO(btracey): Should test the condition number instead of testing that
 350                 // the determinant is exactly zero.
 351                 if lu.Det() == 0 {
 352                         return Condition(math.Inf(1))
 353                 }
 354
 355                 v.reuseAs(n)
 356                 var restore func()
 357                 if v == b {
 358                         v, restore = v.isolatedWorkspace(b)
 359                         defer restore()
 360                 }
 361                 v.CopyVec(b)
 362                 vMat := blas64.General{
 363                         Rows:   n,
 364                         Cols:   1,
 365                         Stride: v.mat.Inc,
 366                         Data:   v.mat.Data,
 367                 }
 368                 t := blas.NoTrans
 369                 if trans {
 370                         t = blas.Trans
 371                 }
 372                 lapack64.Getrs(t, lu.lu.mat, vMat, lu.pivot)
 373                 if lu.cond > ConditionTolerance {
 374                         return Condition(lu.cond)
 375                 }
 376                 return nil
 377         }
 378 }