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.
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"
17 const badSliceLength = "mat: improper slice length"
19 // LU is a type for creating and using the LU factorization of a matrix.
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) {
30 work := getFloats(4*n, false)
32 iwork := getInts(n, false)
35 // This is an approximation. By the definition of a norm,
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)
48 v := lapack64.Gecon(norm, lu.lu.mat, anorm, work, iwork)
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.
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)
63 func (lu *LU) factorize(a Matrix, norm lapack.MatrixNorm) {
69 lu.lu = NewDense(r, r, nil)
75 if cap(lu.pivot) < r {
76 lu.pivot = make([]int, r)
78 lu.pivot = lu.pivot[:r]
79 work := getFloats(r, false)
80 anorm := lapack64.Lange(norm, lu.lu.mat, work)
82 lapack64.Getrf(lu.lu.mat, lu.pivot)
83 lu.updateCond(anorm, norm)
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")
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() {
101 lu.pivot = lu.pivot[:0]
104 func (lu *LU) isZero() bool {
105 return len(lu.pivot) == 0
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
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) {
120 logDiag := getFloats(n, false)
121 defer putFloats(logDiag)
123 for i := 0; i < n; i++ {
128 if lu.pivot[i] != i {
131 logDiag[i] = math.Log(math.Abs(v))
133 return floats.Sum(logDiag), sign
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 {
143 swaps = make([]int, n)
146 panic(badSliceLength)
148 // Perform the inverse of the row swaps in order to find the final
149 // row swap position.
150 for i := range swaps {
153 for i := n - 1; i >= 0; i-- {
155 swaps[i], swaps[v] = swaps[v], swaps[i]
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 {
172 if r, c := y.Dims(); r != n || c != 1 {
177 if cap(lu.pivot) < n {
178 lu.pivot = make([]int, n)
180 lu.pivot = lu.pivot[:n]
182 lu.lu = NewDense(n, n, nil)
186 } else if len(lu.pivot) != n {
189 copy(lu.pivot, orig.pivot)
193 xs := getFloats(n, false)
195 ys := getFloats(n, false)
197 for i := 0; i < n; i++ {
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]
209 for j := 0; j < n; j++ {
210 ujj := lum.Data[j*lum.Stride+j]
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]
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
225 lu.updateCond(-1, CondNorm)
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 {
233 dst = NewTriDense(n, Lower, nil)
235 dst.reuseAs(n, Lower)
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]
243 // Set ones on the diagonal.
244 for i := 0; i < n; i++ {
245 dst.mat.Data[i*dst.mat.Stride+i] = 1
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 {
255 dst = NewTriDense(n, Upper, nil)
257 dst.reuseAs(n, Upper)
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]
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) {
274 for i := 0; i < r; i++ {
275 zero(m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+r])
280 m.mat.Data[i*m.mat.Stride+v] = 1
284 // Solve solves a system of linear equations using the LU decomposition of a matrix.
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
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 {
299 // TODO(btracey): Should test the condition number instead of testing that
300 // the determinant is exactly zero.
302 return Condition(math.Inf(1))
306 bU, _ := untranspose(b)
309 m, restore = m.isolatedWorkspace(bU)
311 } else if rm, ok := bU.(RawMatrixer); ok {
312 m.checkOverlap(rm.RawMatrix())
320 lapack64.Getrs(t, lu.lu.mat, m.mat, lu.pivot)
321 if lu.cond > ConditionTolerance {
322 return Condition(lu.cond)
327 // SolveVec solves a system of linear equations using the LU decomposition of a matrix.
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
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 {
338 if br, bc := b.Dims(); br != n || bc != 1 {
341 switch rv := b.(type) {
344 return lu.Solve(v.asDense(), trans, b)
347 v.checkOverlap(rv.RawVector())
349 // TODO(btracey): Should test the condition number instead of testing that
350 // the determinant is exactly zero.
352 return Condition(math.Inf(1))
358 v, restore = v.isolatedWorkspace(b)
362 vMat := blas64.General{
372 lapack64.Getrs(t, lu.lu.mat, vMat, lu.pivot)
373 if lu.cond > ConditionTolerance {
374 return Condition(lu.cond)