-
Notifications
You must be signed in to change notification settings - Fork 4
/
precond.go
56 lines (49 loc) · 1.2 KB
/
precond.go
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
package rbfscale
import (
"math"
"github.com/unixpickle/num-analysis/linalg"
"github.com/unixpickle/num-analysis/linalg/svd"
)
const precondRadius = 3
// A preconditioner preconditions conjugate gradients for
// solving the interpolation system.
type preconditioner struct {
row linalg.Vector
dim int
}
func newPreconditioner(width, height int, variance float64) *preconditioner {
n := precondRadius*2 + 1
mat := linalg.NewMatrix(n, n)
for i := 0; i < n; i++ {
for j := 0; j < n; j++ {
mag := math.Exp(-math.Pow(float64(i-j), 2) / (2 * variance))
mat.Set(i, j, mag)
}
}
v, d, u := svd.Decompose(mat)
for i := 0; i < n; i++ {
d.Set(i, i, 1/math.Sqrt(float64(d.Get(i, i))))
}
inv := v.Mul(d).Mul(u)
return &preconditioner{
row: inv.Data[precondRadius*inv.Cols : (precondRadius+1)*inv.Cols],
dim: width * height,
}
}
func (p *preconditioner) Dim() int {
return p.dim
}
func (p *preconditioner) Apply(in linalg.Vector) linalg.Vector {
res := make(linalg.Vector, len(in))
for i := range in {
offset := (len(p.row) - 1) / 2
for j := i - offset; j <= i+offset && j < p.dim; j++ {
if j < 0 {
continue
}
coeff := p.row[j-i+offset]
res[i] += in[j] * coeff
}
}
return res
}