Automatic Symbolic Differentiation

Chain rule based automatic differentiation using a symbolic interface

https://github.com/lessthanoptimal/ddogleg/tree/v1.0.0/examples/src/org/ddogleg/example/ExampleChainRuleMinimizationSym.java

 1public static void main( String[] args ) {
 2    var residuals = new Residuals();
 3
 4    UnconstrainedMinimization optimizer = FactoryOptimization.quasiNewtonBfgs(null);
 5    optimizer.setVerbose(System.out, 0);
 6
 7    // f(x) = Σ r_i², and its gradient ∇f = 2·Jᵀ·r, from the residuals and their Jacobian
 8    optimizer.setFunction(new LsToNonLinear(residuals), new LsToNonLinearDeriv<>(residuals, residuals), 0);
 9
10    // [-1.2,  1] is the recommended starting point for testing
11    optimizer.initialize(new double[]{-1.2, 1}, 1e-12, 1e-12);
12    UtilOptimize.process(optimizer, 500);
13
14    double[] found = optimizer.getParameters();
15    System.out.println("Final Error = " + optimizer.getFunctionValue());
16    System.out.printf("x[0]: expected=1.00  found=%5.2f\n", found[0]);
17    System.out.printf("x[1]: expected=1.00  found=%5.2f\n", found[1]);
18}
19
20/// The two Rosenbrock residuals and their Jacobian
21public static class Residuals implements FunctionNtoM, FunctionNtoMxN<DMatrixRMaj> {
22    final Sym symX = Sym.param("x", 2);
23    final SymFunction function;
24    final CRArray<double[]> output = new CRArray<>(new double[2], 0);
25
26    DMatrixRMaj jacobian = new DMatrixRMaj(2, 2);
27
28    public Residuals() {
29        // x1 = x[0],  x2 = x[1]
30        Sym x1 = SymOps.remap(1, symX, 0, -1);
31        Sym x2 = SymOps.remap(1, symX, -1, 0);
32
33        // r0 = 10·(x2 − x1²)
34        Sym r0 = SymOps.residual(10.0, SymOps.mult(x1, x1), x2);
35
36        // r1 = 1 − x1
37        Sym r1 = SymOps.residual(1.0, x1, Sym.constant(1.0));
38
39        // r = [r0, r1]
40        function = new SymGraphBuilder().build(SymOps.concat(r0, r1));
41        function.feed(symX, new double[2]);
42        function.mapParamColumns();
43
44        // every derivative is added into the Jacobian, so it's zeroed before each evaluation
45        function.jacobianTo(( row, col, values, numRows, numCols ) -> {
46            for (int i = 0; i < numRows; i++) {
47                for (int j = 0; j < numCols; j++) {
48                    jacobian.data[(row + i)*2 + col + j] += values[i*numCols + j];
49                }
50            }
51        });
52    }
53
54    /// r(x)
55    @Override public void process( double[] input, double[] output ) {
56        function.feed(symX, input);
57        this.output.wrap(output, 0);
58        function.run(this.output);
59    }
60
61    /// J(x) = ∂r/∂x
62    @Override public void process( double[] input, DMatrixRMaj output ) {
63        jacobian = output;
64        jacobian.zero();
65        function.feed(symX, input);
66        function.jacobianRow = 0;
67        function.runJacobian();
68    }
69
70    @Override public DMatrixRMaj declareMatrixMxN() {return new DMatrixRMaj(2, 2);}
71
72    @Override public int getNumOfInputsN() {return 2;}
73
74    @Override public int getNumOfOutputsM() {return 2;}
75}