mat: add SolveTo and SolveVecTo to SVD
MODIFY
mat/svd.go
mat/svd.go
@@ -10,7 +10,9 @@ import (
"gonum.org/v1/gonum/lapack/lapack64"
)
-// SVD is a type for creating and using the Singular Value Decomposition (SVD)
+const badRcond = "mat: invalid rcond value"
+
+// SVD is a type for creating and using the Singular Value Decomposition
// of a matrix.
type SVD struct {
kind SVDKind
@@ -147,6 +149,26 @@ func (svd *SVD) Kind() SVDKind {
return svd.kind
}
+// Rank returns the rank of A based on the count of singular values greater than
+// rcond scaled by the largest singular value.
+// Rank will panic if the receiver does not contain a successful factorization or
+// rcond is negative.
+func (svd *SVD) Rank(rcond float64) int {
+ if rcond < 0 {
+ panic(badRcond)
+ }
+ if !svd.succFact() {
+ panic(badFact)
+ }
+ s0 := svd.s[0]
+ for i, v := range svd.s {
+ if v <= rcond*s0 {
+ return i
+ }
+ }
+ return len(svd.s)
+}
+
// Cond returns the 2-norm condition number for the factorized matrix. Cond will
// panic if the receiver does not contain a successful factorization.
func (svd *SVD) Cond() float64 {
@@ -249,3 +271,123 @@ func (svd *SVD) VTo(dst *Dense) {
}
dst.Copy(tmp.T())
}
+
+// SolveTo calculates the least-squares solution to a linear matrix equation,
+// AX = B, using a rank-truncated SVD. The result is placed in dst.
+// It returns the residuals calculated from the complete SVD. For this value to be
+// valid the factorization must have been performed with at least SVDFullU.
+// The decomposition must have been factorized computing both the U and V singular
+// vectors, and dst must be either empty or have a row dimension equal to rank.
+// The supplied rank must be a positive end-slice index into svd.Values.
+func (svd *SVD) SolveTo(dst *Dense, b Matrix, rank int) []float64 {
+ if !svd.succFact() {
+ panic(badFact)
+ }
+ if rank < 1 || len(svd.s) < rank {
+ panic("svd: rank out of range")
+ }
+ kind := svd.kind
+ if kind&SVDThinU == 0 && kind&SVDFullU == 0 {
+ panic("svd: u not computed during factorization")
+ }
+ if kind&SVDThinV == 0 && kind&SVDFullV == 0 {
+ panic("svd: v not computed during factorization")
+ }
+
+ u := Dense{
+ mat: svd.u,
+ capRows: svd.u.Rows,
+ capCols: svd.u.Cols,
+ }
+ vt := Dense{
+ mat: svd.vt,
+ capRows: svd.vt.Rows,
+ capCols: svd.vt.Cols,
+ }
+ s := svd.s[:rank]
+
+ _, bc := b.Dims()
+ c := getWorkspace(svd.u.Cols, bc, false)
+ defer putWorkspace(c)
+ c.Mul(u.T(), b)
+
+ y := getWorkspace(rank, bc, false)
+ defer putWorkspace(y)
+ y.DivElem(c.slice(0, rank, 0, bc), repVector{vec: s, cols: bc})
+ dst.Mul(vt.slice(0, rank, 0, svd.vt.Cols).T(), y)
+
+ res := make([]float64, bc)
+ if rank < svd.u.Cols {
+ c = c.slice(len(s), svd.u.Cols, 0, bc)
+ for j := range res {
+ col := c.ColView(j)
+ res[j] = Dot(col, col)
+ }
+ }
+ return res
+}
+
+type repVector struct {
+ vec []float64
+ cols int
+}
+
+func (m repVector) Dims() (r, c int) { return len(m.vec), m.cols }
+func (m repVector) At(i, j int) float64 {
+ if i < 0 || len(m.vec) <= i || j < 0 || m.cols <= j {
+ panic(ErrIndexOutOfRange.string) // Panic with string to prevent mat.Error recovery.
+ }
+ return m.vec[i]
+}
+func (m repVector) T() Matrix { return Transpose{m} }
+
+// SolveVecTo calculates the least-squares solution to a linear matrix equation,
+// Ax = b, using a rank-truncated SVD. The result is placed in dst.
+// It returns the residual calculated from the complete SVD. For this value to be
+// valid the factorization must have been performed with at least SVDFullU.
+// The decomposition must have been factorized computing both the U and V singular
+// vectors, and dst must be either empty or have length equal to rank.
+// The supplied rank must be a positive index into svd.Values.
+func (svd *SVD) SolveVecTo(dst *VecDense, b Vector, rank int) float64 {
+ if !svd.succFact() {
+ panic(badFact)
+ }
+ if rank < 1 || len(svd.s) < rank {
+ panic("svd: rank out of range")
+ }
+ kind := svd.kind
+ if kind&SVDThinU == 0 && kind&SVDFullU == 0 {
+ panic("svd: u not computed during factorization")
+ }
+ if kind&SVDThinV == 0 && kind&SVDFullV == 0 {
+ panic("svd: v not computed during factorization")
+ }
+
+ u := Dense{
+ mat: svd.u,
+ capRows: svd.u.Rows,
+ capCols: svd.u.Cols,
+ }
+ vt := Dense{
+ mat: svd.vt,
+ capRows: svd.vt.Rows,
+ capCols: svd.vt.Cols,
+ }
+ s := svd.s[:rank]
+
+ c := getWorkspaceVec(svd.u.Cols, false)
+ defer putWorkspaceVec(c)
+ c.MulVec(u.T(), b)
+
+ y := getWorkspaceVec(rank, false)
+ defer putWorkspaceVec(y)
+ y.DivElemVec(c.sliceVec(0, rank), NewVecDense(rank, s))
+ dst.MulVec(vt.slice(0, rank, 0, svd.vt.Cols).T(), y)
+
+ var res float64
+ if rank < c.Len() {
+ c = c.sliceVec(rank, c.Len())
+ res = Dot(c, c)
+ }
+ return res
+}

MODIFY
mat/vector.go
mat/vector.go
@@ -121,6 +121,10 @@ func NewVecDense(n int, data []float64) *VecDense {
// SliceVec panics with ErrIndexOutOfRange if the slice is outside the capacity
// of the receiver.
func (v *VecDense) SliceVec(i, k int) Vector {
+ return v.sliceVec(i, k)
+}
+
+func (v *VecDense) sliceVec(i, k int) *VecDense {
if i < 0 || k <= i || v.Cap() < k {
panic(ErrIndexOutOfRange)
}

mat: add SolveTo and SolveVecTo to SVD