private static double[] newtonSolve(Function<double[], double[]> function, double[] initial, double tolerance) {
        RealVector dx = new ArrayRealVector(initial.length);
        dx.set(tolerance + 1);
        int iterations = 0;
        int d = initial.length;
        double[] values = Arrays.copyOf(initial, initial.length);

        while (dx.getNorm() > tolerance) {
            double[] fx = function.apply(values);
            Array2DRowRealMatrix df = new Array2DRowRealMatrix(fx.length, d);
            ArrayRealVector fxVector = new ArrayRealVector(fx);
            for (int i = 0; i < d; i++) {
                double originalValue = values[i];
                values[i] += H;
                double[] fxi = function.apply(values);
                values[i] = originalValue;
                ArrayRealVector fxiVector = new ArrayRealVector(fxi);
                RealVector result = fxiVector.subtract(fxVector);
                result = result.mapDivide(H);
                df.setColumn(i, result.toArray());
            }
            dx = new RRQRDecomposition(df).getSolver().solve(fxVector.mapMultiply(-1));
            // df has size = initial, and fx has size equal to whatever that function produces.
            for (int i = 0; i < values.length; i++) {
                values[i] += dx.getEntry(i);
            }
            iterations++;
            if (iterations % 100 == 0) {
                tolerance *= 10;
            }
        }
        return values;
    }
private static double [ ] newtonSolve ( Function < double [ ] , double [ ] > function , double [ ] initial , double tolerance ) { RealVector dx = new ArrayRealVector ( initial . length ) ; dx . set ( tolerance + 1 ) ; int iterations = 0 ; int d = initial . length ; double [ ] values = Arrays . copyOf ( initial , initial . length ) ; while ( dx . getNorm ( ) > tolerance ) { double [ ] fx = function . apply ( values ) ; Array2DRowRealMatrix df = new Array2DRowRealMatrix ( fx . length , d ) ; ArrayRealVector fxVector = new ArrayRealVector ( fx ) ; for ( int i = 0 ; i < d ; i ++ ) { double originalValue = values [ i ] ; values [ i ] += H ; double [ ] fxi = function . apply ( values ) ; values [ i ] = originalValue ; ArrayRealVector fxiVector = new ArrayRealVector ( fxi ) ; RealVector result = fxiVector . subtract ( fxVector ) ; result = result . mapDivide ( H ) ; df . setColumn ( i , result . toArray ( ) ) ; } dx = new RRQRDecomposition ( df ) . getSolver ( ) . solve ( fxVector . mapMultiply ( - 1 ) ) ; for ( int i = 0 ; i < values . length ; i ++ ) { values [ i ] += dx . getEntry ( i ) ; } iterations ++ ; if ( iterations % 100 == 0 ) { tolerance *= 10 ; } } return values ; }
Finding the best fit using least-squares method for an equation system

@param function Equation system to find fit for. Input: Parameters, Output: Residuals.
@param initial Initial 'guess' for function parameters
@param tolerance How much the function parameters may change before a solution is accepted
@return The parameters to the function that causes the least residuals