Non-Linear Least Squares¶
The following code demonstrates non-linear least squares minimization by finding the best fit line to a set of points. See the optimization manual for an overview of all optimization techniques available in DDogleg. Unconstrained least-squares minimization solves problems which can be described by a function of the form:
\[\min\limits_{\boldsymbol{x} \in \Re^N} f(\boldsymbol{x})=\frac{1}{2}\sum^m_{i=1} r^2_i(\boldsymbol{x})\]
where \(r_i(\boldsymbol{x}) = f_i(\boldsymbol{x}) - y_i\) is a scalar function which outputs the residual for function \(i\), predicted value subtracted the observed value. By definition \(f(\boldsymbol{x}) \ge 0\). In DDogleg you don’t define the \(f(\boldsymbol{x})\) function directly but instead define the set of \(r_j(\boldsymbol{x})\) functions by implementing FunctionNtoM. Implementations of FunctionNtoM take in an array with N elements and output an array with M elements.
In this example we will consider a very simple unconstrained least-squares problem, fitting a line to a set of points. To solve non-linear problems you will need to select which method to use (Levenberg-Mardquardt is usually a good choice), define the function to minimize, and select an initial value. You can also specify a Jacobian, an N by M matrix, which is the residual functions’ derivative. If you do not specify this function then it will be computed for you numerically.
The specific function being optimized is below:
\[t_i = \frac{-p_1(x_i-p_0) + p_0(y_i-p_1)}{p_0^2 + p_i^2}\]\[\begin{split}f_{2i}(\boldsymbol{p}) &= x_i -p_0 + t_i p_1 \\ f_{2i+1}(\boldsymbol{p}) &= y_i -p_1 - t_i p_0 \\\end{split}\]
where \((x_i,y_i)\) is a point being fit to and \(\boldsymbol{p}=(p_0,p_1)\) is the line being fit to. The line has been parameterized using the closest point on the line to the origin. Thus the line’s slope is \(\boldsymbol{p}=(-p_1,p_0)\).
The code below walks you through each step.
Here it is shown how to set up and perform the optimization. DDogleg is a bit unusual in that it allows you to invoke each step of the optimization manually. The advantage of that is that it is very easy to abort your optimization early if you run out of time. For sake of brevity, we use UtilProcess.optimize() to handle the iterations.
1public static void main( String[] args ) {
2 // define a line in 2D space as the tangent from the origin
3 double lineX = -2.1;
4 double lineY = 1.3;
5
6 // randomly generate points along the line
7 var rand = new Random(234);
8 var points = new ArrayList<Point2D>();
9 for (int i = 0; i < 20; i++) {
10 double t = (rand.nextDouble() - 0.5)*10;
11 points.add(new Point2D(lineX + t*lineY, lineY - t*lineX));
12 }
13
14 // Define the function being optimized and create the optimizer
15 var func = new FunctionLineDistanceEuclidean(points);
16 UnconstrainedLeastSquares<DMatrixRMaj> optimizer = FactoryOptimization.levenbergMarquardt(null, true);
17
18 // Send to standard out progress information
19 optimizer.setVerbose(System.out, 0);
20
21 // if no jacobian is specified it will be computed numerically
22 optimizer.setFunction(func, null);
23
24 // provide it an extremely crude initial estimate of the line equation
25 optimizer.initialize(new double[]{-0.5, 0.5}, 1e-12, 1e-12);
26
27 // iterate 500 times or until it converges.
28 // Manually iteration is possible too if more control over is required
29 UtilOptimize.process(optimizer, 500);
30
31 double[] found = optimizer.getParameters();
32
33 // see how accurately it found the solution
34 System.out.println("Final Error = " + optimizer.getFunctionValue());
35
36 // Compare the actual parameters to the found parameters
37 System.out.printf("Actual lineX %5.2f found %5.2f\n", lineX, found[0]);
38 System.out.printf("Actual lineY %5.2f found %5.2f\n", lineY, found[1]);
39}
When you run this example you should see something like the following:
Steps fx change |step| f-test g-test tr-ratio lambda
0 5.117E+01 0.000E+00 0.000E+00 0.000E+00 0.000E+00 0.00 1.00E-04
1 1.685E+01 -3.431E+01 1.793E+00 6.706E-01 2.279E+02 0.709 9.27E-05
2 6.133E-02 -1.679E+01 4.170E-01 9.964E-01 4.849E+01 1.026 3.09E-05
3 9.294E-06 -6.133E-02 6.479E-02 9.998E-01 3.179E-01 1.000 1.03E-05
4 2.136E-14 -9.294E-06 3.236E-04 1.000E+00 3.013E-02 1.000 3.43E-06
5 3.579E-24 -2.136E-14 3.833E-08 1.000E+00 1.929E-07 1.000 1.14E-06
6 2.792E-30 -3.579E-24 5.888E-13 1.000E+00 5.813E-12 1.000 3.81E-07
Converged g-test
Final Error = 2.791828047383737E-30
Actual lineX -2.10 found -2.10
Actual lineY 1.30 found 1.30
The java code for the function being optimized is shown below. It was mathematically described above previously.
1public class FunctionLineDistanceEuclidean implements FunctionNtoM {
2
3 // Data which the line is being fit too
4 List<Point2D> data;
5
6 public FunctionLineDistanceEuclidean( List<Point2D> data ) {
7 this.data = data;
8 }
9
10 /**
11 * Number of parameters used to define the line.
12 */
13 @Override
14 public int getNumOfInputsN() {
15 return 2;
16 }
17
18 /**
19 * Number of output error functions. Two for each point.
20 */
21 @Override
22 public int getNumOfOutputsM() {
23 return data.size()*2;
24 }
25
26 @Override
27 public void process( double[] input, double[] output ) {
28
29 // tangent equation
30 double tanX = input[0], tanY = input[1];
31
32 // convert into parametric equation
33 double lineX = tanX;
34 double lineY = tanY;
35 double slopeX = -tanY;
36 double slopeY = tanX;
37
38 // compute the residual error for each point in the data set
39 for (int i = 0; i < data.size(); i++) {
40 Point2D p = data.get(i);
41
42 double t = slopeX*(p.x - lineX) + slopeY*(p.y - lineY);
43 t /= slopeX*slopeX + slopeY*slopeY;
44
45 double closestX = lineX + t*slopeX;
46 double closestY = lineY + t*slopeY;
47
48 output[i*2] = p.x - closestX;
49 output[i*2 + 1] = p.y - closestY;
50 }
51 }
52}