LevenbergMarquardtOptimizer implements LeastSquaresOptimizer { public Optimum optimize(final LeastSquaresProblem problem) { final int nR = problem.getObservationSize(); final int nC = problem.getParameterSize(); final Incrementor iterationCounter = problem.getIterationCounter(); final Incrementor evaluationCounter = problem.getEvaluationCounter(); final ConvergenceChecker<Evaluation> checker = problem.getConvergenceChecker(); final int solvedCols = FastMath.min(nR, nC); double[] lmDir = new double[nC]; double lmPar = 0; double delta = 0; double xNorm = 0; double[] diag = new double[nC]; double[] oldX = new double[nC]; double[] oldRes = new double[nR]; double[] qtf = new double[nR]; double[] work1 = new double[nC]; double[] work2 = new double[nC]; double[] work3 = new double[nC]; evaluationCounter.incrementCount(); Evaluation current = problem.evaluate(problem.getStart()); double[] currentResiduals = current.getResiduals().toArray(); double currentCost = current.getCost(); double[] currentPoint = current.getPoint().toArray(); boolean firstIteration = true; while (true) { iterationCounter.incrementCount(); final Evaluation previous = current; final InternalData internalData = qrDecomposition(current.getJacobian(), solvedCols); final double[][] weightedJacobian = internalData.weightedJacobian; final int[] permutation = internalData.permutation; final double[] diagR = internalData.diagR; final double[] jacNorm = internalData.jacNorm; double[] weightedResidual = currentResiduals; for (int i = 0; i < nR; i++) { qtf[i] = weightedResidual[i]; } qTy(qtf, internalData); for (int k = 0; k < solvedCols; ++k) { int pk = permutation[k]; weightedJacobian[k][pk] = diagR[pk]; } if (firstIteration) { xNorm = 0; for (int k = 0; k < nC; ++k) { double dk = jacNorm[k]; if (dk == 0) { dk = 1.0; } double xk = dk * currentPoint[k]; xNorm += xk * xk; diag[k] = dk; } xNorm = FastMath.sqrt(xNorm); delta = (xNorm == 0) ? initialStepBoundFactor : (initialStepBoundFactor * xNorm); } double maxCosine = 0; if (currentCost != 0) { for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double s = jacNorm[pj]; if (s != 0) { double sum = 0; for (int i = 0; i <= j; ++i) { sum += weightedJacobian[i][pj] * qtf[i]; } maxCosine = FastMath.max(maxCosine, FastMath.abs(sum) / (s * currentCost)); } } } if (maxCosine <= orthoTolerance) { return new OptimumImpl( current, evaluationCounter.getCount(), iterationCounter.getCount()); } for (int j = 0; j < nC; ++j) { diag[j] = FastMath.max(diag[j], jacNorm[j]); } for (double ratio = 0; ratio < 1.0e-4;) { for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; oldX[pj] = currentPoint[pj]; } final double previousCost = currentCost; double[] tmpVec = weightedResidual; weightedResidual = oldRes; oldRes = tmpVec; lmPar = determineLMParameter(qtf, delta, diag, internalData, solvedCols, work1, work2, work3, lmDir, lmPar); double lmNorm = 0; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; lmDir[pj] = -lmDir[pj]; currentPoint[pj] = oldX[pj] + lmDir[pj]; double s = diag[pj] * lmDir[pj]; lmNorm += s * s; } lmNorm = FastMath.sqrt(lmNorm); if (firstIteration) { delta = FastMath.min(delta, lmNorm); } evaluationCounter.incrementCount(); current = problem.evaluate(new ArrayRealVector(currentPoint)); currentResiduals = current.getResiduals().toArray(); currentCost = current.getCost(); currentPoint = current.getPoint().toArray(); double actRed = -1.0; if (0.1 * currentCost < previousCost) { double r = currentCost / previousCost; actRed = 1.0 - r * r; } for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; double dirJ = lmDir[pj]; work1[j] = 0; for (int i = 0; i <= j; ++i) { work1[i] += weightedJacobian[i][pj] * dirJ; } } double coeff1 = 0; for (int j = 0; j < solvedCols; ++j) { coeff1 += work1[j] * work1[j]; } double pc2 = previousCost * previousCost; coeff1 /= pc2; double coeff2 = lmPar * lmNorm * lmNorm / pc2; double preRed = coeff1 + 2 * coeff2; double dirDer = -(coeff1 + coeff2); ratio = (preRed == 0) ? 0 : (actRed / preRed); if (ratio <= 0.25) { double tmp = (actRed < 0) ? (0.5 * dirDer / (dirDer + 0.5 * actRed)) : 0.5; if ((0.1 * currentCost >= previousCost) || (tmp < 0.1)) { tmp = 0.1; } delta = tmp * FastMath.min(delta, 10.0 * lmNorm); lmPar /= tmp; } else if ((lmPar == 0) || (ratio >= 0.75)) { delta = 2 * lmNorm; lmPar *= 0.5; } if (ratio >= 1.0e-4) { firstIteration = false; xNorm = 0; for (int k = 0; k < nC; ++k) { double xK = diag[k] * currentPoint[k]; xNorm += xK * xK; } xNorm = FastMath.sqrt(xNorm); if (checker != null && checker.converged(iterationCounter.getCount(), previous, current)) { return new OptimumImpl(current, evaluationCounter.getCount(), iterationCounter.getCount()); } } else { currentCost = previousCost; for (int j = 0; j < solvedCols; ++j) { int pj = permutation[j]; currentPoint[pj] = oldX[pj]; } tmpVec = weightedResidual; weightedResidual = oldRes; oldRes = tmpVec; current = previous; } if ((FastMath.abs(actRed) <= costRelativeTolerance && preRed <= costRelativeTolerance && ratio <= 2.0) || delta <= parRelativeTolerance * xNorm) { return new OptimumImpl(current, evaluationCounter.getCount(), iterationCounter.getCount()); } if (FastMath.abs(actRed) <= TWO_EPS && preRed <= TWO_EPS && ratio <= 2.0) { throw new ConvergenceException(LocalizedFormats.TOO_SMALL_COST_RELATIVE_TOLERANCE, costRelativeTolerance); } else if (delta <= TWO_EPS * xNorm) { throw new ConvergenceException(LocalizedFormats.TOO_SMALL_PARAMETERS_RELATIVE_TOLERANCE, parRelativeTolerance); } else if (maxCosine <= TWO_EPS) { throw new ConvergenceException(LocalizedFormats.TOO_SMALL_ORTHOGONALITY_TOLERANCE, orthoTolerance); } } } } LevenbergMarquardtOptimizer();  LevenbergMarquardtOptimizer(
            final double initialStepBoundFactor,
            final double costRelativeTolerance,
            final double parRelativeTolerance,
            final double orthoTolerance,
            final double qrRankingThreshold); LevenbergMarquardtOptimizer withInitialStepBoundFactor(double newInitialStepBoundFactor); LevenbergMarquardtOptimizer withCostRelativeTolerance(double newCostRelativeTolerance); LevenbergMarquardtOptimizer withParameterRelativeTolerance(double newParRelativeTolerance); LevenbergMarquardtOptimizer withOrthoTolerance(double newOrthoTolerance); LevenbergMarquardtOptimizer withRankingThreshold(double newQRRankingThreshold); double getInitialStepBoundFactor(); double getCostRelativeTolerance(); double getParameterRelativeTolerance(); double getOrthoTolerance(); double getRankingThreshold(); Optimum optimize(final LeastSquaresProblem problem);  }
@Override @Test public void testNonInvertible() { try{ LinearProblem problem = new LinearProblem(new double[][] { { 1, 2, -3 }, { 2, 1, 3 }, { -3, 0, -9 } }, new double[] { 1, 1, 1 }); final Optimum optimum = optimizer.optimize( problem.getBuilder().maxIterations(20).build()); Assert.assertTrue(FastMath.sqrt(problem.getTarget().length) * optimum.getRMS() > 0.6); optimum.getCovariances(1.5e-14); fail(optimizer); }catch (SingularMatrixException e){ } }
@Test public void testBevington() { final double[][] dataPoints = { { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, 315, 330, 345, 360, 375, 390, 405, 420, 435, 450, 465, 480, 495, 510, 525, 540, 555, 570, 585, 600, 615, 630, 645, 660, 675, 690, 705, 720, 735, 750, 765, 780, 795, 810, 825, 840, 855, 870, 885, }, { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89, 74, 61, 66, 68, 48, 54, 51, 46, 55, 29, 28, 37, 49, 26, 35, 29, 31, 24, 25, 35, 24, 30, 26, 28, 21, 18, 20, 27, 17, 17, 14, 17, 24, 11, 22, 17, 12, 10, 13, 16, 9, 9, 14, 21, 17, 13, 12, 18, 10, }, }; final double[] start = {10, 900, 80, 27, 225}; final BevingtonProblem problem = new BevingtonProblem(); final int len = dataPoints[0].length; final double[] weights = new double[len]; for (int i = 0; i < len; i++) { problem.addPoint(dataPoints[0][i], dataPoints[1][i]); weights[i] = 1 / dataPoints[1][i]; } final Optimum optimum = optimizer.optimize( builder(problem) .target(dataPoints[1]) .weight(new DiagonalMatrix(weights)) .start(start) .maxIterations(20) .build() ); final RealVector solution = optimum.getPoint(); final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 }; final RealMatrix covarMatrix = optimum.getCovariances(1e-14); final double[][] expectedCovarMatrix = { { 3.38, -3.69, 27.98, -2.34, -49.24 }, { -3.69, 2492.26, 81.89, -69.21, -8.9 }, { 27.98, 81.89, 468.99, -44.22, -615.44 }, { -2.34, -69.21, -44.22, 6.39, 53.80 }, { -49.24, -8.9, -615.44, 53.8, 929.45 } }; final int numParams = expectedSolution.length; for (int i = 0; i < numParams; i++) { final double error = FastMath.sqrt(expectedCovarMatrix[i][i]); Assert.assertEquals("Parameter " + i, expectedSolution[i], solution.getEntry(i), error); } for (int i = 0; i < numParams; i++) { for (int j = 0; j < numParams; j++) { Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]", expectedCovarMatrix[i][j], covarMatrix.getEntry(i, j), FastMath.abs(0.1 * expectedCovarMatrix[i][j])); } } }
@Test public void testCircleFitting2() { final double xCenter = 123.456; final double yCenter = 654.321; final double xSigma = 10; final double ySigma = 15; final double radius = 111.111; final long seed = 59421061L; final RandomCirclePointGenerator factory = new RandomCirclePointGenerator(xCenter, yCenter, radius, xSigma, ySigma, seed); final CircleProblem circle = new CircleProblem(xSigma, ySigma); final int numPoints = 10; for (Vector2D p : factory.generate(numPoints)) { circle.addPoint(p.getX(), p.getY()); } final double[] init = { 90, 659, 115 }; final Optimum optimum = optimizer.optimize( builder(circle).maxIterations(50).start(init).build()); final double[] paramFound = optimum.getPoint().toArray(); final double[] asymptoticStandardErrorFound = optimum.getSigma(1e-14).toArray(); Assert.assertEquals(xCenter, paramFound[0], asymptoticStandardErrorFound[0]); Assert.assertEquals(yCenter, paramFound[1], asymptoticStandardErrorFound[1]); Assert.assertEquals(radius, paramFound[2], asymptoticStandardErrorFound[2]); }
@Test public void testParameterValidator() { final double xCenter = 123.456; final double yCenter = 654.321; final double xSigma = 10; final double ySigma = 15; final double radius = 111.111; final long seed = 3456789L; final RandomCirclePointGenerator factory = new RandomCirclePointGenerator(xCenter, yCenter, radius, xSigma, ySigma, seed); final CircleProblem circle = new CircleProblem(xSigma, ySigma); final int numPoints = 10; for (Vector2D p : factory.generate(numPoints)) { circle.addPoint(p.getX(), p.getY()); } final double[] init = { 90, 659, 115 }; final Optimum optimum = optimizer.optimize(builder(circle).maxIterations(50).start(init).build()); final int numEval = optimum.getEvaluations(); Assert.assertTrue(numEval > 1); final ParameterValidator cheatValidator = new ParameterValidator() { public RealVector validate(RealVector params) { return optimum.getPoint(); } }; final Optimum cheatOptimum = optimizer.optimize(builder(circle).maxIterations(50).start(init).parameterValidator(cheatValidator).build()); final int cheatNumEval = cheatOptimum.getEvaluations(); Assert.assertTrue(cheatNumEval < numEval); }
@Test public void testEvaluationCount() { LeastSquaresProblem lsp = new LinearProblem(new double[][] {{1}}, new double[] {1}) .getBuilder() .checker(new ConvergenceChecker<Evaluation>() { public boolean converged(int iteration, Evaluation previous, Evaluation current) { return true; } }) .build(); Optimum optimum = optimizer.optimize(lsp); Assert.assertThat(optimum.getIterations(), is(1)); Assert.assertThat(optimum.getEvaluations(), is(2)); }