In [1]:

import numpy as np
import matplotlib.pyplot as plt

講義概要

凸最適化問題
最急降下法
ニュートン法
準ニュートン法

凸最適化問題

凸関数¶

任意の $x,x^{\prime}\in\mathcal{X}$、任意の $\theta\in[0,1]$ に対して、以下が成り立つような関数を凸関数(convex function)と呼ぶ。 $$f\left(\left(1-\theta\right)x + \theta x^{\prime}\right)\leq\left(1-\theta\right)f(x) + \theta f\left(x^{\prime}\right)$$
$f$ が一階微分可能の時、$f$ が凸関数になるための必要十分条件は、感覚的には接戦が常に関数の下に来ることで、 $$\forall x,x^{\prime}\in\mathcal{X},\quad f\left(x^{\prime}\right)\geq f(x) + \nabla f(x)^T\left(x^{\prime}-x\right)$$
$f$ が二階微分可能の時、$f$ が凸関数になるための必要十分条件は、ヘッセ行列(Hessian matrix)$\nabla^2f(x)$ が半正定値(positive semi-definite)↔︎$\nabla^2f(x)$ の固有値(eigenvalue)が全て非負
任意の $x,x^{\prime}\in\mathcal{X}$、任意の $\theta\in[0,1]$ に対して、$(1-\theta)x+\theta x^{\prime}$ も $\mathcal{X}$ に属するとき、$\mathcal{X}$ を凸集合(convex set)と呼ぶ。$x$ と $x^{\prime}$ を結ぶ線分上の全ての点が $\mathcal{X}$ に属するということ。

最適化¶

最適化問題(optimization problem)： $\mathcal{X}\subset\mathbb{R}^d$ 上に定義される関数 $f:\mathcal{X}\rightarrow\mathbb{R}$ の最小値 $\underset{\mathbf{x}\in\mathcal{X}}{\min}f(\mathbf{x})$ を求めよ。
最適化アルゴリズム(optimization algorithm)：
- 収束する点列の生成アルゴリズム
- 最適性条件：最適化問題の最適解であるための必要条件
制約なし最適化問題(unconstrained optimization problem)：
- $\mathcal{X}\subset\mathbb{R}^d$ 上に定義される関数 $f:\mathcal{X}\rightarrow\mathbb{R}$ の最小値 $\underset{\mathbf{x}\in\mathcal{X}}{\min}f(\mathbf{x})$ を求めよ。
- $f$ を目的関数(objective function)と呼ぶ。
- 目的関数は微分可能(differentiable)と仮定する。
- 最適性条件：$\mathbf{x}^{\star}$ が最適解のとき $$\nabla f(\mathbf{x}^{\star}) = 0$$
凸最適化問題(convex optimization problem)：
- 関数 $f$ が凸で集合 $\mathcal{X}$ も凸
- 最適解（の集合）が一意に定まる。
- 最適解の必要十分条件：勾配(gradient)がゼロ（凸でない場合は必要条件）

最急降下法

最急降下法(steepest descent method)¶

※ 勾配法(gradient method)とも。

適当に初期値 $\mathbf{x}_0$ を定める。
$k=0,1,2,\ldots$ に対して、収束するまで以下を繰り返す。 $$\mathbf{x}_{k+1} = \mathbf{x}_k - \varepsilon_k\nabla f\left(\mathbf{x}_k\right)\quad (\varepsilon_k>0\text{ : step size})$$

凸最適化問題の場合、大域的最適解(global optimization solution)が求まる
非凸最適化問題の場合、局所最適解(local optimization solution)しか求められず、様々な初期値から何回か学習して最適解を採用する、ことしかできない。

ステップ幅(step size)の選択¶

焼きなまし(annealing)：「最初は大きく、徐々に小さく」。しかし、適切に実装することは容易でない。
正規化(normalization)：$\varepsilon_k = \varepsilon_k^{\prime}/\|\nabla f(\mathbf{x}_k)\|$ 勾配が大きいとき安定するが、勾配が小さいときは不安定

厳密直線探索(exact line search)¶

目的関数の値を最小にするステップ幅を求める。 $$\underset{\varepsilon_k>0}{\min}f\left(\mathbf{x}_k-\varepsilon_k\nabla f\left(\mathbf{x}_k\right)\right)$$
$f$ が単峰関数(unimodal function)の場合は、二分探索法(binary search)や黄金分割探索(golden section search)などの探索によって最小値を探す。
一般に、一次元の非線形最小化問題を解く必要があるため、最適なステップ幅の探索に時間がかかる。
厳密直線探索を用いた最急降下法は、一次収束(linear convergence)する。 $$\text{Q-linear convergence: }\|\mathbf{x}_{k+1}-\mathbf{x}^{\star}\| c\leq \|\mathbf{x}_k-\mathbf{x}^{\star}\|^q,\quad 0<c<1$$

In [2]:

ITERATION = 10
N = 100
fsize = 4

In [3]:

a = 10; b = 1
A = 2*np.diag([a,b])
func = lambda x,y: a*x**2+b*y**2

In [4]:

x1min=-5; x1max=5; x2min=-5; x2max=5
initialX = np.array([x1min+1,x2max-2])

x1 = np.linspace(x1min, x1max, N)
x2 = np.linspace(x2min, x2max, N)
X1,X2 = np.meshgrid(x1,x2)
Z = func(X1,X2)

In [5]:

def find_step_size(method, f, X, A, **kwargs):
    df = A.dot(X)
    if method=="fixed":
        epsilon = 1e-1
    elif method=="exact line search":
        epsilon = np.sum(np.square(df)) / np.sum(np.square(df)*A.diagonal())
    elif method=="backtracking line search":
        epsilon = 1
        alpha = kwargs["alpha"]; beta = kwargs["beta"]
        while True:
            newX = X-epsilon*df
            if f(newX[0],newX[1])-f(X[0],X[1]) <= -alpha*epsilon*np.sum(np.square(df)): break
            epsilon = beta*epsilon
    return epsilon

In [6]:

methods = ["fixed", "exact line search"]
n_methods = len(methods)

In [7]:

fig = plt.figure(figsize=(fsize*n_methods,fsize*2))
for i,method in enumerate(methods):
    X = initialX; pX = initialX
    objective_functions = []
    
    ax_t = fig.add_subplot(2,n_methods,i+1)
    ax_t.contour(X1, X2, Z, 30)
    ax_t.scatter(X[0],X[1],color="red",s=100,marker="*",label="initial")
    
    objective_functions.append(func(X[0],X[1]))
    for it in range(ITERATION):
        epsilon = find_step_size(method, func, X, A, alpha=0.5, beta=0.8)
        df = A.dot(X)
        X = X-epsilon*df
        ax_t.plot((X[0],pX[0]), (X[1],pX[1]), color="black", linestyle=":")
        ax_t.scatter(X[0],X[1],color="black",s=100,marker="*")
        objective_functions.append(func(X[0],X[1]))
        pX=X    
    ax_t.set_title(method, fontsize=18), ax_t.set_xlabel("x1"), ax_t.set_ylabel("x2")
    ax_t.legend()
    
    sharey = ax_b if i!=0 else None
    ax_b = fig.add_subplot(2,n_methods,n_methods+i+1, sharey=sharey)
    ax_b.plot(objective_functions, color="black")
    ax_b.set_xlabel("iteration"), ax_b.set_ylabel("objective function")
plt.tight_layout()
plt.show()

例：$f(x) = 1/2\mathbf{x}^T\mathbf{Ax}$¶

$$ \begin{aligned} \frac{\partial f\left(\mathbf{x}_k-\varepsilon_k\nabla f\left(\mathbf{x}_k\right)\right)}{\partial\varepsilon_k} &= \frac{\partial}{\partial\varepsilon_k}\frac{1}{2}\left(\mathbf{x}_k-\varepsilon_k\nabla f\left(\mathbf{x}_k\right)\right)^T\mathbf{A}\left(\mathbf{x}_k - \varepsilon_k\nabla f\left(\mathbf{x}_k\right)\right)\\ &= \varepsilon_k\nabla f\left(\mathbf{x}_k\right)^T\mathbf{A}\nabla f\left(\mathbf{x}_k\right) - \nabla f\left(\mathbf{x}_k\right)^T\mathbf{Ax}_k\\ &=0\\ \therefore \varepsilon_k^{\star} &= \frac{\|\nabla f\left(\mathbf{x}_k\right)\|^2}{\nabla f\left(\mathbf{x}_k\right)^T\mathbf{A}\nabla f\left(\mathbf{x}_k\right)}\qquad (\nabla f\left(\mathbf{x}_k\right) = \mathbf{Ax}) \end{aligned} $$

バックトラック直線探索(backtracking line search)¶

$\varepsilon_k=1$ に初期化
アルミホ基準(Armijo rule)： $$f\left(\mathbf{x}_k-\varepsilon_k\nabla f\left(\mathbf{x}_k\right) - f\left(\mathbf{x}_k\right)\right) \leq\alpha\|\nabla f\left(\mathbf{x}_k\right)\|^2$$ が成り立つまで $\varepsilon_k$ を $\varepsilon_k\leftarrow\beta\varepsilon_k$ と減衰させる。（$0<\alpha, \beta<1$）
目的関数の減少量の線形予測値の $\alpha$ 倍の減少を保証している。
- $f$ の真の減少量： $$g\left(\varepsilon_k\right) = f\left(\mathbf{x}_k-\varepsilon_k\nabla f\left(\mathbf{x}_k\right)\right) - f\left(\mathbf{x}_k\right)$$
- $f$ の減少量の線形予測値 $$\begin{aligned}\varepsilon_kg^{\prime}(0) &=\varepsilon_k\cdot\left( -\nabla f\left(\mathbf{x}_k\right)^T\nabla f\left(\mathbf{x}_k-0\cdot\nabla f\left(\mathbf{x}_k\right)\right)\right)\\&=-\varepsilon_k\nabla f\left(\mathbf{x}_k\right)^T\nabla f\left(\mathbf{x}_k\right)\\&=-\varepsilon_k\|\nabla f\left(\mathbf{x}_k\right)\|^2 \end{aligned}$$

In [8]:

methods = ["fixed", "exact line search"]
n_methods = len(methods)

In [9]:

fig = plt.figure(figsize=(fsize*n_methods,fsize*2))
for i,method in enumerate(methods):
    X = initialX; pX = initialX
    objective_functions = []
    
    ax_t = fig.add_subplot(2,n_methods,i+1)
    ax_t.contour(X1, X2, Z, 30)
    ax_t.scatter(X[0],X[1],color="red",s=100,marker="*",label="initial")
    
    objective_functions.append(func(X[0],X[1]))
    for it in range(ITERATION):
        epsilon = find_step_size(method, func, X, A, alpha=0.5, beta=0.8)
        df = A.dot(X)
        X = X-epsilon*df
        ax_t.plot((X[0],pX[0]), (X[1],pX[1]), color="black", linestyle=":")
        ax_t.scatter(X[0],X[1],color="black",s=100,marker="*")
        objective_functions.append(func(X[0],X[1]))
        pX=X    
    ax_t.set_title(method, fontsize=18), ax_t.set_xlabel("x1"), ax_t.set_ylabel("x2")
    ax_t.legend()
    
    sharey = ax_b if i!=0 else None
    ax_b = fig.add_subplot(2,n_methods,n_methods+i+1, sharey=sharey)
    ax_b.plot(objective_functions, color="black")
    ax_b.set_xlabel("iteration"), ax_b.set_ylabel("objective function")
plt.tight_layout()
plt.show()

微分不可能な場合¶

目的関数 $f$ が微分不可能(non-differentiable)の場合は、微分の概念を一般化した劣微分を用いる。
凸関数 $f$ の $\mathbf{x}^{\prime}$ での劣勾配(sub-gradient)とは、全ての $\mathbf{x}\in\mathcal{X}$ に対して次式を満たす $\mathbf{\xi}$ $$f(\mathbf{x})\geq f\left(\mathbf{x}^{\prime}\right) + \mathbf{\xi}^T\left(\mathbf{x}-\mathbf{x}^{\prime}\right)$$
- $f$ が微分可能な時、$\mathbf{\xi} = \nabla f\left(\mathbf{x}^{\prime}\right)$
- 上式を満たす $\mathbf{\xi}$ 全体を $\partial f\left(\mathbf{x}^{\prime}\right)$ で表し、劣微分(sub-differential)と呼ぶ。
微分不可能な点では、劣微分のどれかの値を用いる。

ニュートン法

一般にニュートン法という名称は、方程式の解を求めるアルゴリズムを指すが、以下を表すことが多い。

適当に初期値 $\mathbf{x}_0$ を定める。
$k=0,1,2,\ldots$ に対して、収束するまで以下を繰り返す。 $$\mathbf{x}_{k+1} = \mathbf{x}_k - \varepsilon_k\left(\nabla^2f\left(\mathbf{x}_k\right)\right)^{-1}\nabla f\left(\mathbf{x}_k\right)$$

勾配法では一階微分しか用いていないので、二階微分も用いることで反復回数を減らすことが狙い。実際 $k$ が十分大きい時、最適解 $\mathbf{x}^{\star}$ の近傍で二次収束(quadratic convergence)するが、逆行列を求める計算コストが大きい。
二次のテーラー展開(Taylor expansion)を用いて $f$ を現在の解 $\mathbf{x}_k$ の周りで近似 $\left(f(\mathbf{x})\approx f_k(\mathbf{x})\right)$ し、その二次近似を最小にする $\mathbf{x}_k$ に更新する。 $$f_k(\mathbf{x}) = f(\mathbf{x}_k) + \nabla f(\mathbf{x}_k)^T\left(\mathbf{x}-\mathbf{x}_k\right) + \frac{1}{2}\left(\mathbf{x}-\mathbf{x}_k\right)^T\nabla^2f(\mathbf{x}_k)(\mathbf{x}-\mathbf{x}_k)$$
$f_k$ の勾配をゼロと置いた方程式 $$\nabla f_k(\mathbf{x}) = \frac{\partial f_k(\mathbf{x})}{\partial \mathbf{x}}= \nabla f(\mathbf{x}_k) + \nabla^2f(\mathbf{x}_k)\left(\mathbf{x}-\mathbf{x}_k\right) = \mathbf{0}$$ を解けば、$f_k$ の最小解が以下のように求まる $$\mathbf{x}^{\star} = \mathbf{x}_k - \left(\nabla^2 f\left(\mathbf{x}_k\right)\right)^{-1}\nabla f\left(\mathbf{x}_k\right)$$
なお、ヘッセ行列が逆行列を持たない時は、単位行列の定数倍を加えて対処する（正則化, regularization）らしい。 $$\left(\nabla^2f\left(\mathbf{x}_k\right) + \mu\mathbf{I}\right)^{-1},\quad\mu>0$$

準ニュートン法

ニュートン法の近似。ニュートン法ではヘッセ行列の逆行列を求めるのに計算コストがかかるので、勾配ベクトル $\nabla f(\mathbf{x})$ を用いて近似計算する。
BFGSアルゴリズムを用いたヘッセ行列の更新式が有名

なお、レポート課題で詳しく記述しているが、BFGSアルゴリズムのヘッセ行列の更新式は以下で表される。

$$ \mathbf{H}_{k+1} =\mathbf{H}_k+\frac{\mathbf{t}_k\mathbf{t}_k^T}{\mathbf{s}_k^T\mathbf{t}_k}-\frac{\mathbf{H}_k\mathbf{s}_k\mathbf{s}_k^T\mathbf{H}_k}{\mathbf{s}_k^T\mathbf{H}_k\mathbf{s}_k} $$

なお、この時

$$ \begin{aligned} \begin{cases} \mathbf{s}_k &= \mathbf{x}_{k+1} - \mathbf{x}_k\\ \mathbf{t}_k &= \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_k) \end{cases} \end{aligned} $$

まとめ

凸最適化問題

目的関数が凸関数で定義域が凸集合
最適解が一意に定まる

最急降下法

勾配を降下するように値を更新
ステップ幅の選択が重要
目的関数が微分不可能な場合は、劣勾配法や近接勾配法を用いる

ニュートン法

二階微分の情報を利用する
ヘッセ行列の近似を用いる準ニュートン法が実用的

In [ ]:

連続最適化(1)制約なし最適化

講義概要

凸最適化問題

凸関数¶

最適化¶

最急降下法

最急降下法(steepest descent method)¶

ステップ幅(step size)の選択¶

厳密直線探索(exact line search)¶

例：$f(x) = 1/2\mathbf{x}^T\mathbf{Ax}$¶

バックトラック直線探索(backtracking line search)¶

微分不可能な場合¶

ニュートン法

準ニュートン法

まとめ

Table of Contents

Published

Last Updated

Category

Tags

Contact

Other contents