【最適化問題の基礎】ニュートン法とヘッセ行列

前回は最も単純な連続最適化の手法の一つである「最急降下法」について解説しました。続いて、本稿では目的関数の勾配の勾配（2次微分）「ヘッセ行列」の情報を使って最適化する「ニュートン法」について解説します。

---

 

 ニュートン法は2次微分を使う

「最急降下法」では目的関数の勾配、すなわち1次微分を用いていました。これは、目的関数をある点において1次近似する、という数学的な意味があります。多次元になると分かりづらいですが、要するに最急降下法は接線を求める操作とほとんど同じことをして最適化しているのです。

図．最急降下法のイメージ

では、目的関数を1次関数ではなく2次関数で近似することもできるのではないか、という考えに至るのはごく自然です。

関数上のある点における値と勾配（1次微分）の他に、2次微分の情報があれば、その関数に接する2次関数が得られます。ニュートン法ではこうして得られる2次関数（2次近似）を利用して最適化を行う手法です。

図．ニュートン法のイメージ

※「ニュートン法」は「ニュートン・ラフソン法」とも呼ばれることがあります。

 

 ヘッセ行列（ヘシアン）とは？

多変数関数の場合は、$n$変数関数に対して2次微分を成分とする $n \times n$ の正方行列を考えます。これを「ヘッセ行列」または「ヘシアン」（Hessian）と言います。$$\boldsymbol{H} \equiv \left[\begin{array}{cccc}
\dfrac{\partial^{2} f}{\partial x_{1}^{2}} & \dfrac{\partial^{2} f}{\partial x_{1} \partial x_{2}} & \cdots & \dfrac{\partial^{2} f}{\partial x_{1} \partial x_{n}} \\
\dfrac{\partial^{2} f}{\partial x_{2} \partial x_{1}} & \dfrac{\partial^{2} f}{\partial x_{2}^{2}} & \cdots & \dfrac{\partial^{2} f}{\partial x_{2} \partial x_{n}} \\
\vdots & \vdots & \ddots & \vdots \\
\dfrac{\partial^{2} f}{\partial x_{n} \partial x_{1}} & \dfrac{\partial^{2} f}{\partial x_{n} \partial x_{2}} & \cdots & \dfrac{\partial^{2} f}{\partial x_{n}^{2}}
\end{array}\right]$$ヘッセ行列は2次微分であることを明示して$\nabla^2 f(\boldsymbol{x})$と表記されることもあります。ここで、変数$x_{i}$で偏微分してから変数$x_{j}$で偏微分した第2次導関数と、変数$x_{j}$で偏微分してから変数$x_{i}$で偏微分した第2次導関数は一致します。つまり、$$\dfrac{\partial^{2} f}{\partial x_{i} \partial x_{j}}=\dfrac{\partial^{2} f}{\partial x_{j} \partial x_{i}}$$が成り立つので、ヘッセ行列は常に$n$次の対称行列となっています。これより、プログラム中でヘッセ行列の要素を計算する際は、下三角成分か上三角成分のいずれかのみを計算すれば十分です。

---

例えば、2次のヘッセ行列は$$\boldsymbol{H}=\left[\begin{array}{cc}
\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial x^2} & \dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial x \partial y} \\
\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial y \partial x} & \dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial y^2}
\end{array} \right]$$となるので、逆行列は$${\boldsymbol{H}}^{-1}=\dfrac{1}{\det(\boldsymbol{H})}\left[\begin{array}{cc}
\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial y^2} & -\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial x \partial y} \\
-\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial y \partial x} & \dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial x^2}
\end{array} \right]$$で与えられます。ここで$\det(\boldsymbol{H})$はヘシアン$\boldsymbol{H}$の行列式で、この場合$$\det(\boldsymbol{H})=\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial x^2}\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial y^2}-\left(\dfrac{{\partial}^2 f(\boldsymbol{x})}{\partial x \partial y}\right)^2$$で定義されます。$\det(\boldsymbol{H})=0$ のときは逆行列が定義できません。

 

 ニュートン法のコンセプト

ニュートン法は目的関数の2次近似によって最適化を行う手法です。ここで関数の近似について復習しておきましょう。

テイラー展開を思い出して下さい。これは大学教養レベルの解析学で必ず習うものですが、微分可能な関数を多項式で近似するための方法です。知識が怪しい方は「『テイラー展開』の分かりやすい解説」のページで確認して下さい。

---

まず1次元の場合を考えます。$f\left(x+\Delta x\right)$を2次の項までテイラー展開すると、$$f\left(x+\Delta x\right)=f\left(x\right)+f^{\prime}\left(x\right) \Delta x+\frac{1}{2} f^{\prime \prime}\left(x\right) \Delta x^{2}+o\left(\Delta x^{3}\right)$$となります。$f\left(x+\Delta x\right)$の極値ではこれを$\Delta x$で微分した値が$0$となるので、$$f^{\prime}\left(x\right)+f^{\prime \prime}\left(x\right) \Delta x=0$$の解の方向、つまり $\small \Delta x=-\dfrac{f^{\prime}\left(x\right)}{f^{\prime \prime}\left(x\right)}$ だけ座標更新すれば良いことが分かります。

以上の操作を$n$次元に一般化します。$\boldsymbol{\delta}$を列ベクトルとして、座標更新の際に$\boldsymbol{\delta}$だけ点を移動するとします。$f\left(\boldsymbol{x}+\boldsymbol{\delta}\right)$を2次の項までテイラー展開すると、$$f\left(\boldsymbol{x}+\boldsymbol{\delta}\right)=f\left(\boldsymbol{x}\right)+\nabla f\left(\boldsymbol{x}\right)^{\top} \boldsymbol{\delta}+\frac{1}{2} \boldsymbol{\delta}^{\top} \nabla^{2} f\left(\boldsymbol{x}\right) \boldsymbol{\delta}+o\left(\|\boldsymbol{\delta}\|^{2}\right)$$となります。微小項を無視して両辺をベクトル$\boldsymbol{\delta}$で微分すれば（細かい計算は省きますが）$$\dfrac{\mathrm{d}}{\mathrm{d}\boldsymbol{\delta}}f\left(\boldsymbol{x}+\boldsymbol{\delta}\right) \approx 0+\nabla f\left(\boldsymbol{x}\right)+\nabla^{2} f\left(\boldsymbol{x}\right) \boldsymbol{\delta}+0$$となり、1次元の場合と同様に$\dfrac{\mathrm{d}}{\mathrm{d}\boldsymbol{\delta}}f\left(\boldsymbol{x}+\boldsymbol{\delta}\right)$が$0$になる座標に移動すれば良いので、$$\nabla f\left(\boldsymbol{x}\right)+\nabla^{2} f\left(\boldsymbol{x}\right) \boldsymbol{\delta}=0$$ $$\therefore \boldsymbol{\delta}=-(\nabla^{2} f\left(\boldsymbol{x}\right))^{-1}\nabla f\left(\boldsymbol{x}\right)$$を得ます。

ここで、ヘッセ行列$\nabla^{2} f\left(\boldsymbol{x}\right)$が逆行列を持つことを暗に仮定しています。もしヘッセ行列が逆行列を持たなければニュートン法は破綻します。また、ヘッセ行列が正定値行列でなければ「ニュートン方向」は目的関数の降下方向になりません。

※ヘッセ行列の逆行列が上手く計算できないと、ニュートン法は目的関数の値が増加する点を与えたり、計算効率が低下したりします。

※ニュートン法によって更新される新しい座標の方向（ベクトル）を「ニュートン方向」と呼びます。

---

ニュートン法のアルゴリズムは次のようになります。

ステップ幅$\alpha$はニュートン法の理論には登場しませんが、関数によっては正しい極小点に辿り着かない場合も多々あるので用意してあります。なお、幾何学的に解釈すると、更新先の座標$\boldsymbol{x}^{(i+1)}$は、点$\boldsymbol{x}^{(i)}$において目的関数の等高線の接線と曲率を共有するように目的関数を近似した多次元楕円体の中心に相当しています。ニュートン法ではこのようにして極小点を探します。

 

 ニュートン法は「2次収束」する

ニュートン法と最急降下法を比較してみます。$$f(x_1,x_2) = 2\cos(2^{x_1}-x_2^2+1)+e^{(x_1^2+x_2^2)/6}$$という式で表される2変数関数を例に考えてみましょう。初期値$(1.1,\, 0.5)$からステップ幅$1.0$で最適化してみます。

最急降下法だと57回のステップ数で極小点に到達しました。

これに対して、ニュートン法ではステップ数はたったの6回で済みました。これは最急降下法に比べてかなり速い収束速度です。

ニュートン法の収束速度が速いのは偶然ではありません。数式を使って理論的にきちんと説明するのは難しいのですが、最急降下法は「1次収束」するのに対してニュートン法では「2次収束」することが知られています。

点 $\boldsymbol{X}^{*}$に収束する点列 $\{\boldsymbol{x}_{n}\}$ に対して、ある定数$c$（$0<c<1$）と $k \in \mathbb{Z}^{+}$が存在して、任意の $n \geqq k$ に対して$$\|\boldsymbol{x}_{n+1}-\boldsymbol{X}^{*}\| \leqq c\|\boldsymbol{x}_{n}-\boldsymbol{X}^{*}\|$$が成り立つとき、点列 $\{\boldsymbol{x}_{n}\}$ は点 $\boldsymbol{X}^{*}$に1次収束するといいます。例えば最急降下法は1次収束する最適化手法です。

また、点 $\boldsymbol{X}^{*}$に収束する点列 $\{\boldsymbol{x}_{n}\}$ に対して、ある定数$c$（$0<c<1$）と $k \in \mathbb{Z}^{+}$が存在して、任意の $n \geqq k$ に対して$$\|\boldsymbol{x}_{n+1}-\boldsymbol{X}^{*}\| \leqq c\|\boldsymbol{x}_{n}-\boldsymbol{X}^{*}\|^{\color{red}{2}}$$が成り立つとき、点列 $\{\boldsymbol{x}_{n}\}$ は点 $\boldsymbol{X}^{*}$に2次収束するといいます。

ニュートン法は2次収束しますが、これは目的関数を2次近似していることによって近似の精度が上がっているためです。テイラー展開の次数が高ければ高いほど近似の精度が上がるのはよく知られている事実ですね。ただし、ニュートン法は計算すべき項の数が増えるため計算量が膨大になり、大抵は実際の問題に利用できるほどのパフォーマンスにはなりません。次元が増えるほどニュートン法の計算コストが大きくなるということは覚えておきましょう。

※ニュートン法の2次収束性の詳しい証明については「非線形計画法 (応用最適化シリーズ)」や「工学基礎 最適化とその応用 (新・工科系の数学)」などの書籍が参考になります。

---

以下にサンプルコードを置いておくので、適宜参考にしてみてください。2次偏導関数について、$$\dfrac{\partial^{2} f}{\partial x \partial y}=\dfrac{\partial^{2} f}{\partial y \partial x}$$が成り立つので、2次偏導関数「fxy」は1つだけしか用意していません。

import math

def f(x, y):
    return 2*math.cos(2**x+1-y**2)+math.e**((x**2+y**2)/6)

def fx(x, y):
    h = 1e-7
    return (f(x+h, y)-f(x-h, y))/(2*h)

def fy(x, y):
    h = 1e-7
    return (f(x, y+h)-f(x, y-h))/(2*h)

def fxx(x, y):
    h = 1e-7
    return (fx(x+h, y)-fx(x-h, y))/(2*h)

def fxy(x, y):
    h = 1e-7
    return (fx(x, y+h)-fx(x, y-h))/(2*h)

def fyy(x, y):
    h = 1e-7
    return (fy(x, y+h)-fy(x, y-h))/(2*h)

# input section
xc1 = 1.1; yc1 = 0.5     # 初期値
stepsize = 1.0           # ステップ幅
maxtimes = 10000         # 最大試行回数
times = maxtimes         # 試行回数カウント用変数

# Newton-Raphson method ニュートン・ラフソン法
for i in range(1, maxtimes):
    diffx1 = fx(xc1, yc1); diffy1 = fy(xc1, yc1)
    if math.sqrt(diffx1 ** 2 + diffy1 ** 2) < 1e-7:
        times = i
        break
    else:
        det = fxx(xc1, yc1) * fyy(xc1, yc1) - fxy(xc1, yc1) ** 2
        x_element = fx(xc1, yc1) * fyy(xc1, yc1) - fy(xc1, yc1) * fxy(xc1, yc1)
        y_element = - fx(xc1, yc1) * fxy(xc1, yc1) + fxx(xc1, yc1) * fy(xc1, yc1)
        # ステップ幅はスケーリングに利用（オプション）
        xc2 = xc1 - stepsize * x_element / det   # 次のx座標を生成
        yc2 = yc1 - stepsize * y_element / det   # 次のy座標を生成
        xc1 = xc2; yc1 = yc2   # 座標の更新

# 試行回数
print("done! ( i =", times, ")")
# 極小点の座標
print("(", xc1, ",", yc1,")")
# 極小点での値と勾配の大きさ
print("z =", f(xc1, yc1), ", |d| =", math.sqrt(diffx1 ** 2 + diffy1 ** 2))

本来、理論的にはニュートン法にステップ幅というものは登場しませんが、初期値が悪いとあらぬ方向に座標が更新されてしまうことがしばしばあります。これを防止するために一定のステップ幅（縮小因子）を考慮できるようにしてあります。

行列計算に慣れていないと、特にヘッセ行列の扱い方は難しく感じられると思います。上のプログラムは2次元の場合にしか対応していませんが、実際には多次元の空間でニュートン法を用いることもあります。その際は多くの場合、numpy等に付属しているモジュールを使用するのが便利です。

 

 まとめ

ニュートン法は収束が速い最適化手法ですが、時に目的関数の値が増加する方向に向かってしまうことがあります。これはニュートン法が停留点の方向に向かって座標を更新する性質を持っているためです。これを利用すれば、極大・極小の他に「鞍点」と呼ばれる停留点を見つけることが可能になります。次回はニュートン法の応用として鞍点の探索について考えてみます。