public static double minimize(DifferentiableMultivariateFunction func, double[] x, double gtol, int maxIter) {
        if (gtol <= 0.0) {
            throw new IllegalArgumentException("Invalid gradient tolerance: " + gtol);
        }

        if (maxIter <= 0) {
            throw new IllegalArgumentException("Invalid maximum number of iterations: " + maxIter);
        }

        double den, fac, fad, fae, sumdg, sumxi, temp, test;

        int n = x.length;
        double[] dg = new double[n];
        double[] g = new double[n];
        double[] hdg = new double[n];
        double[] xnew = new double[n];
        double[] xi = new double[n];
        double[][] hessin = new double[n][n];

        // Calculate starting function value and gradient and initialize the
        // inverse Hessian to the unit matrix.
        double f = func.g(x, g);
        logger.info(String.format("BFGS: initial function value: %.5f", f));

        for (int i = 0; i < n; i++) {
            hessin[i][i] = 1.0;
            // Initialize line direction.
            xi[i] = -g[i];
        }

        double stpmax = STPMX * max(norm(x), n);

        for (int iter = 1; iter <= maxIter; iter++) {
            // The new function evaluation occurs in line search.
            f = linesearch(func, x, f, g, xi, xnew, stpmax);

            if (iter % 100 == 0) {
                logger.info(String.format("BFGS: the function value after %3d iterations: %.5f", iter, f));
            }

            // update the line direction and current point.
            for (int i = 0; i < n; i++) {
                xi[i] = xnew[i] - x[i];
                x[i] = xnew[i];
            }

            // Test for convergence on x.
            test = 0.0;
            for (int i = 0; i < n; i++) {
                temp = abs(xi[i]) / max(abs(x[i]), 1.0);
                if (temp > test) {
                    test = temp;
                }
            }

            if (test < TOLX) {
                logger.info(String.format("BFGS converges on x after %d iterations: %.5f", iter, f));
                return f;
            }

            System.arraycopy(g, 0, dg, 0, n);

            func.g(x, g);

            // Test for convergence on zero gradient.
            den = max(f, 1.0);
            test = 0.0;
            for (int i = 0; i < n; i++) {
                temp = abs(g[i]) * max(abs(x[i]), 1.0) / den;
                if (temp > test) {
                    test = temp;
                }
            }

            if (test < gtol) {
                logger.info(String.format("BFGS converges on gradient after %d iterations: %.5f", iter, f));
                return f;
            }

            for (int i = 0; i < n; i++) {
                dg[i] = g[i] - dg[i];
            }

            for (int i = 0; i < n; i++) {
                hdg[i] = 0.0;
                for (int j = 0; j < n; j++) {
                    hdg[i] += hessin[i][j] * dg[j];
                }
            }

            fac = fae = sumdg = sumxi = 0.0;
            for (int i = 0; i < n; i++) {
                fac += dg[i] * xi[i];
                fae += dg[i] * hdg[i];
                sumdg += dg[i] * dg[i];
                sumxi += xi[i] * xi[i];
            }

            // Skip update if fac is not sufficiently positive.
            if (fac > sqrt(EPSILON * sumdg * sumxi)) {
                fac = 1.0 / fac;
                fad = 1.0 / fae;

                // The vector that makes BFGS different from DFP.
                for (int i = 0; i < n; i++) {
                    dg[i] = fac * xi[i] - fad * hdg[i];
                }

                // BFGS updating formula.
                for (int i = 0; i < n; i++) {
                    for (int j = i; j < n; j++) {
                        hessin[i][j] += fac * xi[i] * xi[j] - fad * hdg[i] * hdg[j] + fae * dg[i] * dg[j];
                        hessin[j][i] = hessin[i][j];
                    }
                }
            }

            // Calculate the next direction to go.
            Arrays.fill(xi, 0.0);
            for (int i = 0; i < n; i++) {
                for (int j = 0; j < n; j++) {
                    xi[i] -= hessin[i][j] * g[j];
                }
            }
        }

        logger.warn(String.format("BFGS reaches maximum %d iterations: %.5f", maxIter, f));
        return f;

    }
@Test
    public void testLBFGSB3() {
        System.out.println("L-BFGS-B: x_0 ^ 4 + (x_0 - 3)^2 + (x_1 - 4)^2 + 1");

        DifferentiableMultivariateFunction func = new DifferentiableMultivariateFunction() {
            @Override
            public double f(double[] x) {
                double x0 = x[0] - 3;
                double x1 = x[1] - 4;
                return Math.pow(x[0], 4) + x0 * x0 + x1 * x1 + 1;
            }

            @Override
            public double g(double[] x, double[] g) {
                double x0 = x[0] - 3;
                double x1 = x[1] - 4;
                g[0] = 4 * Math.pow(x[0], 3) + 2 * x0;
                g[1] = 2 * x1;
                return Math.pow(x[0], 4) + x0 * x0 + x1 * x1 + 1;
            }
        };

        double[] x = new double[2];
        double[] l = {3.5, 3.5};
        double[] u = {5.0, 5.0};

        double result = BFGS.minimize(func, 5, x, l, u, 1E-8, 500);
        assertEquals(3.5, x[0], 1E-15);
        assertEquals(4.0, x[1], 1E-15);
        assertEquals(151.3125, result, 1E-15);
    }