Solving Unimodal Function by Brent's Method

Introduction

I faced a mathematical problem for solving a root from a unimodal function like quadratic function.
Of course we can solve the root if the function is mono decreasing or mono increasing in the target range. However if the function is in potentially there are 2 solutions or no solutions in the target range.
I know this article is not perfect yet. I will keep improving this article :)

Main Requirements

There are mainly 2 requirements when solving the root.

Minimize the function call as much as possible: because the function needs heavy calculation
Robustness against oscillating: Maybe the below figure helps your understanding. If we solve the root strictly, there can be a lot of solutions in the range because of the small oscillating. However we should choose reasonable accrual solution from candidate points.

Solution

I have tried Newton method, but it didn’t work well because of the above bullet point 2.
I have also tried fitting 3 dimensional function (which can have an inflection point), however it needs a lot of calculation iteration to finding out the fitted function.

Finally I decided to use Brent's method - number of times of calling the function is very low and can handle small oscillating very well.
Okay, now only the problem is the function shape.
The normal Brent's method can be applied the range whose both edges have different sign.
That means we cannot directly apply the normal Brent's method to the function we tries to solve.

My rough strategy is finding local extremum point and apply Brent's method in the range between the given edge point and the local extremum points (if the function value of original edges are both same sign).
Full details of my Strategy is described as below. (only explain convex downward case. convex upward case can be treated similarly)

If one of the edge in the region (might be both edges), just return the edge value as a solution.

Try to apply normal Brent's method (e.g. check the function value of both edges have different sign or not)
If the function value on the edges are different, apply normal Brent's method

If the both signs are same, calculate the function value for golden ration point between the edges as a pre-calculation stage of finding out local extremum point. And the function value for the golden ration point is different sign of both edges, apply normal Brent's method between one of the edge point and the golden ration point.

If the function value for the golden ration point has also same sign as the function value for both edge points, trying to find local extremum point.
If the extremum point has same sign as the the function value for both edge points, we can say "No solution". (because we assume unimodal function)

If the local extremum has different sign from the the function value for both edge points, apply normal Brent's method between one of the edge point and the local extremum point. We can select the solution, e.g. lower or upper solution by selecting target range.

Advantages of My strategy:

If the function is unimodal function in the range you try to find the solution, my advanced Brent's method always can find the solution in certain accuracy.
If there is no solution in the range, my advanced Brent's method can detect the solution does not exist in the range with reasonable reason.
Meets 2 requirements I wrote in the starts of this article - Minimize the function call and Robustness

Code

Thank you for reading above looong introduction.. Okay let's show you my code :)

package com.dukesoftware.utils.solve;

import static java.lang.Math.abs;

import com.dukesoftware.utils.math.function.FunctionUtils;
import com.dukesoftware.utils.math.function.Function;

public class AdvancedBrentMethod implements IRootSolver {

    private final SolutionSelector selector;
    private final BrentsMethod solver = new BrentsMethod(100);
    private final LocalExtremumFinder localExtremumFinder = new LocalExtremumFinder();

    public AdvancedBrentMethod(SolutionSelector selector) {
        this.selector = selector;
    }

    @Override
    public double solve(final Function func, double x1, double x2, double tol) {
        final double flower = func.f(x1);
        final double fupper = func.f(x2);

        boolean isLowerSolution = abs(flower) < tol;
        boolean isUpperSolution = abs(fupper) < tol;
        if (isUpperSolution && isLowerSolution) {
            return selector.selectSolution(x1, x2);
        }
        if (isLowerSolution) {
            return x1;
        }
        if (isUpperSolution) {
            return x2;
        }

        final boolean flowerIsPositive = flower > 0;
        final boolean fupperIsPositive = fupper > 0;
        if ((flowerIsPositive && !fupperIsPositive)
                || (!flowerIsPositive && fupperIsPositive)) {
            return solver.solve(func, x1, x2, tol);
        }

        final double intermediate = x1 + LocalExtremumFinder.GOLDEN_RATIO
                * (x2 - x1);
        final double fintermediate = func.f(x1);
        if (abs(fintermediate) < tol) {
            return intermediate;
        }

        Range range = createRange(func, x1, x2, tol, flower, fupper,
                flowerIsPositive, intermediate, fintermediate);

        return solver.solve(func, range.lower, range.upper, tol);
    }

    private Range createRange(final Function func, double x1, double x2,
            double tol, final double flower, final double fupper,
            final boolean flowerIsPositive, final double intermediate,
            final double fintermediate) {
        final boolean fintermediateIsPositive = fintermediate > 0;
        if (flowerIsPositive && !fintermediateIsPositive || !flowerIsPositive
                && fintermediateIsPositive) {
            return selector.createRange(x1, intermediate, x2);
        }
        final double localMin = localExtremumFinder.find(
                x1,
                x2,
                BrentsMethod.EPS,
                tol,
                10,
                reflectOnYaxisIfConvexUpward(func, flower, fupper,
                        flowerIsPositive, fintermediate,
                        fintermediateIsPositive), intermediate, fintermediate);
        return selector.createRange(x1, localMin, x2);
    }

    private static Function reflectOnYaxisIfConvexUpward(final Function func,
            final double flower, final double fupper,
            final boolean flowerIsPositive, final double fintermediate,
            final boolean fintermediateIsPositive) {
        if (flowerIsPositive && fintermediateIsPositive) {
            if (fintermediate > flower && fintermediate > fupper) {
                throw new UnsolveExeption();
            }
            return FunctionUtils.reflectOnYaxis(func);
        }
        if (fintermediate < flower && fintermediate < fupper) {
            throw new UnsolveExeption();
        }
        return func;
    }

    public static enum SolutionSelector {
        LOWER {

            @Override
            double selectSolution(double x1, double x2) {
                return x1;
            }

            @Override
            Range createRange(double x1, double x2, double x3) {
                return new Range(x1, x2);
            }

        },
        LARGER {

            @Override
            double selectSolution(double x1, double x2) {
                return x2;
            }

            @Override
            Range createRange(double x1, double x2, double x3) {
                return new Range(x2, x3);
            }

        };
        abstract double selectSolution(double x1, double x2);

        abstract Range createRange(double x1, double x2, double x3);

    }

    private final static class Range {
        private final double lower;
        private final double upper;

        private Range(double lower, double upper) {
            this.lower = lower;
            this.upper = upper;
        }

    }

}

package com.dukesoftware.utils.solve;

import com.dukesoftware.utils.math.function.Function;

/**
 * I referred C program from http://people.sc.fsu.edu/~jburkardt/m_src/brent/brent.html.
 */
public class LocalExtremumFinder {

    /**
     * the square of the inverse of the golden ratio.
     */
    public static final double GOLDEN_RATIO = 0.5 * (3.0 - Math.sqrt(5.0));

    public double find(final double a, final double b, double eps, double t,
            int maxIter, Function f) {
        final double x = a + GOLDEN_RATIO * (b - a);
        return find(a, b, eps, t, maxIter, f, x, f.f(x));
    }

    public double find(final double a, final double b, double eps, double t,
            int maxIter, Function f, double x, double fx) {
        double d = 0;
        double fu;
        double m;
        double p;
        double q;
        double r;
        double t2;
        double tol;
        double u;

        double sa = a;
        double sb = b;
        double w = x;
        double v = w;
        double e = 0.0;
        double fw = fx;
        double fv = fw;

        for (int i = 0; i < maxIter; i++) {
            m = 0.5 * (sa + sb);
            tol = eps * Math.abs(x) + t;
            t2 = 2.0 * tol;
            /*
             * Check the stopping criterion.
             */
            if (Math.abs(x - m) <= t2 - 0.5 * (sb - sa)) {
                return x;
            }
            /*
             * Fit a parabola.
             */
            r = 0.0;
            q = r;
            p = q;

            if (tol < Math.abs(e)) {
                r = (x - w) * (fx - fv);
                q = (x - v) * (fx - fw);
                p = (x - v) * q - (x - w) * r;
                q = 2.0 * (q - r);
                if (0.0 < q) {
                    p = -p;
                }
                q = Math.abs(q);
                r = e;
                e = d;
            }

            if (Math.abs(p) < Math.abs(0.5 * q * r) && q * (sa - x) < p
                    && p < q * (sb - x)) {
                /*
                 * Take the parabolic interpolation step.
                 */
                d = p / q;
                u = x + d;
                /*
                 * F must not be evaluated too close to A or B.
                 */
                if ((u - sa) < t2 || (sb - u) < t2) {
                    if (x < m) {
                        d = tol;
                    } else {
                        d = -tol;
                    }
                }
            }
            /*
             * A golden-section step.
             */
            else {
                if (x < m) {
                    e = sb - x;
                } else {
                    e = a - x;
                }
                d = GOLDEN_RATIO * e;
            }
            /*
             * F must not be evaluated too close to X.
             */
            if (tol <= Math.abs(d)) {
                u = x + d;
            } else if (0.0 < d) {
                u = x + tol;
            } else {
                u = x - tol;
            }

            fu = f.f(u);
            /*
             * Update A, B, V, W, and X.
             */
            if (fu <= fx) {
                if (u < x) {
                    sb = x;
                } else {
                    sa = x;
                }
                v = w;
                fv = fw;
                w = x;
                fw = fx;
                x = u;
                fx = fu;
            } else {
                if (u < x) {
                    sa = u;
                } else {
                    sb = u;
                }

                if (fu <= fw || w == x) {
                    v = w;
                    fv = fw;
                    w = u;
                    fw = fu;
                } else if (fu <= fv || v == x || v == w) {
                    v = u;
                    fv = fu;
                }
            }
        }
        throw new UnsolveExeption();
    }
}

package com.dukesoftware.utils.math.function;

/**
 * Function Object
 */
public interface Function {
    double f(double x);
}

package com.dukesoftware.utils.math.function;

public class FunctionUtils {
    public static final Function reflectOnYaxis(final Function func) {
        return new Function() {
            @Override
            public double f(double x) {
                return -func.f(x);
            }
        };
    }
}

Usage

public void testAdvancedBrentsMethod() throws Exception {
    final int num = -3;
    double allowDiff = 0.0001;
    RootSolver solver = new AdvancedBrentMethod(SolutionSelector.LOWER);
    double solution = solver.solve(generateQuadraticfunction(num), -1, 1,allowDiff);
    assertEquals(solution. Math.sqrt(-num), allowDiff);
}

public static Function generateQuadraticfunction(final int c) {
   return new Function() {
        @Override
        public double f(double x) {
            return x * x + c;
        }
    };
}

Reference

http://people.sc.fsu.edu/~jburkardt/m_src/brent/brent.html: I refer the code for finding out local extremum. The nicest one.

Information Technology Tips

このブログを検索