スキップしてメイン コンテンツに移動

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.
  1. Minimize the function call as much as possible: because the function needs heavy calculation
  2. 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)
  1. If one of the edge in the region (might be both edges), just return the edge value as a solution.
  2. Try to apply normal Brent's method (e.g. check the function value of both edges have different sign or not)
  3. If the function value on the edges are different, apply normal Brent's method
  4. 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.
  5. 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.
  6. 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)
  7. 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:
  1. 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.
  2. 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.
  3. 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

コメント

このブログの人気の投稿

Eclipseでコードカバレッジのハイライトを削除する方法

Eclipseには便利なコードカバレッジ表示機能が搭載されていますが、コード内に緑、赤、黄の色付けがされて煩く感じるときもあると思います。 1度カバレッジの色付けが出てしまった後に消す方法の紹介です(方法は簡単)。 下記のキャプチャの青いマーカーで示した「Remove All Sessions」のボタンを押せばすべて消えます。

「特定の文字から始まらない文字列」にマッチする正規表現

「特定の文字から始まらない文字列」 にマッチする正規表現の例です。  以下の例では、Aから始まらない文字列にマッチする正規表現を示しています。 ^(?!A).*$ 私も正規表現の組み方で四苦八苦することがあります。以下の書籍は実践的に様々な正規表現のパターンを例示してくれているので、重宝しています。

ダイソーで買った200円のドライバーセットでHDDを分解

HDDの処分 最近は個人情報の問題もあって、HDDを処分する前にちゃんとデータの消去を気にすることも多くなってきました。消去方法としては大きく分けて下記の3つがあります。 データ消去ソフトでフォーマット HDD内部のプラッタを物理破壊 データ消去を行ってくれる専門の業者や家電量販店(Sofmapやビックカメラで実施していると思います。費用発生。)に持ち込み。 データ消去ソフトでのフォーマットは簡単ですが、欠点として「フォーマットに時間がかかる」「セクタ破損などで中途半端に壊れたディスクのフォーマットができない」などがあります。 またHDD内部のプラッタの物理破壊については、HDDを分解するために、通常のプラスやマイナスドライバーではなく、星形ネジに対応したトルクスドライバーが必要とのこともあって、少し面倒です。 筆者は今回、今後もHDDの廃棄をするだろうなあと思い、思い切って自分で分解して廃棄することにチャレンジしてみました。(家電量販店に持って行くよりも安くできないかというどケチ丸出しですw) HDDの星形ネジ こんなやつです。ちなみに写真はSeagateのST2000DL003というHDDで撮影しました。 トルクスドライバー というわけで、分解のために Amazonでトルクスドライバー を探しました。 調べると T8のもだと使えそう とのことで、いろいろと物色。 セットのものとか T8一本で立派なやつとか 色々あったのですが、HDD壊すだけで800円かぁ(←どケチ)、と思って購入を躊躇。 ネット上で調べると100円ショップのダイソーでも、トルクスドライバーを販売しているとの情報をキャッチ!近所のダイソーに行って、探したところ星形のヘッド交換に対応した精密ドライバーセットがありました。 プラスが10種類、マイナスが8種類、六角が6種類、星形が6種類(今回ほしかったもの)のセットで、何とお値段税抜き200円!、税抜き200円!と安かったので、ダメもとで購入しました。 結論から言うと 買って大正解 でした。 ダイソーの精密ドライバーセット こんな商品です! 星形対応のヘッドを装着するとこんな感じ。ドライバーのグリップもゴムで滑らない様になっていて使いやす

SQLで特定の文字を組み合わせたランダムな文字列を生成

簡易的な方法として「指定した文字列からランダムに1文字選ぶ」を必要な文字の長さ分concat関数でつなげれば実現できます。 1文字ずつ文字を選ぶので、あまり性能もよくない上、セキュリティ的な観点からのランダム性も担保されていないので、あくまで開発中に必要になった時に使う程度が無難だと思います。 下記に英数字大文字小文字を含んだランダムな3文字の文字列を生成するクエリを示します。 # RAND関数で指定した文字列からランダムに1文字選択。 # 下記の例の62の部分はa~z、A~Z、1~9の文字数の合計値を入れた結果 SELECT CONCAT( SUBSTRING('abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ123456789', FLOOR(RAND() * 62 + 1), 1), SUBSTRING('abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ123456789', FLOOR(RAND() * 62 + 1), 1), SUBSTRING('abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ123456789', FLOOR(RAND() * 62 + 1), 1) ) AS random_string;

PHPの配列(array)のメモリ使用量の考察

はじめに 最近PHP上に大量のデータをメモリ上に展開していたのですが、配列(array)の形式(連想配列 or 単純な配列)や配列の要素のデータ構造(数字、配列、文字列など)で大きくメモリ使用量に差が出てくることに気づき、簡単なプログラムを組んで調べてみました。 あくまで筆者の環境での結果なので、細かい数値は参考程度に見てください。 測定環境と方法 OS: Windows 10 PHP 7.4.5 (php-7.4.5-nts-Win32-vc15-x64) 配列に要素を追加するプログラムを書いて、PHPのmemory_get_usage(true)関数を使って実メモリ使用量を計測しました。 計測結果 No. 方式 1MB当たり作成できる 要素数 プログラム 補足 1 キーも値も整数の配列 (整数IDを想定) 28571 // 2,000,000 / 70MB $row = []; for($i = 0; $i < 2000000; $i++) { $row[] = $i; } No.2~6でテストしたプログラム中の要素数は200,000。これだけ一桁多い! 2 キーが文字列、値が整数の連想配列 8333 // 200,000 / 24MB $row = []; for($i = 0; $i < 200000; $i++) { $row[$i.'_key_string'] = $i; } キーの文字列が長い方がメモリ使用量多くなる。 3 キーが整数、値が連想配列の配列 DBから取得してきたデータを想定 2325 // 200,000 / 86MB $row = []; for($i = 0; $i < 200000; $i++) { row[] = ['id' => $i]; } 4 キーが整数、値が連想配列の配列(配列に複数の値を保持) DBから取得してきたデータを想定 2127 // 200,000 /

ADODB.streamオブジェクトを使って文字列とByte配列を相互変換(Excel VBA)

ADODB.streamオブジェクトを使って文字列をByte配列に変換するコードのサンプルです。 ExcelVBAでADODB.streamを使う際には、 1. ExcelのMicrosoft Visual Basic エディタのメニューバーから「ツール->参照設定」とたどる。 2. 表示されたダイアログからMicrosoft ActiveX Data Objectsにチェックを入れる。 という手順が必要です。 文字列からByte配列へ Private Function ADOS_EncodeStringToByte(ByVal cset As String, ByRef strUni As String) As Byte() On Error GoTo e Dim objStm As ADODB.stream: Set objStm = New ADODB.stream objStm.Mode = adModeReadWrite objStm.Open objStm.Type = adTypeText objStm.Charset = cset objStm.WriteText strUni objStm.Position = 0 objStm.Type = adTypeBinary Select Case UCase(cset) Case "UNICODE", "UTF-16" objStm.Position = 2 Case "UTF-8" objStm.Position = 3 End Select ADOS_EncodeStringToByte = objStm.Read() objStm.Close Set objStm = Nothing Exit Function e: Debug.Print "Error occurred while encoding characters" & Err.Description If objStm Is No

MySQLのSQL_CALC_FOUND_ROWS使用上の注意

MySQLのSQL_CALC_FOUND_ROWSを使用する際の無駄なメモリ消費に注意 ページング機能の実装などのために、ヒットした全件数を取得する必要のある場合があるかと思います。 その場合、下記のようにSQL_CALC_FOUND_ROWSを使うと、検索結果に加えて、そのクエリの全ヒット件数が取得できます。 SELECT SQL_CALC_FOUND_ROWS table.column1, table.column2 ...(途中省略)... FROM table WHERE (条件) しかし、SQL_CALC_FOUND_ROWSを使うと、「絞り込んだ結果にヒットする全べての行の結果のセットを作成する」という大きな欠点があります。 これは、LIMIT, OFFSETを指定していても実行されてしまうので、 SELECTで指定するカラム数やデータ量が多い SELECTの結果返ってくる行数が多い 場合、無駄に領域(≒メモリ)を消費してしまうので注意が必要です。 例えば、100万行検索対象としてヒットするクエリの場合、仮にLIMIT 20と指定して最初の20行を取得するようにクエリを書いても、その100万行分の結果セットが作成されてしまうということになります。 対応策 根本的な対応策としては、SQL_CALC_FOUND_ROWSを使わない(結果行数の取得はCOUNTを用いた別クエリで取得する、結果行数をあらかじめサマリーテーブルに保持しておく)ことですが、 SQL_CALC_FOUND_ROWSをどうしても使う必要がある場合は、 SELECT SQL_CALC_FOUND_ROWS primary_id のように最低限のカラムを指定して結果行セットを取得 (LIMIT OFFSET指定前提) 取得したprimary_idを使って必要なデータを取得 して、SQL_CALC_FOUND_ROWSで使用する領域をできるだけ減らすことで、対応するしかないと思います。 世の中ではLIMIT, OFFSETは使わない方がよいとよく書かれていますが、 SQL_CALC_FOUND_ROWSは、書いてしまえばどんなときも検索にヒットする全結果行セットを作成するので、同じくらい使用する際には注意が必要です。

Visual Studio 2010 SP1のアンインストール

Visual Studio 2013に乗り換えるためにVisual Studio 2010をアンインストールしようとしたところで問題発生。。。 先にVisual Studio 2010本体をアンインストールした後、Visual Studio 2010 SP1をアンインストールできなくて困っていました。 Google先生で調べたところ、以下の情報が見つかり、書かれていた通り実施したところ無事Visual Studio 2010 SP1のアンインストールに成功しました。 How to uninstall/remove Visual Studio SP1 アンインストール手順は以下の通りです。 http://www.microsoft.com/en-gb/download/details.aspx?id=23691 からMicrosoft Visual Studio 2010 Service Pack 1 (Installer)をダウンロード VS10sp1-KB983509.exeというファイル名でダウンロードされる(はず)。 コマンドプロンプトから以下のコマンドを実行 (以下の例は、c:\tempにVS10sp1-KB983509.exeがある場合) c:\temp\VS10sp1-KB983509.exe /uninstall /force ダイアログが立ち上がるので、アンインストールを選択して次へ進めばOK!

WinHttp.WinHttpRequestを使ってHttp PostでByte配列やString配列を送信するプログラム(Excel VBA)

WinHttp.WinHttpRequestオブジェクトを使って使ってHttp PostでByte配列のデータを送信するExcel VBAプログラムです。 WinHttp.WinHttpRequestを使う際には、 1. ExcelのMicrosoft Visual Basic エディタのメニューバーから「ツール->参照設定」とたどる。 2. 表示されたダイアログからMicrosoft WinHTTP Serviceにチェックを入れる。 という手順が必要です。 Byte配列をPOST Private Const BOUNDARY As String = "SOMETHING" Private Function httpPostServletByte(url As String, data() As Byte) As Byte() On Error GoTo e Dim WinHttpReq As WinHttp.WinHttpRequest If (WinHttpReq Is Nothing) Then Set WinHttpReq = New WinHttpRequest End If WinHttpReq.Open "POST", url, False WinHttpReq.SetRequestHeader "Content-Type", "application/x-www-form-urlencoded; boundary=" & BOUNDARY WinHttpReq.send data Dim response() As Byte: response = WinHttpReq.ResponseBody httpPostServletByte = response Set WinHttpReq = Nothing Exit Function e: Set WinHttpReq = Nothing Debug.Print "httpPost Error:" & Err.

MySQL: SELECTの結果をUNIONして ORDER BYする際の最適化方法

SELECTの結果をUNIONして ORDER BY する際には下記の点に注意する必要があります。 無駄なメモリ消費 ソートにINDEXが利かない (≒CPU負荷増大) 対応策 可能であればPush-down Condition (各サブクエリ内でORDER BY, LIMIT, OFFSETを適用してからUNION, ORDER BYを実行する)を利用することで、 パフォーマンスを改善できる場合があります。 下記に例を示します。 もともとのクエリ SELECT tmp.* FROM ( SELECT tableA.column1, tableA.column2 FROM tableA WHERE (条件) UNION ALL SELECT tableB.column1, tableB.column2 FROM tableB WHERE (条件) ) AS tmp ORDER BY tmp.column1, tmp.column2 LIMIT 100, 20 Push-down Conditionを用いて書き直したクエリ SELECT tmp.* FROM ( SELECT tableA.column1, tableA.column2 FROM tableA WHERE (条件) ORDER BY tableA.column1, tableA.column2 LIMIT 30 # <- 10 (offset) + 20 (limit) UNION ALL SELECT tableB.column1, tableB.column2 FROM tableB WHERE (条件) ORDER BY tableB.column1, tableB.column2 LIMIT 30 # <- 10 (offset) + 20 (limit) ) AS tmp ORDER BY tmp.column1, tmp.column2 LIMIT 10, 20 ただしこのPush-down Conditionの手法も下記の場合は、効果が半減しますので注意が必要です。 OFFSETの値が大きい場合は、結局全結果セットUNIONと変わらない サブクエリ内のソートで、INDEXが効かない場合