org.apache.commons.math4.fitting.leastsquares.LevenbergMarquardtOptimizer

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

org.apache.commons.math4.fitting.leastsquares.LevenbergMarquardtOptimizer

Marshall Wice
Hello,

I am using the LevenbergMarquardtOptimizer and am receiving very strange results. Running the below code, it appears as though the LevenbergMarquardtOptimizer is diverging rather than converging. In Matlab I run the equivalent code, and the algorithm correctly converges. Is there some bounds that can be set on the algorithm so it does not diverge?

Java Code:
RealVector theta_pr = createZerosVector(9);
double[] gravity = new double[selectedAccData.getColumnDimension()];
Arrays.fill(gravity, Math.pow(9.81744, 2));

LeastSquaresProblem problem = new LeastSquaresBuilder()
        .start(new double[]{0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01})
        .model(costFunc)
        .target(gravity)
        .lazyEvaluation(false)
        .maxEvaluations(150000)
        .maxIterations(6000)
        .build();

LeastSquaresOptimizer.Optimum optimum = new LevenbergMarquardtOptimizer()
        .withCostRelativeTolerance(1.0e-10)
        .optimize(problem);
private MultivariateJacobianFunction costFunc = new MultivariateJacobianFunction() {
    @Override
    public Pair<RealVector, RealMatrix> value(RealVector point) {

        double E1 = point.getEntry(0);
        double E2 = point.getEntry(1);
        double E3 = point.getEntry(2);
        double E4 = point.getEntry(3);
        double E5 = point.getEntry(4);
        double E6 = point.getEntry(5);
        double E7 = point.getEntry(6);
        double E8 = point.getEntry(7);
        double E9 = point.getEntry(8);
       
        RealVector value = new ArrayRealVector(selectedAccData.getColumnDimension());
        RealMatrix jacobian = new Array2DRowRealMatrix(selectedAccData.getColumnDimension(), 9);

        RealMatrix misalignmentMatrix = MatrixUtils.createRealMatrix(new double[][]{{1, -E1, E2}, {0, 1, -E3}, {0 ,0, 1}});
        RealMatrix scalingMatrix = MatrixUtils.createRealDiagonalMatrix(new double[]{E4, E5, E6});
        RealMatrix dMatrix = MatrixUtils.createRealDiagonalMatrix(new double[]{E7, E8, E9});
        RealMatrix oMatrix = createOnesMatrix(3,selectedAccData.getColumnDimension());
        RealMatrix a_bar = misalignmentMatrix.multiply(scalingMatrix).multiply(selectedAccData).subtract(dMatrix.multiply(oMatrix));

        for (int i = 0; i < selectedAccData.getColumnDimension(); ++i) {

            double modelI = (Math.pow(a_bar.getEntry(0,i),2)+Math.pow(a_bar.getEntry(1,i),2)+Math.pow(a_bar.getEntry(2,i),2));
            value.setEntry(i, modelI);

            double a1 = selectedAccData.getEntry(0,i);
            double a2 = selectedAccData.getEntry(1,i);
            double a3 = selectedAccData.getEntry(2,i);

            jacobian.setEntry(i, 0, 2*E5*a2*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3));
            jacobian.setEntry(i, 1, -2*E6*a3*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3));
            jacobian.setEntry(i, 2, 2*E6*a3*(E8 - E5*a2 + E3*E6*a3));
            jacobian.setEntry(i, 3, -2*a1*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3));
            jacobian.setEntry(i, 4, 2*E1*a2*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3) - 2*a2*(E8 - E5*a2 + E3*E6*a3));
            jacobian.setEntry(i, 5, 2*E3*a3*(E8 - E5*a2 + E3*E6*a3) - 2*E2*a3*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3) - 2*a3*(E9 - E6*a3));
            jacobian.setEntry(i, 6, 2*E7 - 2*E4*a1 + 2*E1*E5*a2 - 2*E2*E6*a3);
            jacobian.setEntry(i, 7, 2*E8 - 2*E5*a2 + 2*E3*E6*a3);
            jacobian.setEntry(i, 8, 2*E9 - 2*E6*a3);
        }

        return new Pair<>(value, jacobian);
    }
};

Equivalent Matlab Code:
theta_pr = [0, 0, 0, 0, 0, 0, 0, 0, 0];
ObjectiveFunction = @(theta_pr) costFunc(theta_pr, selectedAccData);
options = optimset('MaxFunEvals', 150000, 'MaxIter', 6000, 'TolFun', 10^(-10));
[theta_pr rsnorm] = lsqnonlin(ObjectiveFunction, theta_pr, [], [], options);
function [ res_vector ] = costFunc( E, a_hat )

misalignmentMatrix = [1, -E(1), E(2); 0, 1, -E(3); 0, 0, 1];
scalingMtrix = diag([E(4), E(5), E(6)]);

a_bar = misalignmentMatrix*scalingMtrix*(a_hat) - (diag([E(7), E(8), E(9)])*ones(3, length(a_hat)));

residuals = zeros(length(a_bar(1,:)), 1);

for i = 1:length(a_bar(1,:))
residuals(i,1) = 9.81744^2 - (a_bar(1,i)^2 + a_bar(2,i)^2 + a_bar(3,i)^2)
end

        res_vector = residuals;

end

If this is not the correct forum to be asking this, or if there is a better way to communicate and share my code, please let me know.

Thank you,

Marshall
Reply | Threaded
Open this post in threaded view
|

Re: org.apache.commons.math4.fitting.leastsquares.LevenbergMarquardtOptimizer

Gilles Sadowski-2
Hi.

Le lun. 10 juin 2019 à 11:01, Marshall Wice <[hidden email]> a écrit :
>
> Hello,
>
> I am using the LevenbergMarquardtOptimizer and am receiving very strange results. Running the below code, it appears as though the LevenbergMarquardtOptimizer is diverging rather than converging.

Sorry that you encounter problems with that code. :-(

> In Matlab I run the equivalent code, and the algorithm correctly converges. Is there some bounds that can be set on the algorithm so it does not diverge?

Bounds are not supported in this implementation.
Users can set them indirectly through a "logit" transform on the model
parameters.

>
> Java Code:
> RealVector theta_pr = createZerosVector(9);
> double[] gravity = new double[selectedAccData.getColumnDimension()];
> Arrays.fill(gravity, Math.pow(9.81744, 2));
>
> LeastSquaresProblem problem = new LeastSquaresBuilder()
>         .start(new double[]{0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01})
>         .model(costFunc)
>         .target(gravity)
>         .lazyEvaluation(false)
>         .maxEvaluations(150000)
>         .maxIterations(6000)
>         .build();
>
> LeastSquaresOptimizer.Optimum optimum = new LevenbergMarquardtOptimizer()
>         .withCostRelativeTolerance(1.0e-10)
>         .optimize(problem);
> private MultivariateJacobianFunction costFunc = new MultivariateJacobianFunction() {
>     @Override
>     public Pair<RealVector, RealMatrix> value(RealVector point) {
>
>         double E1 = point.getEntry(0);
>         double E2 = point.getEntry(1);
>         double E3 = point.getEntry(2);
>         double E4 = point.getEntry(3);
>         double E5 = point.getEntry(4);
>         double E6 = point.getEntry(5);
>         double E7 = point.getEntry(6);
>         double E8 = point.getEntry(7);
>         double E9 = point.getEntry(8);
>
>         RealVector value = new ArrayRealVector(selectedAccData.getColumnDimension());
>         RealMatrix jacobian = new Array2DRowRealMatrix(selectedAccData.getColumnDimension(), 9);
>
>         RealMatrix misalignmentMatrix = MatrixUtils.createRealMatrix(new double[][]{{1, -E1, E2}, {0, 1, -E3}, {0 ,0, 1}});
>         RealMatrix scalingMatrix = MatrixUtils.createRealDiagonalMatrix(new double[]{E4, E5, E6});
>         RealMatrix dMatrix = MatrixUtils.createRealDiagonalMatrix(new double[]{E7, E8, E9});
>         RealMatrix oMatrix = createOnesMatrix(3,selectedAccData.getColumnDimension());
>         RealMatrix a_bar = misalignmentMatrix.multiply(scalingMatrix).multiply(selectedAccData).subtract(dMatrix.multiply(oMatrix));
>
>         for (int i = 0; i < selectedAccData.getColumnDimension(); ++i) {
>
>             double modelI = (Math.pow(a_bar.getEntry(0,i),2)+Math.pow(a_bar.getEntry(1,i),2)+Math.pow(a_bar.getEntry(2,i),2));
>             value.setEntry(i, modelI);
>
>             double a1 = selectedAccData.getEntry(0,i);
>             double a2 = selectedAccData.getEntry(1,i);
>             double a3 = selectedAccData.getEntry(2,i);
>
>             jacobian.setEntry(i, 0, 2*E5*a2*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3));
>             jacobian.setEntry(i, 1, -2*E6*a3*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3));
>             jacobian.setEntry(i, 2, 2*E6*a3*(E8 - E5*a2 + E3*E6*a3));
>             jacobian.setEntry(i, 3, -2*a1*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3));
>             jacobian.setEntry(i, 4, 2*E1*a2*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3) - 2*a2*(E8 - E5*a2 + E3*E6*a3));
>             jacobian.setEntry(i, 5, 2*E3*a3*(E8 - E5*a2 + E3*E6*a3) - 2*E2*a3*(E7 - E4*a1 + E1*E5*a2 - E2*E6*a3) - 2*a3*(E9 - E6*a3));
>             jacobian.setEntry(i, 6, 2*E7 - 2*E4*a1 + 2*E1*E5*a2 - 2*E2*E6*a3);
>             jacobian.setEntry(i, 7, 2*E8 - 2*E5*a2 + 2*E3*E6*a3);
>             jacobian.setEntry(i, 8, 2*E9 - 2*E6*a3);
>         }
>
>         return new Pair<>(value, jacobian);
>     }
> };

The above calls methods that are not defined in this code excerpt.

Bette is to file a report on the issue-tracking system:
   https://issues.apache.org/jira/projects/MATH

There you should attached a minimal working example that shows
the problem, preferably in the form of a unit test, like those that
are already provided:
    https://gitbox.apache.org/repos/asf?p=commons-math.git;a=tree;f=src/test/java/org/apache/commons/math4/fitting/leastsquares

> Equivalent Matlab Code:
> theta_pr = [0, 0, 0, 0, 0, 0, 0, 0, 0];
> ObjectiveFunction = @(theta_pr) costFunc(theta_pr, selectedAccData);
> options = optimset('MaxFunEvals', 150000, 'MaxIter', 6000, 'TolFun', 10^(-10));
> [theta_pr rsnorm] = lsqnonlin(ObjectiveFunction, theta_pr, [], [], options);
> function [ res_vector ] = costFunc( E, a_hat )
>
> misalignmentMatrix = [1, -E(1), E(2); 0, 1, -E(3); 0, 0, 1];
> scalingMtrix = diag([E(4), E(5), E(6)]);
>
> a_bar = misalignmentMatrix*scalingMtrix*(a_hat) - (diag([E(7), E(8), E(9)])*ones(3, length(a_hat)));
>
> residuals = zeros(length(a_bar(1,:)), 1);
>
> for i = 1:length(a_bar(1,:))
> residuals(i,1) = 9.81744^2 - (a_bar(1,i)^2 + a_bar(2,i)^2 + a_bar(3,i)^2)
> end
>
>         res_vector = residuals;
>
> end
>
> If this is not the correct forum to be asking this, or if there is a better way to communicate and share my code, please let me know.

Cf. above.

Thanks for the report,
Gilles

>
> Thank you,
>
> Marshall

---------------------------------------------------------------------
To unsubscribe, e-mail: [hidden email]
For additional commands, e-mail: [hidden email]