1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
// Copyright ©2015 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 (
"gonum.org/v1/gonum/blas"
"gonum.org/v1/gonum/blas/blas64"
"gonum.org/v1/gonum/lapack/lapack64"
)
// Solve solves the linear least squares problem
//
// minimize over x |b - A*x|_2
//
// where A is an m×n matrix, b is a given m element vector and x is n element
// solution vector. Solve assumes that A has full rank, that is
//
// rank(A) = min(m,n)
//
// If m >= n, Solve finds the unique least squares solution of an overdetermined
// system.
//
// If m < n, there is an infinite number of solutions that satisfy b-A*x=0. In
// this case Solve finds the unique solution of an underdetermined system that
// minimizes |x|_2.
//
// Several right-hand side vectors b and solution vectors x can be handled in a
// single call. Vectors b are stored in the columns of the m×k matrix B. Vectors
// x will be stored in-place into the n×k receiver.
//
// If A does not have full rank, a Condition error is returned. See the
// documentation for Condition for more information.
func (m *Dense) Solve(a, b Matrix) error {
ar, ac := a.Dims()
br, bc := b.Dims()
if ar != br {
panic(ErrShape)
}
m.reuseAsNonZeroed(ac, bc)
// TODO(btracey): Add special cases for SymDense, etc.
aU, aTrans := untranspose(a)
bU, bTrans := untranspose(b)
switch rma := aU.(type) {
case RawTriangular:
side := blas.Left
tA := blas.NoTrans
if aTrans {
tA = blas.Trans
}
switch rm := bU.(type) {
case RawMatrixer:
if m != bU || bTrans {
if m == bU || m.checkOverlap(rm.RawMatrix()) {
tmp := getDenseWorkspace(br, bc, false)
tmp.Copy(b)
m.Copy(tmp)
putDenseWorkspace(tmp)
break
}
m.Copy(b)
}
default:
if m != bU {
m.Copy(b)
} else if bTrans {
// m and b share data so Copy cannot be used directly.
tmp := getDenseWorkspace(br, bc, false)
tmp.Copy(b)
m.Copy(tmp)
putDenseWorkspace(tmp)
}
}
rm := rma.RawTriangular()
blas64.Trsm(side, tA, 1, rm, m.mat)
work := getFloat64s(3*rm.N, false)
iwork := getInts(rm.N, false)
cond := lapack64.Trcon(CondNorm, rm, work, iwork)
putFloat64s(work)
putInts(iwork)
if cond > ConditionTolerance {
return Condition(cond)
}
return nil
}
switch {
case ar == ac:
if a == b {
// x = I.
if ar == 1 {
m.mat.Data[0] = 1
return nil
}
for i := 0; i < ar; i++ {
v := m.mat.Data[i*m.mat.Stride : i*m.mat.Stride+ac]
zero(v)
v[i] = 1
}
return nil
}
var lu LU
lu.Factorize(a)
return lu.SolveTo(m, false, b)
case ar > ac:
var qr QR
qr.Factorize(a)
return qr.SolveTo(m, false, b)
default:
var lq LQ
lq.Factorize(a)
return lq.SolveTo(m, false, b)
}
}
// SolveVec solves the linear least squares problem
//
// minimize over x |b - A*x|_2
//
// where A is an m×n matrix, b is a given m element vector and x is n element
// solution vector. Solve assumes that A has full rank, that is
//
// rank(A) = min(m,n)
//
// If m >= n, Solve finds the unique least squares solution of an overdetermined
// system.
//
// If m < n, there is an infinite number of solutions that satisfy b-A*x=0. In
// this case Solve finds the unique solution of an underdetermined system that
// minimizes |x|_2.
//
// The solution vector x will be stored in-place into the receiver.
//
// If A does not have full rank, a Condition error is returned. See the
// documentation for Condition for more information.
func (v *VecDense) SolveVec(a Matrix, b Vector) error {
if _, bc := b.Dims(); bc != 1 {
panic(ErrShape)
}
_, c := a.Dims()
// The Solve implementation is non-trivial, so rather than duplicate the code,
// instead recast the VecDenses as Dense and call the matrix code.
if rv, ok := b.(RawVectorer); ok {
bmat := rv.RawVector()
if v != b {
v.checkOverlap(bmat)
}
v.reuseAsNonZeroed(c)
m := v.asDense()
// We conditionally create bm as m when b and v are identical
// to prevent the overlap detection code from identifying m
// and bm as overlapping but not identical.
bm := m
if v != b {
b := VecDense{mat: bmat}
bm = b.asDense()
}
return m.Solve(a, bm)
}
v.reuseAsNonZeroed(c)
m := v.asDense()
return m.Solve(a, b)
}
|
