488 lines
12 KiB
Go
488 lines
12 KiB
Go
// Copyright ©2013 The Gonum Authors. All rights reserved.
|
||
// Use of this source code is governed by a BSD-style
|
||
// license that can be found in the LICENSE file.
|
||
|
||
package mat
|
||
|
||
import (
|
||
"math"
|
||
|
||
"gonum.org/v1/gonum/blas"
|
||
"gonum.org/v1/gonum/blas/blas64"
|
||
"gonum.org/v1/gonum/floats"
|
||
"gonum.org/v1/gonum/lapack"
|
||
"gonum.org/v1/gonum/lapack/lapack64"
|
||
)
|
||
|
||
const (
|
||
badSliceLength = "mat: improper slice length"
|
||
badLU = "mat: invalid LU factorization"
|
||
)
|
||
|
||
// LU is a square n×n matrix represented by its LU factorization with partial
|
||
// pivoting.
|
||
//
|
||
// The factorization has the form
|
||
//
|
||
// A = P * L * U
|
||
//
|
||
// where P is a permutation matrix, L is lower triangular with unit diagonal
|
||
// elements, and U is upper triangular.
|
||
//
|
||
// Note that this matrix representation is useful for certain operations, in
|
||
// particular for solving linear systems of equations. It is very inefficient at
|
||
// other operations, in particular At is slow.
|
||
type LU struct {
|
||
lu *Dense
|
||
swaps []int
|
||
piv []int
|
||
cond float64
|
||
ok bool // Whether A is nonsingular
|
||
}
|
||
|
||
var _ Matrix = (*LU)(nil)
|
||
|
||
// Dims returns the dimensions of the matrix A.
|
||
func (lu *LU) Dims() (r, c int) {
|
||
if lu.lu == nil {
|
||
return 0, 0
|
||
}
|
||
return lu.lu.Dims()
|
||
}
|
||
|
||
// At returns the element of A at row i, column j.
|
||
func (lu *LU) At(i, j int) float64 {
|
||
n, _ := lu.Dims()
|
||
if uint(i) >= uint(n) {
|
||
panic(ErrRowAccess)
|
||
}
|
||
if uint(j) >= uint(n) {
|
||
panic(ErrColAccess)
|
||
}
|
||
|
||
i = lu.piv[i]
|
||
var val float64
|
||
for k := 0; k < min(i, j+1); k++ {
|
||
val += lu.lu.at(i, k) * lu.lu.at(k, j)
|
||
}
|
||
if i <= j {
|
||
val += lu.lu.at(i, j)
|
||
}
|
||
return val
|
||
}
|
||
|
||
// T performs an implicit transpose by returning the receiver inside a
|
||
// Transpose.
|
||
func (lu *LU) T() Matrix {
|
||
return Transpose{lu}
|
||
}
|
||
|
||
// updateCond updates the stored condition number of the matrix. anorm is the
|
||
// norm of the original matrix. If anorm is negative it will be estimated.
|
||
func (lu *LU) updateCond(anorm float64, norm lapack.MatrixNorm) {
|
||
n := lu.lu.mat.Cols
|
||
work := getFloat64s(4*n, false)
|
||
defer putFloat64s(work)
|
||
iwork := getInts(n, false)
|
||
defer putInts(iwork)
|
||
if anorm < 0 {
|
||
// This is an approximation. By the definition of a norm,
|
||
// |AB| <= |A| |B|.
|
||
// Since A = L*U, we get for the condition number κ that
|
||
// κ(A) := |A| |A^-1| = |L*U| |A^-1| <= |L| |U| |A^-1|,
|
||
// so this will overestimate the condition number somewhat.
|
||
// The norm of the original factorized matrix cannot be stored
|
||
// because of update possibilities.
|
||
u := lu.lu.asTriDense(n, blas.NonUnit, blas.Upper)
|
||
l := lu.lu.asTriDense(n, blas.Unit, blas.Lower)
|
||
unorm := lapack64.Lantr(norm, u.mat, work)
|
||
lnorm := lapack64.Lantr(norm, l.mat, work)
|
||
anorm = unorm * lnorm
|
||
}
|
||
v := lapack64.Gecon(norm, lu.lu.mat, anorm, work, iwork)
|
||
lu.cond = 1 / v
|
||
}
|
||
|
||
// Factorize computes the LU factorization of the square matrix A and stores the
|
||
// result in the receiver. The LU decomposition will complete regardless of the
|
||
// singularity of a.
|
||
//
|
||
// The L and U matrix factors can be extracted from the factorization using the
|
||
// LTo and UTo methods. The matrix P can be extracted as a row permutation using
|
||
// the RowPivots method and applied using Dense.PermuteRows.
|
||
func (lu *LU) Factorize(a Matrix) {
|
||
lu.factorize(a, CondNorm)
|
||
}
|
||
|
||
func (lu *LU) factorize(a Matrix, norm lapack.MatrixNorm) {
|
||
m, n := a.Dims()
|
||
if m != n {
|
||
panic(ErrSquare)
|
||
}
|
||
if lu.lu == nil {
|
||
lu.lu = NewDense(n, n, nil)
|
||
} else {
|
||
lu.lu.Reset()
|
||
lu.lu.reuseAsNonZeroed(n, n)
|
||
}
|
||
lu.lu.Copy(a)
|
||
lu.swaps = useInt(lu.swaps, n)
|
||
lu.piv = useInt(lu.piv, n)
|
||
work := getFloat64s(n, false)
|
||
anorm := lapack64.Lange(norm, lu.lu.mat, work)
|
||
putFloat64s(work)
|
||
lu.ok = lapack64.Getrf(lu.lu.mat, lu.swaps)
|
||
lu.updatePivots(lu.swaps)
|
||
lu.updateCond(anorm, norm)
|
||
}
|
||
|
||
func (lu *LU) updatePivots(swaps []int) {
|
||
// Replay the sequence of row swaps in order to find the row permutation.
|
||
for i := range lu.piv {
|
||
lu.piv[i] = i
|
||
}
|
||
n, _ := lu.Dims()
|
||
for i := n - 1; i >= 0; i-- {
|
||
v := swaps[i]
|
||
lu.piv[i], lu.piv[v] = lu.piv[v], lu.piv[i]
|
||
}
|
||
}
|
||
|
||
// isValid returns whether the receiver contains a factorization.
|
||
func (lu *LU) isValid() bool {
|
||
return lu.lu != nil && !lu.lu.IsEmpty()
|
||
}
|
||
|
||
// Cond returns the condition number for the factorized matrix.
|
||
// Cond will panic if the receiver does not contain a factorization.
|
||
func (lu *LU) Cond() float64 {
|
||
if !lu.isValid() {
|
||
panic(badLU)
|
||
}
|
||
return lu.cond
|
||
}
|
||
|
||
// Reset resets the factorization so that it can be reused as the receiver of a
|
||
// dimensionally restricted operation.
|
||
func (lu *LU) Reset() {
|
||
if lu.lu != nil {
|
||
lu.lu.Reset()
|
||
}
|
||
lu.swaps = lu.swaps[:0]
|
||
lu.piv = lu.piv[:0]
|
||
}
|
||
|
||
func (lu *LU) isZero() bool {
|
||
return len(lu.swaps) == 0
|
||
}
|
||
|
||
// Det returns the determinant of the matrix that has been factorized. In many
|
||
// expressions, using LogDet will be more numerically stable.
|
||
// Det will panic if the receiver does not contain a factorization.
|
||
func (lu *LU) Det() float64 {
|
||
if !lu.ok {
|
||
return 0
|
||
}
|
||
det, sign := lu.LogDet()
|
||
return math.Exp(det) * sign
|
||
}
|
||
|
||
// LogDet returns the log of the determinant and the sign of the determinant
|
||
// for the matrix that has been factorized. Numerical stability in product and
|
||
// division expressions is generally improved by working in log space.
|
||
// LogDet will panic if the receiver does not contain a factorization.
|
||
func (lu *LU) LogDet() (det float64, sign float64) {
|
||
if !lu.isValid() {
|
||
panic(badLU)
|
||
}
|
||
|
||
_, n := lu.lu.Dims()
|
||
logDiag := getFloat64s(n, false)
|
||
defer putFloat64s(logDiag)
|
||
sign = 1.0
|
||
for i := 0; i < n; i++ {
|
||
v := lu.lu.at(i, i)
|
||
if v < 0 {
|
||
sign *= -1
|
||
}
|
||
if lu.swaps[i] != i {
|
||
sign *= -1
|
||
}
|
||
logDiag[i] = math.Log(math.Abs(v))
|
||
}
|
||
return floats.Sum(logDiag), sign
|
||
}
|
||
|
||
// RowPivots returns the row permutation that represents the permutation matrix
|
||
// P from the LU factorization
|
||
//
|
||
// A = P * L * U.
|
||
//
|
||
// If dst is nil, a new slice is allocated and returned. If dst is not nil and
|
||
// the length of dst does not equal the size of the factorized matrix, RowPivots
|
||
// will panic. RowPivots will panic if the receiver does not contain a
|
||
// factorization.
|
||
func (lu *LU) RowPivots(dst []int) []int {
|
||
if !lu.isValid() {
|
||
panic(badLU)
|
||
}
|
||
_, n := lu.lu.Dims()
|
||
if dst == nil {
|
||
dst = make([]int, n)
|
||
}
|
||
if len(dst) != n {
|
||
panic(badSliceLength)
|
||
}
|
||
copy(dst, lu.piv)
|
||
return dst
|
||
}
|
||
|
||
// Pivot returns the row pivots of the receiver.
|
||
//
|
||
// Deprecated: Use RowPivots instead.
|
||
func (lu *LU) Pivot(dst []int) []int {
|
||
return lu.RowPivots(dst)
|
||
}
|
||
|
||
// RankOne updates an LU factorization as if a rank-one update had been applied to
|
||
// the original matrix A, storing the result into the receiver. That is, if in
|
||
// the original LU decomposition P * L * U = A, in the updated decomposition
|
||
// P * L' * U' = A + alpha * x * yᵀ.
|
||
// RankOne will panic if orig does not contain a factorization.
|
||
func (lu *LU) RankOne(orig *LU, alpha float64, x, y Vector) {
|
||
if !orig.isValid() {
|
||
panic(badLU)
|
||
}
|
||
|
||
// RankOne uses algorithm a1 on page 28 of "Multiple-Rank Updates to Matrix
|
||
// Factorizations for Nonlinear Analysis and Circuit Design" by Linzhong Deng.
|
||
// http://web.stanford.edu/group/SOL/dissertations/Linzhong-Deng-thesis.pdf
|
||
_, n := orig.lu.Dims()
|
||
if r, c := x.Dims(); r != n || c != 1 {
|
||
panic(ErrShape)
|
||
}
|
||
if r, c := y.Dims(); r != n || c != 1 {
|
||
panic(ErrShape)
|
||
}
|
||
if orig != lu {
|
||
if lu.isZero() {
|
||
lu.swaps = useInt(lu.swaps, n)
|
||
lu.piv = useInt(lu.piv, n)
|
||
if lu.lu == nil {
|
||
lu.lu = NewDense(n, n, nil)
|
||
} else {
|
||
lu.lu.reuseAsNonZeroed(n, n)
|
||
}
|
||
} else if len(lu.swaps) != n {
|
||
panic(ErrShape)
|
||
}
|
||
copy(lu.swaps, orig.swaps)
|
||
lu.updatePivots(lu.swaps)
|
||
lu.lu.Copy(orig.lu)
|
||
}
|
||
|
||
xs := getFloat64s(n, false)
|
||
defer putFloat64s(xs)
|
||
ys := getFloat64s(n, false)
|
||
defer putFloat64s(ys)
|
||
for i := 0; i < n; i++ {
|
||
xs[i] = x.AtVec(i)
|
||
ys[i] = y.AtVec(i)
|
||
}
|
||
|
||
// Adjust for the pivoting in the LU factorization
|
||
for i, v := range lu.swaps {
|
||
xs[i], xs[v] = xs[v], xs[i]
|
||
}
|
||
|
||
lum := lu.lu.mat
|
||
omega := alpha
|
||
for j := 0; j < n; j++ {
|
||
ujj := lum.Data[j*lum.Stride+j]
|
||
ys[j] /= ujj
|
||
theta := 1 + xs[j]*ys[j]*omega
|
||
beta := omega * ys[j] / theta
|
||
gamma := omega * xs[j]
|
||
omega -= beta * gamma
|
||
lum.Data[j*lum.Stride+j] *= theta
|
||
for i := j + 1; i < n; i++ {
|
||
xs[i] -= lum.Data[i*lum.Stride+j] * xs[j]
|
||
tmp := ys[i]
|
||
ys[i] -= lum.Data[j*lum.Stride+i] * ys[j]
|
||
lum.Data[i*lum.Stride+j] += beta * xs[i]
|
||
lum.Data[j*lum.Stride+i] += gamma * tmp
|
||
}
|
||
}
|
||
lu.updateCond(-1, CondNorm)
|
||
}
|
||
|
||
// LTo extracts the lower triangular matrix from an LU factorization.
|
||
//
|
||
// If dst is empty, LTo will resize dst to be a lower-triangular n×n matrix.
|
||
// When dst is non-empty, LTo will panic if dst is not n×n or not Lower.
|
||
// LTo will also panic if the receiver does not contain a successful
|
||
// factorization.
|
||
func (lu *LU) LTo(dst *TriDense) *TriDense {
|
||
if !lu.isValid() {
|
||
panic(badLU)
|
||
}
|
||
|
||
_, n := lu.lu.Dims()
|
||
if dst.IsEmpty() {
|
||
dst.ReuseAsTri(n, Lower)
|
||
} else {
|
||
n2, kind := dst.Triangle()
|
||
if n != n2 {
|
||
panic(ErrShape)
|
||
}
|
||
if kind != Lower {
|
||
panic(ErrTriangle)
|
||
}
|
||
}
|
||
// Extract the lower triangular elements.
|
||
for i := 1; i < n; i++ {
|
||
copy(dst.mat.Data[i*dst.mat.Stride:i*dst.mat.Stride+i], lu.lu.mat.Data[i*lu.lu.mat.Stride:i*lu.lu.mat.Stride+i])
|
||
}
|
||
// Set ones on the diagonal.
|
||
for i := 0; i < n; i++ {
|
||
dst.mat.Data[i*dst.mat.Stride+i] = 1
|
||
}
|
||
return dst
|
||
}
|
||
|
||
// UTo extracts the upper triangular matrix from an LU factorization.
|
||
//
|
||
// If dst is empty, UTo will resize dst to be an upper-triangular n×n matrix.
|
||
// When dst is non-empty, UTo will panic if dst is not n×n or not Upper.
|
||
// UTo will also panic if the receiver does not contain a successful
|
||
// factorization.
|
||
func (lu *LU) UTo(dst *TriDense) {
|
||
if !lu.isValid() {
|
||
panic(badLU)
|
||
}
|
||
|
||
_, n := lu.lu.Dims()
|
||
if dst.IsEmpty() {
|
||
dst.ReuseAsTri(n, Upper)
|
||
} else {
|
||
n2, kind := dst.Triangle()
|
||
if n != n2 {
|
||
panic(ErrShape)
|
||
}
|
||
if kind != Upper {
|
||
panic(ErrTriangle)
|
||
}
|
||
}
|
||
// Extract the upper triangular elements.
|
||
for i := 0; i < n; i++ {
|
||
copy(dst.mat.Data[i*dst.mat.Stride+i:i*dst.mat.Stride+n], lu.lu.mat.Data[i*lu.lu.mat.Stride+i:i*lu.lu.mat.Stride+n])
|
||
}
|
||
}
|
||
|
||
// SolveTo solves a system of linear equations
|
||
//
|
||
// A * X = B if trans == false
|
||
// Aᵀ * X = B if trans == true
|
||
//
|
||
// using the LU factorization of A stored in the receiver. The solution matrix X
|
||
// is stored into dst.
|
||
//
|
||
// If A is singular or near-singular a Condition error is returned. See the
|
||
// documentation for Condition for more information. SolveTo will panic if the
|
||
// receiver does not contain a factorization.
|
||
func (lu *LU) SolveTo(dst *Dense, trans bool, b Matrix) error {
|
||
if !lu.isValid() {
|
||
panic(badLU)
|
||
}
|
||
|
||
_, n := lu.lu.Dims()
|
||
br, bc := b.Dims()
|
||
if br != n {
|
||
panic(ErrShape)
|
||
}
|
||
|
||
if !lu.ok {
|
||
return Condition(math.Inf(1))
|
||
}
|
||
|
||
dst.reuseAsNonZeroed(n, bc)
|
||
bU, _ := untranspose(b)
|
||
if dst == bU {
|
||
var restore func()
|
||
dst, restore = dst.isolatedWorkspace(bU)
|
||
defer restore()
|
||
} else if rm, ok := bU.(RawMatrixer); ok {
|
||
dst.checkOverlap(rm.RawMatrix())
|
||
}
|
||
|
||
dst.Copy(b)
|
||
t := blas.NoTrans
|
||
if trans {
|
||
t = blas.Trans
|
||
}
|
||
lapack64.Getrs(t, lu.lu.mat, dst.mat, lu.swaps)
|
||
if lu.cond > ConditionTolerance {
|
||
return Condition(lu.cond)
|
||
}
|
||
return nil
|
||
}
|
||
|
||
// SolveVecTo solves a system of linear equations
|
||
//
|
||
// A * x = b if trans == false
|
||
// Aᵀ * x = b if trans == true
|
||
//
|
||
// using the LU factorization of A stored in the receiver. The solution matrix x
|
||
// is stored into dst.
|
||
//
|
||
// If A is singular or near-singular a Condition error is returned. See the
|
||
// documentation for Condition for more information. SolveVecTo will panic if the
|
||
// receiver does not contain a factorization.
|
||
func (lu *LU) SolveVecTo(dst *VecDense, trans bool, b Vector) error {
|
||
if !lu.isValid() {
|
||
panic(badLU)
|
||
}
|
||
|
||
_, n := lu.lu.Dims()
|
||
if br, bc := b.Dims(); br != n || bc != 1 {
|
||
panic(ErrShape)
|
||
}
|
||
|
||
switch rv := b.(type) {
|
||
default:
|
||
dst.reuseAsNonZeroed(n)
|
||
return lu.SolveTo(dst.asDense(), trans, b)
|
||
case RawVectorer:
|
||
if dst != b {
|
||
dst.checkOverlap(rv.RawVector())
|
||
}
|
||
|
||
if !lu.ok {
|
||
return Condition(math.Inf(1))
|
||
}
|
||
|
||
dst.reuseAsNonZeroed(n)
|
||
var restore func()
|
||
if dst == b {
|
||
dst, restore = dst.isolatedWorkspace(b)
|
||
defer restore()
|
||
}
|
||
dst.CopyVec(b)
|
||
vMat := blas64.General{
|
||
Rows: n,
|
||
Cols: 1,
|
||
Stride: dst.mat.Inc,
|
||
Data: dst.mat.Data,
|
||
}
|
||
t := blas.NoTrans
|
||
if trans {
|
||
t = blas.Trans
|
||
}
|
||
lapack64.Getrs(t, lu.lu.mat, vMat, lu.swaps)
|
||
if lu.cond > ConditionTolerance {
|
||
return Condition(lu.cond)
|
||
}
|
||
return nil
|
||
}
|
||
}
|