vendor/gonum.org/v1/gonum/mat/qr.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"
  13         "gonum.org/v1/gonum/lapack/lapack64"
  14 )
  15
  16 // QR is a type for creating and using the QR factorization of a matrix.
  17 type QR struct {
  18         qr   *Dense
  19         tau  []float64
  20         cond float64
  21 }
  22
  23 func (qr *QR) updateCond(norm lapack.MatrixNorm) {
  24         // Since A = Q*R, and Q is orthogonal, we get for the condition number κ
  25         //  κ(A) := |A| |A^-1| = |Q*R| |(Q*R)^-1| = |R| |R^-1 * Q^T|
  26         //        = |R| |R^-1| = κ(R),
  27         // where we used that fact that Q^-1 = Q^T. However, this assumes that
  28         // the matrix norm is invariant under orthogonal transformations which
  29         // is not the case for CondNorm. Hopefully the error is negligible: κ
  30         // is only a qualitative measure anyway.
  31         n := qr.qr.mat.Cols
  32         work := getFloats(3*n, false)
  33         iwork := getInts(n, false)
  34         r := qr.qr.asTriDense(n, blas.NonUnit, blas.Upper)
  35         v := lapack64.Trcon(norm, r.mat, work, iwork)
  36         putFloats(work)
  37         putInts(iwork)
  38         qr.cond = 1 / v
  39 }
  40
  41 // Factorize computes the QR factorization of an m×n matrix a where m >= n. The QR
  42 // factorization always exists even if A is singular.
  43 //
  44 // The QR decomposition is a factorization of the matrix A such that A = Q * R.
  45 // The matrix Q is an orthonormal m×m matrix, and R is an m×n upper triangular matrix.
  46 // Q and R can be extracted using the QTo and RTo methods.
  47 func (qr *QR) Factorize(a Matrix) {
  48         qr.factorize(a, CondNorm)
  49 }
  50
  51 func (qr *QR) factorize(a Matrix, norm lapack.MatrixNorm) {
  52         m, n := a.Dims()
  53         if m < n {
  54                 panic(ErrShape)
  55         }
  56         k := min(m, n)
  57         if qr.qr == nil {
  58                 qr.qr = &Dense{}
  59         }
  60         qr.qr.Clone(a)
  61         work := []float64{0}
  62         qr.tau = make([]float64, k)
  63         lapack64.Geqrf(qr.qr.mat, qr.tau, work, -1)
  64
  65         work = getFloats(int(work[0]), false)
  66         lapack64.Geqrf(qr.qr.mat, qr.tau, work, len(work))
  67         putFloats(work)
  68         qr.updateCond(norm)
  69 }
  70
  71 // Cond returns the condition number for the factorized matrix.
  72 // Cond will panic if the receiver does not contain a successful factorization.
  73 func (qr *QR) Cond() float64 {
  74         if qr.qr == nil || qr.qr.IsZero() {
  75                 panic("qr: no decomposition computed")
  76         }
  77         return qr.cond
  78 }
  79
  80 // TODO(btracey): Add in the "Reduced" forms for extracting the n×n orthogonal
  81 // and upper triangular matrices.
  82
  83 // RTo extracts the m×n upper trapezoidal matrix from a QR decomposition.
  84 // If dst is nil, a new matrix is allocated. The resulting dst matrix is returned.
  85 func (qr *QR) RTo(dst *Dense) *Dense {
  86         r, c := qr.qr.Dims()
  87         if dst == nil {
  88                 dst = NewDense(r, c, nil)
  89         } else {
  90                 dst.reuseAs(r, c)
  91         }
  92
  93         // Disguise the QR as an upper triangular
  94         t := &TriDense{
  95                 mat: blas64.Triangular{
  96                         N:      c,
  97                         Stride: qr.qr.mat.Stride,
  98                         Data:   qr.qr.mat.Data,
  99                         Uplo:   blas.Upper,
 100                         Diag:   blas.NonUnit,
 101                 },
 102                 cap: qr.qr.capCols,
 103         }
 104         dst.Copy(t)
 105
 106         // Zero below the triangular.
 107         for i := r; i < c; i++ {
 108                 zero(dst.mat.Data[i*dst.mat.Stride : i*dst.mat.Stride+c])
 109         }
 110
 111         return dst
 112 }
 113
 114 // QTo extracts the m×m orthonormal matrix Q from a QR decomposition.
 115 // If dst is nil, a new matrix is allocated. The resulting Q matrix is returned.
 116 func (qr *QR) QTo(dst *Dense) *Dense {
 117         r, _ := qr.qr.Dims()
 118         if dst == nil {
 119                 dst = NewDense(r, r, nil)
 120         } else {
 121                 dst.reuseAsZeroed(r, r)
 122         }
 123
 124         // Set Q = I.
 125         for i := 0; i < r*r; i += r + 1 {
 126                 dst.mat.Data[i] = 1
 127         }
 128
 129         // Construct Q from the elementary reflectors.
 130         work := []float64{0}
 131         lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, -1)
 132         work = getFloats(int(work[0]), false)
 133         lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, dst.mat, work, len(work))
 134         putFloats(work)
 135
 136         return dst
 137 }
 138
 139 // Solve finds a minimum-norm solution to a system of linear equations defined
 140 // by the matrices A and b, where A is an m×n matrix represented in its QR factorized
 141 // form. If A is singular or near-singular a Condition error is returned.
 142 // See the documentation for Condition for more information.
 143 //
 144 // The minimization problem solved depends on the input parameters.
 145 //  If trans == false, find X such that ||A*X - b||_2 is minimized.
 146 //  If trans == true, find the minimum norm solution of A^T * X = b.
 147 // The solution matrix, X, is stored in place into m.
 148 func (qr *QR) Solve(m *Dense, trans bool, b Matrix) error {
 149         r, c := qr.qr.Dims()
 150         br, bc := b.Dims()
 151
 152         // The QR solve algorithm stores the result in-place into the right hand side.
 153         // The storage for the answer must be large enough to hold both b and x.
 154         // However, this method's receiver must be the size of x. Copy b, and then
 155         // copy the result into m at the end.
 156         if trans {
 157                 if c != br {
 158                         panic(ErrShape)
 159                 }
 160                 m.reuseAs(r, bc)
 161         } else {
 162                 if r != br {
 163                         panic(ErrShape)
 164                 }
 165                 m.reuseAs(c, bc)
 166         }
 167         // Do not need to worry about overlap between m and b because x has its own
 168         // independent storage.
 169         x := getWorkspace(max(r, c), bc, false)
 170         x.Copy(b)
 171         t := qr.qr.asTriDense(qr.qr.mat.Cols, blas.NonUnit, blas.Upper).mat
 172         if trans {
 173                 ok := lapack64.Trtrs(blas.Trans, t, x.mat)
 174                 if !ok {
 175                         return Condition(math.Inf(1))
 176                 }
 177                 for i := c; i < r; i++ {
 178                         zero(x.mat.Data[i*x.mat.Stride : i*x.mat.Stride+bc])
 179                 }
 180                 work := []float64{0}
 181                 lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, x.mat, work, -1)
 182                 work = getFloats(int(work[0]), false)
 183                 lapack64.Ormqr(blas.Left, blas.NoTrans, qr.qr.mat, qr.tau, x.mat, work, len(work))
 184                 putFloats(work)
 185         } else {
 186                 work := []float64{0}
 187                 lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, x.mat, work, -1)
 188                 work = getFloats(int(work[0]), false)
 189                 lapack64.Ormqr(blas.Left, blas.Trans, qr.qr.mat, qr.tau, x.mat, work, len(work))
 190                 putFloats(work)
 191
 192                 ok := lapack64.Trtrs(blas.NoTrans, t, x.mat)
 193                 if !ok {
 194                         return Condition(math.Inf(1))
 195                 }
 196         }
 197         // M was set above to be the correct size for the result.
 198         m.Copy(x)
 199         putWorkspace(x)
 200         if qr.cond > ConditionTolerance {
 201                 return Condition(qr.cond)
 202         }
 203         return nil
 204 }
 205
 206 // SolveVec finds a minimum-norm solution to a system of linear equations,
 207 //  Ax = b.
 208 // See QR.Solve for the full documentation.
 209 func (qr *QR) SolveVec(v *VecDense, trans bool, b Vector) error {
 210         r, c := qr.qr.Dims()
 211         if _, bc := b.Dims(); bc != 1 {
 212                 panic(ErrShape)
 213         }
 214
 215         // The Solve implementation is non-trivial, so rather than duplicate the code,
 216         // instead recast the VecDenses as Dense and call the matrix code.
 217         bm := Matrix(b)
 218         if rv, ok := b.(RawVectorer); ok {
 219                 bmat := rv.RawVector()
 220                 if v != b {
 221                         v.checkOverlap(bmat)
 222                 }
 223                 b := VecDense{mat: bmat, n: b.Len()}
 224                 bm = b.asDense()
 225         }
 226         if trans {
 227                 v.reuseAs(r)
 228         } else {
 229                 v.reuseAs(c)
 230         }
 231         return qr.Solve(v.asDense(), trans, bm)
 232
 233 }