[最优化技术] 3-7 坐标轮换法 – H_Elden

座標回転法

簡単な説明

\(\qquad\)前回の記事では、修正ニュートン法を紹介しました。この記事では、勾配情報を必要としない直接検索方法について説明します。座標降下法。変数回転法、交互法、次元削減法とも呼ばれます。

\(\qquad\)最急降下法と同様、これも反復探索法です。違いは、最急降下法では負の勾配方向が使用されるのに対し、座標回転法では負の勾配方向が使用されることです。座標軸方向(基底ベクトル) を検索方向として指定します。

\(\qquad\)この方法の基本的な考え方は、一度に 1 つの座標軸方向に沿った 1 次元検索のみを実行することです。その方向の最小値点を見つけたら、次の座標軸方向に切り替えて探索します。収束条件が満たされるまで、各座標軸方向に沿って順番に探索します。この方法は、多次元の問題を一連の 1 次元の問題に変換するため、「次元削減法」とも呼ばれます。

原理

のために \(n\) メタファンクション \(f(X)\)探索方向は座標系の基底ベクトルとして取得されます。 \(e_i\)：

\[S_i^{(k)} = e_i = (0, \dots, 1, \dots, 0)^T \quad (1 \text{ 在第 } i \text{ 个位置})
\]

1 回の反復で、アルゴリズムは次のように順次実行されます。 \(e_1, e_2, \dots, e_n\) 方向 \(n\) 一次元検索。いいえ。 \(k\) 車輪 \（私\） 検索の反復は次のとおりです。

\[X_i^{(k)} = X_{i-1}^{(k)} + \alpha_i^{(k)} S_i^{(k)}, \quad i=1,2,\dots,n
\]

で \(\alpha_i^{(k)}\) エッジ用 \(S_i^{(k)}\) 方向に 1 次元探索を実行することで得られる最適なステップ サイズ。

注: 座標回転法の探索方向は固定です。\(S_i^{(k)}\) と反復ラウンド数 \(k\) 何もすることはありません。

具体的な手順

与えられた初期点に対して \(X_0^{(0)}\) と収束精度 \(\バレプシロン\)：

① 始点を決める \(X_0^{(0)}\)収束精度 \(\バレプシロン\)ラウンド数を注文する \(k=1\)；

② 基底ベクトルを探索方向とし、\(S_i^{(k)} = e_i, \; i=1, 2, \dots, n\)(二項関数の場合 \(n=2\)）；

③ 現時点から \(X_0^{(k)}\) 出発する：

• 平行 \(S_1^{(k)}\) 点を取得する方向に一次元探索を実行します。 \(X_1^{(k)}\)；
• から \(X_1^{(k)}\) 出発、一緒に \(S_2^{(k)}\) 点を取得する方向に一次元探索を実行します。 \(X_2^{(k)}\)；
• \(\ドット\)
• から \(X_{n-1}^{(k)}\) 出発、一緒に \(S_n^{(k)}\) このラウンドの終点を取得する方向に 1 次元検索を実行します。 \(X_n^{(k)}\)。

今すぐ：\(X_0^{(k)} \xrightarrow{S_1} X_1^{(k)} \xrightarrow{S_2} X_2^{(k)} \dots \xrightarrow{S_n} X_n^{(k)}\)。

④収束判定：if \(\|X_n^{(k)} – X_0^{(k)}\| \le \varepsilon\)それで終わり、最適解は \(X_n^{(k)}\);それ以外の場合は、 \(X_0^{(k+1)} = X_n^{(k)}\)、\(k=k+1\)手順③を繰り返します。

記号の説明：

• \(X_i^{(k)}\) 意味：最初に \(k\) 反復ラウンドでは、 \（私\） 座標軸方向に一次元探索して得られる点。どちらも初めてです \（私\） 探索の終点もまた、 \(i+1\) 検索の開始点。
• のために \(n\) メタ関数、反復の各ラウンドに必要 \(n\) 一次元検索。

ブロック図

記号の説明: 写真内 \(x\)、\(d_i^k\)、\(\alpha_i^k\) それぞれ本文に対応 \(X\)、\(S_i^{(k)}\)、\（\アルファ\）;写真のラウンド数はから始まります。 \(k=0\) 計算から始めると、基準は次のとおりです。 \(<\バレプシロン\); からのテキスト \(k=1\) 計算から始めると、基準は次のとおりです。 \(\le\バレプシロン\)。

計算例

ここでは、前の例を使用して最初の 2 ラウンドの反復を示し、収束プロセスを示します。

トピック

\[f(X) = x_1^2 + 2x_2^2 – 4x_1 – 2x_1 x_2,\;\; X^{(0)} = [1,1]^T,\;\バレプシロン = 0.01 \]

覚えて \(X_0^{(1)} = X^{(0)}\)。

答え

記述の便宜上、各 1 次元検索で最適化された変数は次のように一律に記録されます。 \（\アルファ\）(この検索で​​は自由座標であり、上記の最適なステップ サイズではありません)、対応する 1 変数関数は次のように記録されます。 \(\varphi(\alpha)\)。

反復 1 (\(k=1\))：

初期点 \(X_0^{(1)} = [1, 1]^T\)。

1) 沿って \(x_1\) 軸方向検索（\(S_1 = [1, 0]^T\))：

修理済み \(x_2 = 1\)関数は約になります \（\アルファ\） 次の単項関数:

\[\varphi(\alpha) = f(\alpha, 1) = \alpha^2 + 2(1)^2 – 4\alpha – 2\alpha(1) = \alpha^2 – 6\alpha + 2
\]

微分を求めて 0 とします。\(\varphi'(\alpha) = 2\alpha – 6 = 0 \implies \alpha = 3\)。
新しいポイントを獲得する \(X_1^{(1)} = [3, 1]^T\)。

2) 沿って \(x_2\) 軸方向検索（\(S_2 = [0, 1]^T\))：

から \(X_1^{(1)} = [3, 1]^T\) オフにする、固定する \(x_1 = 3\)も使用します \（\アルファ\） 現在の自由変数を表します。

\[\varphi(\alpha) = f(3, \alpha) = 3^2 + 2\alpha^2 – 4(3) – 2(3)\alpha = 2\alpha^2 – 6\alpha – 3
\]

微分を求めて 0 とします。\(\varphi'(\alpha) = 4\alpha – 6 = 0 \implies \alpha = 1.5\)。
このラウンドのエンドポイントを取得する \(X_2^{(1)} = [3, 1.5]^T\)。

3) 収束判定：

\[\|X_2^{(1)} – X_0^{(1)}\| = \|[3, 1.5]^T – [1, 1]^T\| = \|[2, 0.5]^T\| = \sqrt{4 + 0.25} \約 2.06 > \バレプシロン \]

収束条件が満たされない場合は、次のラウンドに進みます。

2 回目の反復 (\(k=2\))：

から \(X_0^{(2)} = X_2^{(1)} = [3, 1.5]^T\) 出発する。

1) 沿って \(x_1\) 軸方向検索（\(S_1 = [1, 0]^T\))：

修理済み \(x_2 = 1.5\)：

\[\varphi(\alpha) = f(\alpha, 1.5) = \alpha^2 + 2(1.5)^2 – 4\alpha – 2\alpha(1.5) = \alpha^2 – 7\alpha + 4.5
\]

導関数を求めます \(\varphi'(\alpha) = 2\alpha – 7 = 0 \implies \alpha = 3.5\)、それで \(X_1^{(2)} = [3.5, 1.5]^T\)。

2) 沿って \(x_2\) 軸方向検索（\(S_2 = [0, 1]^T\))：

修理済み \(x_1 = 3.5\)：

\[\varphi(\alpha) = f(3.5, \alpha) = 2\alpha^2 – 7\alpha – 1.75
\]

導関数を求めます \(\varphi'(\alpha) = 4\alpha – 7 = 0 \implies \alpha = 1.75\)、それで \(X_2^{(2)} = [3.5, 1.75]^T\)。

3) 収束判定：

\[\|X_2^{(2)} – X_0^{(2)}\| = \|[0.5, 0.25]^T\| = \sqrt{0.3125} \約 0.56 > \バレプシロン \]

それでも満足できないので、繰り返し続けます。

正確な 1 次元検索では、各ラウンドの結果は次のようになります (ラウンド 1 と 2 の計算については上記を参照)。

第8ラウンドの変位 \(0.008735 \le \バレプシロン\)反復を停止して取得します

\[X^* \approx [3.992188,\ 1.996094]^T,\qquad f(X^*) \およそ -7.999969 \]

すでに理論上の最適解に非常に近づいています \([4,2]^T\) そして \(f^*=-8\)。

メソッド分析

1. シンプルな方向、グラデーションは不要: 外側の検索には関数値のみが必要であり、偏導関数を計算する必要はありません。 1 次元検索では、関数の値のみに依存する黄金分割法などの方法を使用できます。これは、導関数を見つけるのが難しい関数に非常に適しています。

2. 多次元問題は効率が低く、収束が遅い: 一度に 1 つの変数のみが変更され、他の変数は固定され、検索パスは直角にジグザグになります。図 4-21(b) に示すように、等値線が長軸と短軸が座標軸に平行でない楕円である場合、複数回の反復が必要となり、効率が大幅に低下します。

3. 特別な状況では、1ラウンドが実施されます：等高線の長軸と短軸が座標軸に平行な同心円または楕円、つまり次のような形状の場合 \(F(X) = ax_1^2 + bx_2^2 + \text{線形項}\) 関数 (変数間の交差項なし) \(x_1 x_2\))、反復は 1 ラウンドのみ (\(n\) 探索）極限点に到達する。バイナリ関数の場合、図 4-21(a) に示すように、2 つの 1 次元検索になります。このセクションの例には以下が含まれることに注意してください。 \(x_1 x_2\) 交差項と等高線は斜めの楕円であり、(b) の場合に属するため、複数回の反復が必要です。

4. 故障状態（リッジの問題）：座標軸と斜めに交わる等高線上に「尾根」が現れる場合、座標軸方向の探索は稜線上で止まり、尾根方向に沿って探索を続けることができないため、最適な点を見つけることができません。このとき、図 4-21(c) に示すように、軸変換（座標変換）が必要になります。

コード例

コアループのみが以下に残ります。完全なプログラムの関数定義、一次元検索サブルーチン、および出力を以下に示します。

// 核心循环：每轮依次沿 n 个坐标轴方向作一维搜索
for (int k = 0; k < max_rounds; ++k) {
    vector x_start = x;

    for (int i = 0; i < n; ++i) {
        vector d(n, 0.0);
        d[i] = 1.0;

        auto phi = [&](double alpha) {
            vector new_x(x.size());
            for (size_t j = 0; j < x.size(); ++j)
                new_x[j] = x[j] + alpha * d[j];
            return func(new_x);
        };

        // 省略异常处理，完整写法见下
        auto bracket = bracket_minimum(phi, 0.0, 0.1);
        double alpha_star = golden_section(phi, bracket.first, bracket.second, eps);

        for (size_t j = 0; j < x.size(); ++j)
            x[j] += alpha_star * d[j];
    }

    double dist = 0.0;
    for (int i = 0; i < n; ++i)
        dist += (x[i] - x_start[i]) * (x[i] - x_start[i]);
    if (sqrt(dist) <= eps) break;
}

完全なプログラム

#include 
#include 
#include 
#include 
#include 
#include 

using namespace std;

// 进退法求初始区间 (复用)
pair bracket_minimum(function func, double a0, double h) {
    const int MAX_ITER = 1e4;
    if (h <= 0) throw invalid_argument("Step size h must be positive.");
    if (func(a0) <= func(a0 + h)) {
        double x0 = a0 + h, x1 = a0;
        h *= 2;
        double x2 = x1 - h;
        int iter = 0;
        while (func(x2) < func(x1)) {
            if (iter++ > MAX_ITER) throw runtime_error("Max iterations exceeded");
            x0 = x1; x1 = x2; h *= 2; x2 = x1 - h;
        }
        return make_pair(x2, x0);
    } else {
        double x0 = a0;
        double x1 = a0 + h;
        h *= 2;
        double x2 = x1 + h;
        int iter = 0;
        while (func(x2) < func(x1)) {
            if (iter++ > MAX_ITER) throw runtime_error("Max iterations exceeded");
            x0 = x1; x1 = x2; h *= 2; x2 = x1 + h;
        }
        return make_pair(x0, x2);
    }
}

// 黄金分割法 (复用)
double golden_section(function func, double l, double r, double eps) {
    double a = l, b = r, aa, bb;
    while (b - a > eps) {
        aa = b - 0.618 * (b - a);
        bb = a + 0.618 * (b - a);
        if (func(aa) < func(bb)) b = bb;
        else a = aa;
    }
    return (a + b) / 2;
}

// 坐标轮换法
vector coordinate_descent(
    function&)> func,
    const vector& x0,
    double eps,
    int max_rounds = 1000
) {
    vector x = x0;
    int n = x0.size();
    
    cout << "\n迭代过程:" << endl;
    cout << "轮数\t点 X\t\t\t函数值 f(X)" << endl;
    cout << "--------------------------------------------------------" << endl;

    for (int k = 0; k < max_rounds; ++k) {
        vector x_start = x;

        for (int i = 0; i < n; ++i) {
            vector d(n, 0.0);
            d[i] = 1.0;

            auto phi = [&](double alpha) {
                vector new_x(x.size());
                for (size_t j = 0; j < x.size(); ++j) new_x[j] = x[j] + alpha * d[j];
                return func(new_x);
            };

            pair bracket;
            try { bracket = bracket_minimum(phi, 0.0, 0.1); } 
            catch (...) { bracket = make_pair(0.0, 10.0); }

            double alpha_star = golden_section(phi, bracket.first, bracket.second, eps);

            for (size_t j = 0; j < x.size(); ++j) x[j] += alpha_star * d[j];
        }

        cout << k + 1 << "\t";
        for (size_t i = 0; i < x.size(); ++i) {
            printf("%.6lf", x[i]);
            cout << (i == x.size() - 1 ? "" : ", ");
        }
        printf("\t%.6lf\n", func(x));

        double dist = 0.0;
        for (int i = 0; i < n; ++i) dist += (x[i] - x_start[i]) * (x[i] - x_start[i]);
        if (sqrt(dist) <= eps) {
            cout << "收敛！" << endl;
            break;
        }
    }
    return x;
}

int main() {
    cout << "坐标轮换法求解多元函数极小值" << endl;
    cout << "目标函数: f(X) = x1^2 + 2*x2^2 - 4*x1 - 2*x1*x2" << endl;
    cout << "理论最优解: X* = [4.0, 2.0]^T, f(X*) = -8.0" << endl << endl;

    auto f = [](vector x) {
        return x[0] * x[0] + 2 * x[1] * x[1] - 4 * x[0] - 2 * x[0] * x[1];
    };

    vector x0 = {1.0, 1.0};
    double eps = 0.01;

    vector result = coordinate_descent(f, x0, eps);

    cout << "\n最终结果:" << endl;
    cout << "最优解 X* = [";
    for (size_t i = 0; i < result.size(); ++i) cout << result[i] << (i == result.size() - 1 ? "" : ", ");
    cout << "]^T" << endl;
    cout << "目标函数值 f(X*) = " << f(result) << endl;

    return 0;
}

\(\qquad\)プログラムは黄金分割法を使用します (\(\バレプシロン=0.01\)) 1 次元検索の場合、実際の測定の 1 回目は次のとおりです。 \((3.002196,1.501746)\)手で計算した \((3,1.5)\) 少し異なります。ラウンド 8 は次のように収束します。 \((3.995681,1.997798)\)、\(f=-7.999991\)最終的なホイール変位は約 \(0.0097 \the 0.01\)。この違いは、1 次元検索の精度によって生じます。

要約する

座標回転法は、最も直接的な次元削減の最適化法です。

アドバンテージ：

• 原理は非常にシンプルで、勾配を計算する必要はありません。
• それぞれの 1 次元検索で所定の精度を達成できます。解析手法を使えば理論的には正確に解くことができます。
• 等値線が長軸と短軸が座標軸に平行な円または楕円である問題の場合、効率は非常に高く、2 回の探索で最適な点に到達できます。

欠点がある：

• 収束速度は通常遅く、等高線の形状に強く影響されます。
• 輪郭が斜めの楕円である場合、複数回の反復が必要となり、効率が大幅に低下します。
• 尾根が現れると探索が尾根上で停止し、最適な点を見つけることができず、軸変換が必要になります。

\(\qquad\)座標回転の次元ごとの検索の考え方は、その後の方法の基礎にもなります。パウエルの方法では、座標回転の原理を使用して検索方向を構築します。

\(\qquad\)次の記事で紹介します共役方向法。座標回転法と最急降下法から始まり、二次関数を最大化するための共役方向のセットを構築します。 \(n\) 段階的に収束する。次の記事で紹介します共役勾配法。

この研究ノートの参考資料:

[1] 白慶春、孫敬民、梁英春。機械最適化設計 第 7 版[M]。北京: Machinery Industry Press、2024 年。ISBN: 978-7-111-75103-8 オンライン リンク

[2] 武漢理工大学「最適化技術B」コースウェア、講師：Yan Bin氏