Pbenner Autodiff Save

Documentation

Autodiff is a numerical optimization and linear algebra library for the Go / Golang programming language. It implements basic automatic differentation for many mathematical routines. The documentation of this package can be found here.

Scalars

Autodiff defines three different scalar types. A Scalar contains a single mutable value that can be the result of a mathematical operation, whereas the value of a ConstScalar is constant and fixed when the scalar is created. Automatic differentiation is implemented by MagicScalar types that allow to compute first and second order derivatives. Autodiff supports the following scalars types:

The ConstScalar, Scalar and MagicScalar interfaces define the following operations:

The Scalar and MagicScalar interfaces define the following operations:

The Scalar and MagicScalar interfaces define the following mathematical operations:

Vectors and Matrices

Autodiff implements dense and sparse vectors and matrices that support basic linear algebra operations. The following vector and matrix types are provided by autodiff:

Autodiff defines three vector interfaces ConstVector, Vector, and MagicVector:

Vectors support the following mathematical operations:

Autodiff defines three matrix interfaces ConstMatrix, Matrix, and MagicMatrix:

Matrices support the following linear algebra operations:

Methods, such as VaddV and MaddM, are generic and accept vector or matrix types that implement the respective ConstVector or ConstMatrix interface. However, opertions on interface types are much slower than on concrete types, which is why most vector and matrix types in autodiff also implement methods that operate on concrete types. For instance, DenseFloat64Vector implements a method called VADDV that takes as arguments two objects of type DenseFloat64Vector. Methods that operate on concrete types are always named in capital letters.

Algorithms

The algorithms package contains more complex linear algebra and optimization routines:

Basic usage

Import the autodiff library with

  import . "github.com/pbenner/autodiff"

A scalar holding the value 1.0 can be defined in several ways, i.e.

  a := NullScalar(Real64Type)
  a.SetFloat64(1.0)
  b := NewReal64(1.0)
  c := NewFloat64(1.0)

a and b are both MagicScalars, however a has type Scalar whereas b has type *Real64 which implements the Scalar interface. Variable c is of type Float64 which cannot carry any derivatives. Basic operations such as additions are defined on all Scalars, i.e.

  a.Add(a, b)

which stores the result of adding a and b in a. The ConstFloat64 type allows to define float64 constants without allocation of additional memory. For instance

  a.Add(a, ConstFloat64(1.0))

adds a constant value to a where a type cast is used to define the constant 1.0.

To differentiate a function f(x,y) = x y^3 + 4, we define

  f := func(x, y ConstScalar) MagicScalar {
    // compute f(x,y) = x*y^3 + 4
    z := NewReal64()
    z.Pow(y, ConstFloat64(3.0))
    z.Mul(z, x)
    z.Add(z, ConstFloat64(4.0))
    return z
  }

that accepts as arguments two ConstScalar variables and returns a MagicScalar. We first define two MagicReal variables

  x := NewReal64(2)
  y := NewReal64(4)

that store the value at which the derivatives should be evaluated. Afterwards, x and y must be activated with

  Variables(2, x, y)

where the first argument sets the order of the derivatives. Variables() should only be called once, as it allocates memory for the given magic variables. In this case, derivatives up to second order are computed. After evaluating f, i.e.

  z := f(x, y)

the function value at (x,y) = (2, 4) can be retrieved with z.GetFloat64(). The first and second partial derivatives can be accessed with z.GetDerivative(i) and z.GetHessian(i, j), where the arguments specify the index of the variable. For instance, the derivative of f with respect to x is returned by z.GetDerivative(0), whereas the derivative with respect to y by z.GetDerivative(1).

Basic linear algebra

Vectors and matrices can be created with

  v := NewDenseFloat64Vector([]float64{1,2})
  m := NewDenseFloat64Matrix([]float64{1,2,3,4}, 2, 2)

  v_ := NewDenseReal64Vector([]float64{1,2})
  m_ := NewDenseReal64Matrix([]float64{1,2,3,4}, 2, 2)

where v has length 2 and m is a 2x2 matrix. With

  v := NullDenseFloat64Vector(2)
  m := NullDenseFloat64Matrix(2, 2)

all values are initially set to zero. Vector and matrix elements can be accessed with the At, MagicAt or ConstAt methods, which return a reference to the scalar implementing either a Scalar, MagicScalar or ConstScalar, i.e.

  m.At(1,1).Add(v.ConstAt(0), v.ConstAt(1))

adds the first two values in v and stores the result in the lower right element of the matrix m. Autodiff supports basic linear algebra operations, for instance, the vector matrix product can be computed with

  w := NullDenseFloat64Vector(2)
  w.MdotV(m, v)

where the result is stored in w. Other operations, such as computing the eigenvalues and eigenvectors of a matrix, require importing the respective package from the algorithm library, i.e.

  import "github.com/pbenner/autodiff/algorithm/eigensystem"

  lambda, _, _ := eigensystem.Run(m)

Examples

Gradient descent

Compare vanilla gradient descent with resilient backpropagation

  import . "github.com/pbenner/autodiff"
  import   "github.com/pbenner/autodiff/algorithm/gradientDescent"
  import   "github.com/pbenner/autodiff/algorithm/rprop"

  f := func(x_ ConstVector) MagicScalar {
    x := x_.ConstAt(0)
    // x^4 - 3x^3 + 2
    r := NewReal64()
    s := NewReal64()
    r.Pow(x.ConstAt(0), ConstFloat64(4.0)
    s.Mul(ConstFloat64(3.0), s.Pow(x, ConstFloat64(3.0)))
    r.Add(ConstFloat64(2.0), r.Add(r, s))
    return r
  }
  x0 := NewDenseFloat64Vector([]float64{8})
  // vanilla gradient descent
  xn1, _ := gradientDescent.Run(f, x0, 0.0001, gradientDescent.Epsilon{1e-8})
  // resilient backpropagation
  xn2, _ := rprop.Run(f, x0, 0.0001, 0.4, rprop.Epsilon{1e-8})

Matrix inversion

Compute the inverse r of a matrix m by minimizing the Frobenius norm ||mb - I||

  import . "github.com/pbenner/autodiff"
  import   "github.com/pbenner/autodiff/algorithm/rprop"

  // define matrix r
  m := NewDenseFloat64Matrix([]float64{1,2,3,4}, 2, 2)
  // create identity matrix I
  I := NullDenseFloat64Matrix(2, 2)
  I.SetIdentity()

  // magic variables for computing the Frobenius norm and its derivative
  t := NewDenseReal64Matrix(2, 2)
  s := NewReal64()
  // objective function
  f := func(x ConstVector) MagicScalar {
    t.Set(x)
    s.Mnorm(t.MsubM(t.MmulM(m, t), I))
    return s
  }
  r, _ := rprop.Run(f, r.GetValues(), 0.01, 0.1, rprop.Epsilon{1e-12})

Newton's method

Find the root of a function f with initial value x0 = (1,1)

  import . "github.com/pbenner/autodiff"
  import   "github.com/pbenner/autodiff/algorithm/newton"

  t := NullReal64()

  f := func(x ConstVector) MagicVector {
    x1 := x.ConstAt(0)
    x2 := x.ConstAt(1)
    y  := NullDenseReal64Vector(2)
    y1 := y.At(0)
    y2 := y.At(1)
    // y1 = x1^2 + x2^2 - 6
    t .Pow(x1, ConstFloat64(2.0))
    y1.Add(y1, t)
    t .Pow(x2, ConstFloat64(2.0))
    y1.Add(y1, t)
    y1.Sub(y1, ConstFloat64(6.0))
    // y2 = x1^3 - x2^2
    t .Pow(x1, ConstFloat64(3.0))
    y2.Add(y2, t)
    t .Pow(x2, ConstFloat64(2.0))
    y2.Sub(y2, t)

    return y
  }

  x0    := NewDenseFloat64Vector([]float64{1,1})
  xn, _ := newton.RunRoot(f, x0, newton.Epsilon{1e-8})

Minimize Rosenbrock's function

Compare Adam, Newton's method, BFGS and Rprop for minimizing Rosenbrock's function

  import   "fmt"

  import . "github.com/pbenner/autodiff"
  import   "github.com/pbenner/autodiff/algorithm/adam"
  import   "github.com/pbenner/autodiff/algorithm/rprop"
  import   "github.com/pbenner/autodiff/algorithm/bfgs"
  import   "github.com/pbenner/autodiff/algorithm/newton"

  f := func(x ConstVector) (MagicScalar, error) {
    // f(x1, x2) = (a - x1)^2 + b(x2 - x1^2)^2
    // a = 1
    // b = 100
    // minimum: (x1,x2) = (a, a^2)
    a := ConstFloat64(  1.0)
    b := ConstFloat64(100.0)
    c := ConstFloat64(  2.0)
    s := NullReal64()
    t := NullReal64()
    s.Pow(s.Sub(a, x.ConstAt(0)), c)
    t.Mul(b, t.Pow(t.Sub(x.ConstAt(1), t.Mul(x.ConstAt(0), x.ConstAt(0))), c))
    s.Add(s, t)
    return s, nil
  }
  hook_adam := func(x, gradient ConstVector, hessian ConstMatrix, y ConstScalar) bool {
    fmt.Println("x       :", x)
    fmt.Println("gradient:", gradient)
    fmt.Println("y       :", y)
    fmt.Println()
    return false
  }
  hook_rprop := func(gradient, step []float64, x ConstVector, y ConstScalar) bool {
    fmt.Println("x       :", x)
    fmt.Println("gradient:", gradient)
    fmt.Println("y       :", y)
    fmt.Println()
    return false
  }
  hook_bfgs := func(x, gradient ConstVector, y ConstScalar) bool {
    fmt.Println("x       :", x)
    fmt.Println("gradient:", gradient)
    fmt.Println("y       :", y)
    fmt.Println()
    return false
  }
  hook_newton := func(x, gradient ConstVector, hessian ConstMatrix, y ConstScalar) bool {
    fmt.Println("x       :", x)
    fmt.Println("gradient:", gradient)
    fmt.Println("y       :", y)
    fmt.Println()
    return false
  }

  x0 := NewDenseFloat64Vector([]float64{-0.5, 2})

  adam.Run(f, x0,
    adam.StepSize{0.1},
    adam.Hook{hook_adam},
    adam.Epsilon{1e-10})

  rprop.Run(f, x0, 0.05, []float64{1.2, 0.8},
    rprop.Hook{hook_rprop},
    rprop.Epsilon{1e-10})

  bfgs.Run(f, x0,
    bfgs.Hook{hook_bfgs},
    bfgs.Epsilon{1e-10})

  newton.RunMin(f, x0,
    newton.HookMin{hook_newton},
    newton.Epsilon{1e-8},
    newton.HessianModification{"LDL"})

Constrained optimization

Maximize the function f(x, y) = x + y subject to x^2 + y^2 = 1 by finding the critical point of the corresponding Lagrangian

  import . "github.com/pbenner/autodiff"
  import   "github.com/pbenner/autodiff/algorithm/newton"

  z := NullReal64()
  t := NullReal64()
  // define the Lagrangian
  f := func(x_ ConstVector) (MagicScalar, error) {
    // z = x + y + lambda(x^2 + y^2 - 1)
    x      := x_.ConstAt(0)
    y      := x_.ConstAt(1)
    lambda := x_.ConstAt(2)
    z.Reset()
    t.Pow(x, ConstFloat64(2.0))
    z.Add(z, t)
    t.Pow(y, ConstFloat64(2.0))
    z.Add(z, t)
    z.Sub(z, ConstFloat64(1.0))
    z.Mul(z, lambda)
    z.Add(z, y)
    z.Add(z, x)

    return z, nil
  }
  // initial value
  x0    := NewDenseFloat64Vector([]float64{3,  5, 1})
  // run Newton's method
  xn, _ := newton.RunCrit(
      f, x0,
      newton.Epsilon{1e-8})