オペレーションズ・リサーチ

最急降下法をPythonで【原理から解説】

オペレーションズ・リサーチ

[PR]当サイトはアフィリエイト広告を利用しています。

今回は最急降下法について解説します。

最急降下法は数理最適化において、非線形の無制約最適化問題を解く方法です。最近では機械学習の資格試験の問題で出題されたりするようですね。

なお数理最適化に関する枠組みに関してはこちらをご覧ください。

●エクセルソルバーとは【使い方と数理最適化の基礎】

最適化のテキストによく登場する最急降下法。世に多くの解説がありますが、複雑な内容なので理解が追い付かないという場合もあるかと思います。

今回はなるべくかみ砕いて最急降下法を解説して、学習者の方の手助けをするということを目標にしていきます。最終的に最急降下法を数値微分を用いて汎用的に解くプログラムを作成して簡単な例題を解いてみます。

なお数学的な厳密性は保証していないので、概念の理解ができればよいというところを目指します。

数学的な解説が欲しい時は例えば以下の書籍などを参考にしてください。

文献(1):梅谷俊治, しっかり学ぶ数理最適化 モデルからアルゴリズムまで, 株式会社講談社, 2021.

最急降下法の仕組み

最急降下法のアルゴリズムを示します。

STEP1:初期点\(\boldsymbol{x}_{0}\)を決め、\(k=0\)とする

STEP2:\(\boldsymbol{x}^{(k)}\)が最適解に十分近ければ終了。そうでないならSTEP3へ進む

STEP3:探索方向を最急降下方向\(-\nabla f(\boldsymbol{x}^{(k)})\)とする。またステップ幅\(\alpha^{(k)}\)を以下のように定める。

\(f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)を最小にするような\(\alpha\)。ただし\(\alpha>0\)

STEP4:\(\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}-\alpha^{(k)}\nabla f(\boldsymbol{x}^{(k)})\)とおく

STEP5:\(k=k+1\)としてSTEP2に戻る

このアルゴリズムに対して順に解説していきます。

最急降下方向とは何か

最急降下方向として\(-\nabla f(\boldsymbol{x}^{(k)})\)を採用しました。ようするに関数を偏微分した関数に現在の座標\(\boldsymbol{x}^{(k)}\)を入れた値をマイナスにしたわけです。変数の数を\(p\)とすると以下のように定義できます。

\[\nabla f(\boldsymbol{x})=(\frac{\partial f(\boldsymbol{x})}{\partial x_{1}},\frac{\partial f(\boldsymbol{x})}{\partial x_{2}},\ldots,\frac{\partial f(\boldsymbol{x})}{\partial x_{j}},\ldots,\frac{\partial f(\boldsymbol{x})}{\partial x_{p}})^{T}\]

高校数学の微分積分の範囲で微分係数を求める問題を解いたかと思います。1階の微分を求めて、値が正なら増加、負なら減少というのをやりました。

つまり微分係数が正の範囲なら変数\(x\)は0、1、5と順方向に進めば増加します。負なら0、1、5と順方向に進めば減少します。逆に言うと微分係数が正の範囲で5、1、0と逆方向に進むと関数は減少します。そして負なら5、1、0と進むと増加します。

これは微分係数の符号が変わらなければずっと成立します。

上の移動方向で増加か減少が変わるという話を思い出してください。

  • 微分係数が正なら、逆方向に進めば関数は減少する
  • 微分係数が負なら順方向に進めば減少する

今\(f'(x_{1})=2\)と\(f'(x_{2})=-2\)を考えてみましょう。これに定数\(\alpha=0.2\)を掛けてみましょう。

\(\alpha f'(x_{1})=0.4\)、\(\alpha f'(x_{2})=-0.4\)となります。

つまり微分係数の正負は変わりません。

これを踏まえて\(x^{(k+1)}=x^{(k)}-\alpha^{(k)}f'(x^{(k)})\)を考えます。

\(x_{1}^{(k+1)}=x_{1}^{(k)}-\alpha^{(k)}f'(x_{1}^{(k)})\)

この場合、\(\alpha^{(k)}f'(x_{1}^{(k)})=0.4\)は\(x_{1}\)の順方向への増加です。つまり微分係数が正の時、この逆方向に進めば関数は減少します。

次に\(x_{2}^{(k+1)}=x_{2}^{(k)}-\alpha^{(k)}f'(x_{2}^{(k)})\)の場合。

\(\alpha^{(k)}f'(x_{2}^{(k)})=-0.4\)は\(x_{2}\)の逆方向への移動です。つまり微分係数が負の場合、逆方向5、1、0と進めば関数は増加し、この逆の順方向へ進めば関数は減少します。

つまり微分係数にマイナスを掛けた方向に進めば関数は減少します。これが最急降下方向\(-\nabla f(\boldsymbol{x}^{(k)})\)の意味となります。定数\(\alpha\)が掛けられていますが、正の数なので最急降下方向を変更するものではありません。

移動した点で再び最急降下方向に進めば、移動するたびに目的関数値は小さくなっていきます。

最終的に最急降下方向が0になる。つまり微分した値が0になる点に到達したらアルゴリズムは終了です。そこでは微分した値が0、つまりどの方向に進んでも関数値が変わらないということで、それ以上関数が減少しないからです。それを最適解とします。

なお局所解と大域解という区別もあります。詳しくは以下でまとめています。最急降下法で見つかるのは局所解です。

●エクセルソルバーとは【使い方と数理最適化の基礎】

ステップ幅は加減が大事

上で\(\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}-\alpha^{(k)}\nabla f(\boldsymbol{x}^{(k)})\)と次の探索点を定めました。

単純に考えると、最急降下方向にものすごく進めば、それだけ関数はものすごく減少するように思います。しかしそうはなりません。次の関数を考えます。

\(f(x)=x^{2}\)を考えましょう。\(f'(x)=2x\)ですね。

点\(x_{1}=1\)で\(f(x_{1})=1\)、\(f'(x_{1})=2\)です。普通に考えて微分した値が正なので、逆方向に進めば関数は減少します。\(\alpha=0.2\)とすると、\(\alpha f'(x_{1})=0.4\)なので、

\(x_{1}^{(k+1)}=x_{1}^{(k)}-\alpha^{(k)}f'(x_{1}^{(k)})=1-0.4=0.6\)
となり、\(f(x^{(k+1)})=0.36\)となり確かに減少しました。

しかしながら関数が減少する方向を大きくとる、つまり\(\alpha=5\)とすると、\(\alpha f'(x_{1})=10\)なので

\(x_{1}^{(k+1)}=x_{1}^{(k)}-\alpha^{(k)}f'(x_{1}^{(k)})=1-10=-9\)
となり、\(f(x^{(k+1)})=81\)となり大幅に増加しました。

このように適切にステップ幅\(\alpha\)を選ばないとうまく減少する方向に進みません。

これを定めているのが\(f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)の部分です。

要するにあるステップ幅で進んだ時に一番目的関数値を減少させる(最小にする)ステップ幅を決めなさいという話になります。

ステップ幅は大きすぎると上のようにうまく減少しないです。

またものすごく小さくすれば確かに最急降下方向で目的関数は減少します。微分した位置から少しだけ移動しても微分した値はほとんど変化しないので、正負が変わることはないからです。つまり微分した値が正ならそのまま正だし、負ならそのまま負なので、最急降下方向の考え方が移動先でも成立しています。

しかし小さすぎるステップ幅だと、移動距離が少なすぎて、いつまでたっても関数の値が変化しません。

大きすぎず、小さすぎず、適切なステップ幅を決める必要があります。

\(f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)を最小化するというのは適切なステップ幅を求めるための問題なのです。一番減少する適切なステップ幅を求めるということです。

適切なステップ幅の条件【Armijoの条件】

一般に適切なステップ幅を決めることは容易ではありません。

\(f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)を最小化するのは簡単ではないということです。

そこで次のアルミホ(Armijo)の条件というものが利用されます。

\(g(\alpha)=f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)

ある定数\(0<\tau <1\)に対して

\(g(\alpha)\leq g(0)+\tau g'(0)\alpha\)

を満たす\(\alpha\)を選ぶ。

これの図形的な意味は以下となります。

\(\alpha_{1}\)では青い線より\(g(\alpha)\)が大きいので\(g(\alpha)\leq g(0)+\tau g'(0)\alpha\)は満たしていません。\(\alpha_{2}\)でも同様です。

しかし\(\alpha_{3}\)では満たすのでアルミホの条件が成立します。よって\(\alpha_{3}\)を採用します。

この図形の意味ですが、\(g(0)+\tau g'(0)\alpha\)というのは\(g(0)\)のときの接線の傾きを緩和させた直線です。

つまり最小値を求めたい\(g(\alpha)\)がとりあえず小さくなる方向に\(\alpha\)を進めて、ある程度小さくなるし、それなりに大きく減少する範囲を選んでいます。

\(\tau\)を傾きに掛けることで接線の傾きを緩やかにして、ある程度\(g(\alpha)\)が減少する範囲をサーチライトのように見ているわけです。

途中で\(g(\alpha)\)が増加しても、緩和させた接線以下という条件で増加した部分を除外できます。

上の例では\(\alpha =0\)のときの\(g^{\prime\prime}(\alpha)\)は正でしたが、これが負になるような次のような場合も考えられます。

こういうときにもアルミホの条件は機能していますね。

なお、こういう下のような\(g(\alpha)\)も式だけ見ると考えられそうです。

しかし、探索方向に進むと必ず目的関数は減少するという最急降下方向を選んでいるということを思い出してください。

\(g(\alpha)=f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)

\(\alpha\)の方向に少しだけ進めば必ず目的関数は減少します。つまり上のように最初の\(\alpha =0\)での接線の傾きが正になるような状況を考える必要はありません。

ウルフ条件も備えればより良いが、バックトラック法でなんとかなる

上のアルミホの条件では基本的に小さい\(\alpha\)を選べば条件を満たしやすいです。

しかしながら目的関数を一番減少させるステップ幅という視点では\(\alpha\)はなるべく大きい方が移動量が大きくなって良いわけです。

このような\(\alpha\)を選ぶ条件としてウルフ(Wolfe)の条件というのがあります。

\(g(\alpha)=f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)

ある定数\(0<\tau_{1}<\tau_{2}<1\)に対して

\(g(\alpha)\leq g(0)+\tau_{1}g'(0)\alpha\)

\(g'(\alpha)\geq \tau_{2}g'(0)\)

を満たす\(\alpha\)を選ぶ。

詳細は専門書に譲りますが、要するに\(\alpha\)が小さすぎるとウルフ条件が成立しなくなっています。

しかしながらアルミホの条件とウルフ条件をどちらも満たすように\(\alpha\)を効率的に求めるのは容易ではないのでバックトラック法という方法でアルミホの条件を満たしつつそれなりに大きな\(\alpha\)を求める方法を利用するのが一般的です。

バックトラック法は以下となります。

\(\alpha^{(0)}=1\)とする。

ある定数\(0<\beta<1\)を用いる。

STEP1:\(\alpha^{(i)}\)がアルミホの条件を満たしているか判定する。満たしていたら終了、そうでないなら次の式で\(\alpha\)を更新する。

\(\alpha^{(i+1)}=\beta\alpha^{(i)}\)

STEP2:STEP1へ戻る。

要するに大きな\(\alpha\)からアルミホの条件を満たしているか確認して、満たしていないならどんどん\(\alpha\)を小さくしていくというものです。

ウルフの条件を考えずに大きめの\(\alpha\)が選べるので一般的にはこちらを利用します。

微分【勾配】をどう求めるか

一般にある座標\(\boldsymbol{x}^{(k)}\)における微分の値を求める方法は二つあります。

  • 数式を手計算で微分して新たな関数\(f'(\boldsymbol{x})\)を求めて、それに\(\boldsymbol{x}^{(k)}\)を代入する
  • 数値微分の方法で近似的な\(f'(\boldsymbol{x}^{(k)})\)を得る

数値微分は過去に計算してみたことがあります。

●偏微分を数値計算で求めたい

今回も利用するので、計算式は後の節で掲載します。

関数の形が単純なものなら手計算で勾配を求めれば済みますが、そうではない場合でも汎用的に勾配を求めたい、という需要もあります。

最適化ソルバーはなるべくたくさんの非線形関数を解きたいので、手計算した微分の式をいちいちセッティングしてください、ではあまりよろしくありません。

そこで数値微分を利用します。

数値微分

数値微分での偏導関数の近似値は以下のようになります。前進差分、中心差分、後進差分と三種類あるのですが、精度が良いのは中心差分です。しかしながら計算回数が増えるので、今回は前進差分を採用します。

多変数の偏導関数の差分近似は以下のように定義されます。\(h\)は十分小さい値とします。\(p\)は変数の個数です。\(\boldsymbol{x}=\boldsymbol{a}\)とします。

\[\left. \frac{\partial f(\boldsymbol{x})}{\partial x_{j}} \right|_{\boldsymbol{x}=\boldsymbol{a}}\approx\frac{f(a_{1},a_{2},\ldots,a_{j}+h,a_{j+1},\ldots,a_{p})-f(a_{1},a_{2},\ldots,a_{j},a_{j+1},\ldots,a_{p})}{h}\]

要するにある位置\(\boldsymbol{x}\)で、関数\(f\)の偏微分する変数の座標だけ微小量\(h\)増加させた\(f\)を計算し、そこから元の位置の\(f\)を引きます。その値を微小量\(h\)で割ります。

例えば2変数の場合は以下のようになります。

\[\left. \frac{\partial f(\boldsymbol{x})}{\partial x_{1}} \right|_{\boldsymbol{x}=\boldsymbol{a}}\approx\frac{f(a_{1}+h,a_{2})-f(a_{1},a_{2})}{h}\]
\[\left. \frac{\partial f(\boldsymbol{x})}{\partial x_{2}} \right|_{\boldsymbol{x}=\boldsymbol{a}}\approx\frac{f(a_{1},a_{2}+h)-f(a_{1},a_{2})}{h}\]

最急降下法のプログラム

以上を踏まえて最急降下法のプログラムを作成すると以下のようになります。

\(f(\boldsymbol{x})=x_{1}^{2}+x_{2}^{2}\)

初期点\(\boldsymbol{x}=(10, -10)\)とします。

最適解は\(\boldsymbol{x}^{*}=(0, 0)^{T}\)です。

def calc(x):
    sum = 0
    for i in range(len(x)):
        sum = sum + x[i]*x[i]
    return sum
def bibun(step,x,k):
    high = []
    for i in range(len(x)):
        if(i != k):
            high.append(x[i])
        else:
            high.append(x[i] + step)
    return (calc(high) - calc(x))/(step)
step = 0.0001
e = 0.0001
x = [10, -10]
k = 1000
tau = 0.5
beta = 0.5
for i in range(k):
    judge = False
    for j in range(len(x)):
        if abs(bibun(step, x, j)) > e:
            judge = True
            break
    if judge == False:
        break
    a = 1
    tempx = []
    for j in range(len(x)):
        tempx.append(x[j] - a * bibun(step, x, j))
    tempy = calc(tempx)
    tempgx = []
    for j in range(len(x)):
        tempgx.append(x[j] - step * bibun(step, x, j))
    tempgy = calc(x) + tau * ((calc(tempgx) - calc(x)) / step) * a
    while tempy > tempgy :
        a = a * beta
        for j in range(len(x)):
            tempx[j] = (x[j] - a * bibun(step, x, j))
        tempy = calc(tempx)
        for j in range(len(x)):
            tempgx[j] = (x[j] - step * bibun(step, x, j))
        tempgy = calc(x) + tau * ((calc(tempgx) - calc(x)) / step) * a
    for j in range(len(x)):
        x[j] = x[j] - a * bibun(step, x, j)
print("y:", calc(x))
print("x:",x)    

実行すると

y: 4.9999937219794095e-09

x: [-4.999994530408003e-05, -4.9999991915683495e-05]

きちんと最適解が求められています。

以下プログラムの解説です。

関数の偏導関数を求める

def bibun(step,x,k):
    high = []
    for i in range(len(x)):
        if(i != k):
            high.append(x[i])
        else:
            high.append(x[i] + step)
    return (calc(high) - calc(x))/(step)

何番目の変数かをkに格納して、上で述べた偏微分の定義通りに偏導関数を求める関数です。

最適解かどうかの判定

for i in range(k):
    judge = False
    for j in range(len(x)):
        if abs(bibun(step, x, j)) > e:
            judge = True
            break
    if judge == False:
        break

変数judgeは現在の変数\(\boldsymbol{x}\)での偏導関数の値をそれぞれの変数で計算し、もし偏導関数の大きさが大きいものがあればTrueになります。すべての偏導関数の値が十分小さければFalseのままで、Falseならif文により繰り返しを抜けて探索終了です。

g(α)を求める

アルミホの条件を計算するため\(g(\alpha)\)を以下のように求めます。

    a = 1
    tempx = []
    for j in range(len(x)):
        tempx.append(x[j] - a * bibun(step, x, j))
    tempy = calc(tempx)

\(g(\alpha)=f(\boldsymbol{x}^{(k)}-\alpha\nabla f(\boldsymbol{x}^{(k)}))\)

この部分の計算です。

g(0)+τg'(0)αを求める

続いて\(g(0)+\tau g'(0)\alpha\)を計算します。

    tempgx = []
    for j in range(len(x)):
        tempgx.append(x[j] - step * bibun(step, x, j))
    tempgy = calc(x) + tau * ((calc(tempgx) - calc(x)) / step) * a

\(g(0)=f(\boldsymbol{x}^{(k)})\)

なのでcalc(x)に相当します。

\(g'(0)=\frac{g(0+h)-g(0)}{h}\)

\(g(0+h)=f(\boldsymbol{x}^{(k)}-h\nabla f(\boldsymbol{x}^{(k)}))\)

なので、\(\boldsymbol{x}^{(k)}-h\nabla f(\boldsymbol{x}^{(k)})\)の部分がtempgxとなります。

そして((calc(tempgx) – calc(x)) / step)の部分で\(g'(0)\)を計算しています。

アルミホの条件を満たすまでバックトラック法を行う

以下の部分でバックトラック法を行い\(\alpha\)を求めます。

    while tempy > tempgy :
        a = a * beta
        for j in range(len(x)):
            tempx[j] = (x[j] - a * bibun(step, x, j))
        tempy = calc(tempx)
        for j in range(len(x)):
            tempgx[j] = (x[j] - step * bibun(step, x, j))
        tempgy = calc(x) + tau * ((calc(tempgx) - calc(x)) / step) * a

\(\alpha\)を\(\beta\)倍して小さくし、アルミホの条件に出てくる左辺と右辺を計算します。その後条件をwhile文のところで判定して、満たしていないなら繰り返して\(\alpha\)を小さくしていきます。

探索点の位置を更新して繰り返す

探索点を更新して繰り返しの最初に戻ります。

    for j in range(len(x)):
        x[j] = x[j] - a * bibun(step, x, j)

他の例題で解いてみる

最初に挙げた文献(1)にある次の問題もプログラム中のcalcを変更して解いてみました。

  • 目的関数:
    \(f(\boldsymbol{x})=(x_{1}-2)^{4}+(x_{1}-2x_{2})^{2}\)
    →最小化
  • 変数:\(x_{1}, x_{2}\)
  • 初期点:\(\boldsymbol{x}^{0}=(0, 3)^{T}\)
  • 最適解:\(\boldsymbol{x}^{*}=(2, 1)^{T}\)
def calc(x):
    return (x[0] - 2)**4 + (x[0] - 2 * x[1])**2

出力は以下となりました。

y: 6.171466812735292e-08

x: [1.9849199272911, 0.99240996364555]

ほぼ最適解ですね。

まとめ【最急降下法はアルミホの条件と数値微分がカギ】

今回は最急降下法について理論から実際のプログラムの作成まで解説しました。

  • アルミホの条件
  • 数値微分

ここが理解できるかというのがカギです。

数式で解説していったので、難解かもしれませんが、一変数で考えれば高校数学くらいの知識でなんとかなるのではないでしょうか。

今回の解説が資格試験や大学の講義で最急降下法を学ぶ必要のある方の参考になりましたら幸いです。

一変数で考えて具体的な数値を入れてみるのも理解を助けます

当ブログ(シルルスのコードおきば)では数理最適化関係の記事を他にも執筆しています。参考になりましたら幸いです。

●エクセルソルバーとは【使い方と数理最適化の基礎】

●偏微分を数値計算で求めたい

●ソルバーをVBAで繰り返す【初期値を変えて複数解を自動探索】

●エクセルソルバーで複数解を求める方法【初期値を変えるしかないです】

●エクセルソルバーと組み合わせ最適化問題【解法例と基礎知識】

●ソルバーとユーザー定義関数の連携


USBメモリ 128GB USB2.0 日本製【翌日配達送料無料】 KIOXIA(旧東芝メモリー)TransMemory U202 キャップ式 ホワイト 海外パッケージ LU202W128GG4