r3.Mat: add Hessian
MODIFY
spatial/r3/mat.go
spatial/r3/mat.go
@@ -197,3 +197,53 @@ func (m *Mat) Det() float64 {
detc := m.At(1, 0)*m.At(2, 1) - m.At(1, 1)*m.At(2, 0)
return a*deta - b*detb + c*detc
}
+
+// Hessian sets the receiver to the Hessian matrix of the scalar field at the point p,
+// approximated using finite differences with the given step sizes.
+// The field is evaluated at points in the area surrounding p by adding
+// at most 2 components of step to p. Hessian expects the field's second partial
+// derivatives are all continuous for correct results.
+func (m *Mat) Hessian(p, step Vec, field func(Vec) float64) {
+ dx := Vec{X: step.X}
+ dy := Vec{Y: step.Y}
+ dz := Vec{Z: step.Z}
+ fp := field(p)
+ fxp := field(Add(p, dx))
+ fxm := field(Sub(p, dx))
+ fxx := (fxp - 2*fp + fxm) / (step.X * step.X)
+
+ fyp := field(Add(p, dy))
+ fym := field(Sub(p, dy))
+ fyy := (fyp - 2*fp + fym) / (step.Y * step.Y)
+
+ aux := Add(dx, dy)
+ fxyp := field(Add(p, aux))
+ fxym := field(Sub(p, aux))
+ fxy := (fxyp - fxp - fyp + 2*fp - fxm - fym + fxym) / (2 * step.X * step.Y)
+
+ fzp := field(Add(p, dz))
+ fzm := field(Sub(p, dz))
+ fzz := (fzp - 2*fp + fzm) / (step.Z * step.Z)
+
+ aux = Add(dx, dz)
+ fxzp := field(Add(p, aux))
+ fxzm := field(Sub(p, aux))
+ fxz := (fxzp - fxp - fzp + 2*fp - fxm - fzm + fxzm) / (2 * step.X * step.Z)
+
+ aux = Add(dy, dz)
+ fyzp := field(Add(p, aux))
+ fyzm := field(Sub(p, aux))
+ fyz := (fyzp - fyp - fzp + 2*fp - fym - fzm + fyzm) / (2 * step.Y * step.Z)
+
+ m.Set(0, 0, fxx)
+ m.Set(0, 1, fxy)
+ m.Set(0, 2, fxz)
+
+ m.Set(1, 0, fxy)
+ m.Set(1, 1, fyy)
+ m.Set(1, 2, fyz)
+
+ m.Set(2, 0, fxz)
+ m.Set(2, 1, fyz)
+ m.Set(2, 2, fzz)
+}

MODIFY
spatial/r3/mat_test.go
spatial/r3/mat_test.go
@@ -230,3 +230,63 @@ func BenchmarkQuat(b *testing.B) {
}
}
}
+
+var scalarFields = []struct {
+ // This is the scalar field function.
+ field func(p Vec) float64
+ hessian func(p Vec) *Mat
+}{
+ {
+ field: func(p Vec) float64 {
+ // 4*x^3 + 5*y^2 + 3*z^4
+ z2 := p.Z * p.Z
+ return 4*p.X*p.X*p.X + 5*p.Y*p.Y + 3*z2*z2
+ },
+ hessian: func(p Vec) *Mat {
+ return NewMat([]float64{
+ 24 * p.X, 0, 0,
+ 0, 10, 0,
+ 0, 0, 36 * p.Z * p.Z,
+ })
+ },
+ },
+ {
+ field: func(p Vec) float64 {
+ // cos(x) * sin(z) * y^4
+ y2 := p.Y * p.Y
+ return math.Cos(p.X) * math.Sin(p.Z) * y2 * y2
+ },
+ hessian: func(p Vec) *Mat {
+ y3 := p.Y * p.Y * p.Y
+ y4 := y3 * p.Y
+ sx, cx := math.Sincos(p.X)
+ sz, cz := math.Sincos(p.Z)
+ return NewMat([]float64{
+ -y4 * cx * sz, -4 * y3 * sx * sz, -y4 * sx * cz,
+ -4 * y3 * sx * sz, 12 * p.Y * p.Y * cx * sz, 4 * y3 * cx * cz,
+ -y4 * sx * cz, 4 * y3 * cx * cz, -y4 * cx * sz,
+ })
+ },
+ },
+}
+
+func TestMatHessian(t *testing.T) {
+ const (
+ tol = 1e-5
+ h = 8e-4
+ )
+ step := Vec{X: h, Y: h, Z: h}
+ rnd := rand.New(rand.NewSource(1))
+ for _, test := range scalarFields {
+ for i := 0; i < 30; i++ {
+ p := randomVec(rnd)
+ got := NewMat(nil)
+ got.Hessian(p, step, test.field)
+ want := test.hessian(p)
+ if !mat.EqualApprox(got, want, tol) {
+ t.Errorf("matrices not equal within tol\ngot: %v\nwant: %v",
+ mat.Formatted(got), mat.Formatted(want))
+ }
+ }
+ }
+}

r3.Mat: add Hessian