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