LevenbergMarquardtOptimizer unable to perform Q.R decomposition on the 107x2 jacobian matrix

StackOverflow https://stackoverflow.com/questions/19116987

  •  30-06-2022
  •  | 
  •  

Question

I really need some help to get the LevenbergMarquardtOptimizer up and running. Please see the following compileable code (needs apache common math 3). No matter how I try to return the gradients an Exception will raise ... hmmm get stuck ... 

package tryout;

import java.sql.Date;
import java.util.Calendar;
import java.util.Arrays;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math.util.FastMath;
import org.apache.commons.math3.analysis.ParametricUnivariateFunction;
import org.apache.commons.math3.analysis.differentiation.DerivativeStructure;
import org.apache.commons.math3.fitting.CurveFitter;
import org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer;

import static org.apache.commons.math.util.FastMath.pow;
import static org.apache.commons.math.util.FastMath.cos;
import static org.apache.commons.math.util.FastMath.sin;

public class Sornette2NoLogPriceLevenberg {
    static String[] T = new String[]{"2010-04-01","2010-04-05","2010-04-06","2010-04-07","2010-04-08","2010-04-09","2010-04-12","2010-04-13","2010-04-14","2010-04-15","2010-04-16","2010-04-19","2010-04-20","2010-04-21","2010-04-22","2010-04-23","2010-04-26","2010-04-27","2010-04-28","2010-04-29","2010-04-30","2010-05-03","2010-05-04","2010-05-05","2010-05-06","2010-05-07","2010-05-10","2010-05-11","2010-05-12","2010-05-13","2010-05-14","2010-05-18","2010-05-19","2010-05-20","2010-05-21","2010-05-24","2010-05-25","2010-05-26","2010-05-27","2010-05-28","2010-06-01","2010-06-02","2010-06-03","2010-06-04","2010-06-07","2010-06-08","2010-06-09","2010-06-10","2010-06-11","2010-06-14","2010-06-15","2010-06-16","2010-06-17","2010-06-18","2010-06-21","2010-06-22","2010-06-23","2010-06-24","2010-06-25","2010-06-28","2010-06-29","2010-06-30","2010-07-01","2010-07-02","2010-07-06","2010-07-07","2010-07-08","2010-07-09","2010-07-12","2010-07-13","2010-07-14","2010-07-15","2010-07-16","2010-07-19","2010-07-20","2010-07-21","2010-07-22","2010-07-23","2010-07-26","2010-07-27","2010-07-28","2010-07-29","2010-07-30","2010-08-02","2010-08-03","2010-08-04","2010-08-05","2010-08-06","2010-08-09","2010-08-10","2010-08-11","2010-08-12","2010-08-13","2010-08-16","2010-08-17","2010-08-18","2010-08-19","2010-08-20","2010-08-23","2010-08-24","2010-08-25","2010-08-26","2010-08-27","2010-08-30","2010-08-31","2010-09-01","2010-09-02","2010-09-03","2010-09-07","2010-09-08","2010-09-09","2010-09-10","2010-09-13","2010-09-14","2010-09-15","2010-09-16","2010-09-17","2010-09-20","2010-09-21","2010-09-22","2010-09-23","2010-09-24","2010-09-27","2010-09-28","2010-09-29","2010-09-30","2010-10-01","2010-10-04","2010-10-05","2010-10-06","2010-10-07","2010-10-08","2010-10-11","2010-10-12","2010-10-13","2010-10-14","2010-10-15","2010-10-18","2010-10-19","2010-10-20","2010-10-21","2010-10-22","2010-10-25","2010-10-26","2010-10-27","2010-10-28","2010-10-29","2010-11-01","2010-11-02","2010-11-03","2010-11-04","2010-11-05","2010-11-08","2010-11-09","2010-11-10","2010-11-11","2010-11-12","2010-11-15","2010-11-16","2010-11-17","2010-11-18","2010-11-19","2010-11-22","2010-11-23","2010-11-24","2010-11-26","2010-11-29","2010-11-30","2010-12-01","2010-12-02","2010-12-03","2010-12-06","2010-12-07","2010-12-08","2010-12-09","2010-12-10","2010-12-13","2010-12-14","2010-12-15","2010-12-16","2010-12-17","2010-12-20","2010-12-21","2010-12-22","2010-12-23","2010-12-27","2010-12-28","2010-12-29","2010-12-30","2010-12-31","2011-01-03","2011-01-04","2011-01-05","2011-01-06","2011-01-07","2011-01-10","2011-01-11","2011-01-12","2011-01-13","2011-01-14","2011-01-18","2011-01-19","2011-01-20","2011-01-21","2011-01-24","2011-01-25","2011-01-26","2011-01-27","2011-01-28","2011-01-31","2011-02-01","2011-02-02","2011-02-03","2011-02-04","2011-02-07","2011-02-08","2011-02-09","2011-02-10","2011-02-11","2011-02-14","2011-02-15","2011-02-16","2011-02-17","2011-02-18","2011-02-22","2011-02-23","2011-02-24","2011-02-25","2011-02-28","2011-03-01","2011-03-02","2011-03-03","2011-03-04","2011-03-07","2011-03-08","2011-03-09","2011-03-10","2011-03-11","2011-03-14","2011-03-15","2011-03-16","2011-03-17","2011-03-18","2011-03-21","2011-03-22","2011-03-23","2011-03-24","2011-03-25","2011-03-28","2011-03-29","2011-03-30","2011-03-31","2011-04-01","2011-04-04","2011-04-05","2011-04-06","2011-04-07","2011-04-08","2011-04-11","2011-04-12","2011-04-13","2011-04-14","2011-04-15","2011-04-18","2011-04-19","2011-04-20","2011-04-21","2011-04-25","2011-04-26","2011-04-27","2011-04-28","2011-04-29","2011-05-02","2011-05-03","2011-05-04","2011-05-05","2011-05-06","2011-05-09","2011-05-10","2011-05-11","2011-05-12","2011-05-13","2011-05-16","2011-05-17","2011-05-18","2011-05-19","2011-05-20","2011-05-23","2011-05-24","2011-05-25","2011-05-26","2011-05-27","2011-05-31","2011-06-01","2011-06-02","2011-06-03","2011-06-06","2011-06-07","2011-06-08","2011-06-09","2011-06-10","2011-06-13","2011-06-14","2011-06-15","2011-06-16","2011-06-17","2011-06-20","2011-06-21","2011-06-22","2011-06-23","2011-06-24","2011-06-27","2011-06-28","2011-06-29","2011-06-30","2011-07-01","2011-07-05","2011-07-06","2011-07-07","2011-07-08","2011-07-11","2011-07-12","2011-07-13","2011-07-14","2011-07-15","2011-07-18","2011-07-19","2011-07-20","2011-07-21","2011-07-22","2011-07-25","2011-07-26","2011-07-27","2011-07-28","2011-07-29","2011-08-01","2011-08-02","2011-08-03","2011-08-04","2011-08-05","2011-08-08","2011-08-09","2011-08-10","2011-08-11","2011-08-12","2011-08-15","2011-08-16","2011-08-17","2011-08-18","2011-08-19","2011-08-22","2011-08-23","2011-08-24","2011-08-25","2011-08-26","2011-08-29","2011-08-30","2011-08-31","2011-09-01","2011-09-02","2011-09-06","2011-09-07","2011-09-08","2011-09-09","2011-09-12","2011-09-13","2011-09-14","2011-09-15","2011-09-16","2011-09-19","2011-09-20","2011-09-21","2011-09-22","2011-09-23","2011-09-26","2011-09-27","2011-09-28","2011-09-29","2011-09-30","2011-10-03","2011-10-04","2011-10-05","2011-10-06","2011-10-07","2011-10-10","2011-10-11","2011-10-12","2011-10-13","2011-10-14","2011-10-17","2011-10-18","2011-10-19","2011-10-20","2011-10-21","2011-10-24","2011-10-25","2011-10-26","2011-10-27","2011-10-28","2011-10-31","2011-11-01","2011-11-02","2011-11-03","2011-11-04","2011-11-07","2011-11-08","2011-11-09","2011-11-10","2011-11-11","2011-11-14","2011-11-15","2011-11-16","2011-11-17","2011-11-18","2011-11-21","2011-11-22","2011-11-23","2011-11-25","2011-11-28","2011-11-29","2011-11-30","2011-12-01","2011-12-02","2011-12-05","2011-12-06","2011-12-07","2011-12-08","2011-12-09","2011-12-12","2011-12-13","2011-12-14","2011-12-15","2011-12-16","2011-12-19","2011-12-20","2011-12-21","2011-12-22","2011-12-23","2011-12-27","2011-12-28","2011-12-29","2011-12-30","2012-01-03","2012-01-04","2012-01-05","2012-01-06","2012-01-09","2012-01-10","2012-01-11","2012-01-12","2012-01-13","2012-01-17","2012-01-18","2012-01-19","2012-01-20","2012-01-23","2012-01-24","2012-01-25","2012-01-26","2012-01-27","2012-01-30","2012-01-31","2012-02-01","2012-02-02","2012-02-03","2012-02-06","2012-02-07","2012-02-08","2012-02-09","2012-02-10","2012-02-13","2012-02-14","2012-02-15","2012-02-16","2012-02-17","2012-02-21","2012-02-22","2012-02-23","2012-02-24","2012-02-27","2012-02-28","2012-02-29","2012-03-01","2012-03-02","2012-03-05","2012-03-06","2012-03-07","2012-03-08","2012-03-09","2012-03-12","2012-03-13","2012-03-14","2012-03-15","2012-03-16","2012-03-19","2012-03-20","2012-03-21","2012-03-22","2012-03-23","2012-03-26","2012-03-27","2012-03-28","2012-03-29","2012-03-30","2012-04-02","2012-04-03","2012-04-04","2012-04-05","2012-04-09","2012-04-10","2012-04-11","2012-04-12","2012-04-13","2012-04-16","2012-04-17","2012-04-18","2012-04-19","2012-04-20","2012-04-23","2012-04-24","2012-04-25","2012-04-26","2012-04-27","2012-04-30","2012-05-01","2012-05-02","2012-05-03","2012-05-04","2012-05-07","2012-05-08","2012-05-09","2012-05-10","2012-05-11","2012-05-14","2012-05-15","2012-05-16","2012-05-17","2012-05-18","2012-05-21","2012-05-22","2012-05-23","2012-05-24","2012-05-25","2012-05-29","2012-05-30","2012-05-31","2012-06-01","2012-06-04","2012-06-05","2012-06-06","2012-06-07","2012-06-08","2012-06-11","2012-06-12","2012-06-13","2012-06-14","2012-06-15","2012-06-18","2012-06-19","2012-06-20","2012-06-21","2012-06-22","2012-06-25","2012-06-26","2012-06-27","2012-06-28","2012-06-29","2012-07-02","2012-07-03","2012-07-05","2012-07-06","2012-07-09","2012-07-10","2012-07-11","2012-07-12","2012-07-13","2012-07-16","2012-07-17","2012-07-18","2012-07-19","2012-07-20","2012-07-23","2012-07-24","2012-07-25","2012-07-26","2012-07-27"};
    static double[] p = new double[]{229.49,231.94,232.97,234,233.36,235.15,235.64,235.78,238.95,242.09,240.61,240.29,237.88,252.11,259.16,263.4,262.1,254.85,254.42,261.27,253.92,259.04,251.58,248.96,239.49,229.39,247.02,249.48,254.9,251.27,246.85,245.43,241.52,231.23,235.67,239.99,238.49,237.41,246.4,249.83,253.67,256.71,255.9,248.93,244.05,242.49,236.52,243.63,246.55,247.3,252.56,259.91,264.41,266.55,262.75,266.33,263.53,261.62,259.38,260.94,249.14,244.63,241.66,240.16,241.81,251.57,251.01,252.49,250.23,244.89,245.79,244.55,243.04,238.84,244.98,247.26,251.91,252.81,252.16,256.83,253.8,251.03,250.19,254.66,254.74,255.76,254.52,252.95,254.57,252.29,243.32,244.88,242.26,240.84,245.05,246.12,243.02,242.79,239.05,233.34,236.22,233.69,234.99,235.84,236.43,243.46,245.25,251.67,250.73,255.7,255.85,256.18,259.71,260.7,262.8,268.98,267.81,275.46,275.98,279.85,280.99,284.3,283.17,278.99,279.48,275.96,274.77,270.99,281.01,281.25,281.28,286,287.25,290.35,291.9,294.01,306.1,309.27,301,302.01,301.02,299.03,300.36,299.59,299.38,296.86,292.72,295.83,300.87,304.21,309.53,308.43,309.87,307.4,309.3,307.96,299.58,298.61,293.31,292.25,299.96,298.31,304.76,300.26,306.16,306.35,308.17,302.61,307.72,309.42,308.73,311.36,309.48,312.2,310.98,311.76,312.84,311.5,311.57,312.43,311.81,313.37,315.3,316.24,314.72,315.77,316.54,316.36,314.78,313.71,320.52,322.2,324.83,324.57,326.89,333.05,332.26,334.97,336.19,338.92,331.3,329.54,323.55,317.75,328.19,332.03,334.41,333.79,326.88,330.01,335.56,334.87,334.01,336.99,342.22,345.45,348.33,344.81,347.06,349.32,350.02,353.16,348.47,340.94,329.32,333.22,333.47,338.6,343.52,339.72,342.46,349.69,350.12,345.61,346,342.8,337.15,342.33,343.86,335.95,320.95,325.46,321.59,329.99,331.84,329.88,335.5,341.89,340.82,341.33,339.06,338.94,335.1,331.83,329.59,328.76,328.8,325.86,321.72,323.28,326.9,323.3,318.47,322.74,328.59,333.01,341.07,343.32,340.8,340.54,337.23,340.52,336.78,338.64,339.98,337.23,337.15,338.06,339.86,337.7,337.06,331.15,324.15,326.91,330.54,331.18,326.02,325.22,323.07,327.54,325.81,328.15,338.28,336.03,336.6,334.01,328.76,322.93,323.12,322.39,316.96,317.64,323.32,317.78,316.24,311.47,306.67,316.37,313.76,322.14,317.39,322.93,326.06,324.87,326.46,333.84,339.84,342.11,347.4,349.84,344.28,344.04,348.19,347.95,354.9,363.54,366.51,376.28,376.66,382.51,387.56,392.34,381.81,381.07,379.76,385.86,378.24,381.8,367.01,363.37,343.52,363.74,353.71,363.44,366.64,372.89,370.04,370,356,346.26,346.66,363.35,365.85,363.46,373.05,379.27,379.29,374.27,370.57,363.78,369.32,373.39,373.6,367.12,369.51,374.06,378.61,382.17,389.51,400.33,402.1,400.83,390.79,393.2,392.1,388.3,386.11,379.85,370.85,364.32,362.28,367.87,367.01,359.65,378.14,389.3,391.15,397.22,410.42,408.46,410.65,387.68,384.46,382.09,394.63,386.85,389.6,393.58,393.84,393.67,385.63,386.5,392.01,389.25,388.76,395.08,384.43,374.65,374.06,368.85,378.16,374.21,367.05,364.65,358.88,366.18,356.92,353.59,365.8,362.96,371.71,377.28,379,382.22,380.22,378.41,379.94,382.82,381.09,378.14,369.75,368.54,370.56,371.72,385.08,385.57,387.61,392.26,395.37,391.59,394,393.88,399.94,402.09,406.56,410.81,410.15,411.62,410.95,409.82,408.29,413.04,417.33,416.01,408.76,415.68,408.87,434.4,432.43,435,440.58,443.95,443.67,442.63,447.06,451.24,455.96,463.6,479.63,479.88,488.81,495.48,484.01,488.43,488.34,500.72,498.96,502.22,508.07,511.33,520.71,527.55,529.53,530.22,518.53,515.71,516.12,527.11,530.21,536.85,552.51,573.4,569.49,569.5,584.6,589.33,585.96,582.89,579.69,590.32,597.61,600.67,593.12,583.09,601.65,612.05,607.17,616.29,618.77,611.19,609.01,605.68,588.62,564.21,592.97,591.64,571.32,557.25,556.01,544.9,593.26,591.02,586.45,567.95,566.15,569.9,565.85,549.74,553.85,552.59,553.56,554.86,551.16,542.9,537.99,531.09,515.57,515.82,545.87,541.68,554.9,549.8,546.86,556.56,563.27,561.87,545.59,548.8,547.38,555.78,556.03,564.39,555.49,560.35,556.46,555.84,558.37,569.7,571.29,569.66,561.81,566.12,555.1,556.33,558.73,553.43,567.97,576.26,582.96,593.2,589.25,597.04,591.52,587.84,582.46,588.37,590.25,590.28,589.62,597.46,587.71,587.26,584.43,559.19,559.1,569.1};
    static int N = 507; // T.length; // 507 is just right before all time high
    static int start = 400;

    public static void main(String[] args) throws Exception {
        System.out.println(T.length == p.length);
        if (args.length>0) start = Integer.parseInt(args[0]);
        if (args.length>1) N = Integer.parseInt(args[1]);

        double tc = decimalYear(T[N+2]); 
        System.out.println("tc: " + tc);
        double m = 0.5;
        double w = 1d;

        // Fitting
        System.out.println("Curve fitting ... ");
        LogPeriodic lppl = new LogPeriodic(tc);

        final CurveFitter fitter = new CurveFitter(new LevenbergMarquardtOptimizer());
        for (int i=start;i<N;i++) fitter.addObservedPoint(decimalYear(T[i]), p[i]);

        double[] initialguess = new double[]{
            0.5,                  // m
            9d                    // w
        };

        double[] fitted = fitter.fit(lppl, initialguess);
        System.out.println("fitted m,w: " + Arrays.toString(fitted));

        // print result as csv
        StringBuilder sb = new StringBuilder();
        for (int i=start;i<T.length;i++) {
            sb.append(T[i]).append(";");
            sb.append(log(p[i])).append(";");
            sb.append(p[i]).append(";");
            sb.append(lppl.value(decimalYear(T[i]), fitted)).append(";");
            sb.append(i>=N?1:0).append(";"); // is predicted price
            sb.append("\n");
        }

        System.out.println(sb.toString());
    }


    public static double decimalYear(String date) {
        Calendar cal = Calendar.getInstance();
        cal.setTime(Date.valueOf(date));
        int doy = cal.get(Calendar.DAY_OF_YEAR);
        int year = cal.get(Calendar.YEAR);
        return doy/365.2425 + year;
        // return doy/(((year & 3) == 0 && ((year % 25) != 0 || (year & 15) == 0)) ? 366 : 365) + year;
    }

    // helper to switch between log and log10 
    public static double log(double d) { return FastMath.log(d); } // Math.log10

    // http://commons.apache.org/proper/commons-math/userguide/linear.html
    public static double[] getABCC(double tc, double m, double w) {
        double sum_fi=0,sum_gi=0,sum_fi2=0,sum_fi_gi=0,sum_gi2=0,sum_yi=0,sum_yi_fi=0,sum_yi_gi=0,sum_hi=0,sum_fi_hi=0,sum_gi_hi=0,sum_hi2=0,sum_yi_hi=0;

        for (int i=start;i<N;i++) {
            double t = decimalYear(T[i]);
            double fi = pow((tc - t),m);
            double gi = fi * cos(w * log(tc - t));
            double hi = fi * sin(w * log(tc - t));
            // double yi = log(p[i]);
            double yi = p[i];

            sum_fi += fi;
            sum_gi += gi;
            sum_fi2 += fi * fi;
            sum_fi_gi += fi * gi;
            sum_gi2 += gi * gi;
            sum_yi += yi;  
            sum_yi_fi += yi * fi;
            sum_yi_gi += yi * gi;
            sum_hi += hi;
            sum_fi_hi += fi * hi;
            sum_gi_hi += gi * hi;
            sum_hi2 += hi * hi;
            sum_yi_hi += yi * hi;
        }

        RealMatrix coefficients = new Array2DRowRealMatrix(new double[][] {
            { N-start+0, sum_fi,     sum_gi,        sum_hi },     // interessanterweise passt der chart nicht so gut bei N+1 übereinander, der zeitpunkt war aber präzierser, ist das bei allen so?
            { sum_fi,    sum_fi2,    sum_fi_gi,     sum_fi_hi }, 
            { sum_gi,    sum_fi_gi,  sum_gi2,       sum_gi_hi },
            { sum_hi,    sum_fi_hi,  sum_gi_hi,     sum_hi2 }},
            false
        );

        DecompositionSolver solver = new LUDecomposition(coefficients).getSolver();

        //Next create a RealVector array to represent the constant vector B and use solve(RealVector) to solve the system
        RealVector constants = new ArrayRealVector(new double[] { sum_yi, sum_yi_fi, sum_yi_gi, sum_yi_hi }, false); // das hinter dem = in der matrix gleichung
        RealVector solution = solver.solve(constants);

        return new double[] {solution.getEntry(0),solution.getEntry(1),solution.getEntry(2),solution.getEntry(3)};
    }

    // 12.6: http://commons.apache.org/proper/commons-math/userguide/optimization.html
    public static class LogPeriodic implements ParametricUnivariateFunction {
        public final double tc;
        public double A,B,C1,C2,w,m,o;
        public double[] best = new double[3];
        public double[] info;

        public LogPeriodic(double tc) {
            this.tc = tc;
        }

        @Override
        public double value(double t, double... parameters) {
            m = parameters[0];
            w = parameters[1];

            // solve the linear parameters
            double[] ABCC = getABCC(tc, m, w);
            A = ABCC[0];
            B = ABCC[1];
            C1 = ABCC[2];
            C2 = ABCC[3];

            double a=pow((tc - t),m);
            return A + B * a + C1 * a * cos(w * log(tc - t)) + C2 * a * sin(w * log(tc - t)); 
        }

        @Override
        public double[] gradient(double t, double... parameters) {
            int params = 2;
            int order = 2;

            DerivativeStructure m = new DerivativeStructure(params,order,0, parameters[0]);
            DerivativeStructure w = new DerivativeStructure(params,order,1, parameters[1]);

            double[] ABCC = getABCC(tc, parameters[0], parameters[1]);
            A = ABCC[0];
            B = ABCC[1];
            C1 = ABCC[2];
            C2 = ABCC[3];

            //    return A + B * pow((tc - t),m) + C1 * pow((tc - t),m) * cos(w * log(tc - t)) + C2 * pow((tc - t),m) * sin(w * log(tc - t)); 
            /*       ==  A 
             *         + B * pow((tc - t),m)
             *         + C1 * pow((tc - t),m) * cos(w * log(tc - t)) 
             *         + C2 * pow((tc - t),m) * sin(w * log(tc - t)); 
             */
            DerivativeStructure f = new DerivativeStructure(params, order, A);
            f = f.add(B).multiply(tc-t).pow(m);
            f = f.add(C1).multiply(tc-t).pow(m).multiply(w.multiply(log(tc-t)).cos());
            f = f.add(C2).multiply(tc-t).pow(m).multiply(w.multiply(log(tc-t)).sin());

            // throws java.lang.ArrayIndexOutOfBoundsException: 2 at org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer.qrDecomposition(LevenbergMarquardtOptimizer.java:862)
            // return f.getAllDerivatives(); 

            // throws unable to perform Q.R decomposition on the 107x2 jacobian matrix
            return new double[] {
                f.getPartialDerivative(1,0),
                f.getPartialDerivative(1,1),
            };
        }
    }
}

EDIT - Stack:

Exception in thread "main" org.apache.commons.math3.exception.ConvergenceException: illegal state: unable to perform Q.R decomposition on the 107x2 jacobian matrix
    at org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer.qrDecomposition(LevenbergMarquardtOptimizer.java:884)
    at org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer.doOptimize(LevenbergMarquardtOptimizer.java:331)
    at org.apache.commons.math3.optim.nonlinear.vector.jacobian.LevenbergMarquardtOptimizer.doOptimize(LevenbergMarquardtOptimizer.java:113)
    at org.apache.commons.math3.optim.BaseOptimizer.optimize(BaseOptimizer.java:143)
    at org.apache.commons.math3.optim.BaseMultivariateOptimizer.optimize(BaseMultivariateOptimizer.java:66)
    at org.apache.commons.math3.optim.nonlinear.vector.MultivariateVectorOptimizer.optimize(MultivariateVectorOptimizer.java:87)
    at org.apache.commons.math3.optim.nonlinear.vector.JacobianMultivariateVectorOptimizer.optimize(JacobianMultivariateVectorOptimizer.java:83)
    at org.apache.commons.math3.optim.nonlinear.vector.jacobian.AbstractLeastSquaresOptimizer.optimize(AbstractLeastSquaresOptimizer.java:197)
    at org.apache.commons.math3.fitting.CurveFitter.fit(CurveFitter.java:172)
    at org.apache.commons.math3.fitting.CurveFitter.fit(CurveFitter.java:136)
    at tryout.Sornette2NoLogPriceLevenberg.main(Sornette2NoLogPriceLevenberg.java:62)
Java Result: 1
Was it helpful?

Solution 2

finally I was able to get things working by using another library: http://finmath.net There is a LevenbergMarquardt class you need to extend. And you do not need to provide any derivatives.

public static class LogPeriodic extends LevenbergMarquardt {
        public final double tc;
        public double A,B,C1,C2,w,m,o;
        public double[] best = new double[3];
        public double[] info;

        public LogPeriodic(int threads, double tc) {
            super(threads);
            this.tc = tc;
        }

        public double value(double t, double... parameters) {
            m = parameters[0];
            w = parameters[1];

            // solve the linear parameters
            double[] ABCC = getABCC(tc, m, w);
            A = ABCC[0];
            B = ABCC[1];
            C1 = ABCC[2];
            C2 = ABCC[3];

            double a=pow((tc - t),m);
            return A + B * a + C1 * a * cos(w * log(tc - t)) + C2 * a * sin(w * log(tc - t)); 
        }

        @Override
        public void setValues(double[] parameters, double[] values) {
            int j=0;
            for (int i=start;i<N;i++) {
                double t = decimalYear(T[i]);
                double a=pow((tc - t),m);
                values[j++] = value(decimalYear(T[i]), parameters);
            }
        }

    }

then you can simply

    LogPeriodic lppl = new LogPeriodic(1,tc);

    int j=0;
    double[] curve = new double[N-start];
    double[] weight = new double[N-start];
    for (int i=start;i<N;i++) {
        curve[j] = p[i];
        weight[j++] = 1;
    }

    double[] initialguess = new double[]{
        0.5,                  // m
        9d                    // w
    };

    lppl.setInitialParameters(initialguess);
    lppl.setWeights(weight);
    lppl.setMaxIteration(100);
    lppl.setTargetValues(curve);
    lppl.setErrorTolerance(1d);
    lppl.run();

    double[] fitted = lppl.getBestFitParameters();
    System.out.println("fitted m,w: " + Arrays.toString(fitted) + " in iterations " + lppl.getIterations());

Outputs:

Curve fitting ... 
fitted m,w: [-0.0895644617344549, 8.162571129643649] in iterations 5

OTHER TIPS

It looks like you have solved your problem, but i'll post my suggestion anyway as it might be usefull to others.

I did get this exception when the jacobian matrix contained infinite or NaN values. And indeed, when trying to run your code, the partial derivatives return NaNs.

There is surely a problem in your derivatives calculation as the code 

DerivativeStructure f = new DerivativeStructure(params, order, A);
f = f.add(B).multiply(tc-t).pow(m);
f = f.add(C1).multiply(tc-t).pow(m).multiply(w.multiply(log(tc-t)).cos());
f = f.add(C2).multiply(tc-t).pow(m).multiply(w.multiply(log(tc-t)).sin());

does not calculate the formula that is in the comments. The order of operations is wrong. For example the second line does ((A+B)*(tc-t))^m not A + B * (tc-t)^m.

Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top