Python
Haskell
数学
最適化
python3

Haskell (と Python3)で最急降下法

初投稿です.
先日最適化理論の講義のレポートで最急降下法のプログラムをPythonで書いたので,
Haskellでも書いてみました.

HaskellはPythonに比べてまだまだ書き方を理解していないので,
こうすればいいよ!というのがあればご教示ください.

最急降下法

次のような制約なし最適化を考えます.

\min_{\mathbf{x} \in \mathbb{R}^n} \; f( \mathbf{x} )

つまり,やりたいことは,$ f(\mathbf{x}) $の最小値と最小解をもとめることです
($f$は凸関数と仮定しています).
例として,$ f(x, y) = (x-1)^2 + (y-2)^2 $を考えてみます.
この関数は,$ (x, y) = (1, 2) $で,最小値$ 0 $を取ります.
図で書くと,次のようになります.
等高線.png

この関数の$(a, b)$における勾配$ \nabla f(a, b) $は,図で書くと次のようになります(緑色).
等高線2.png
最急降下法とは,
「$(a, b)$から$ - \nabla f(a, b) $の方にある程度進む」
という戦略をとる最適化の手法です.
この「ある程度」というのは,適当ではありません.
関数値が下がる方向に進みすぎても,進まなすぎてもダメです(下の図を参考程度に).
Picture1.png

この,「どれくらい進むか」という進み具合を表す幅を「ステップ幅」といい,あとでプログラムを書くときに説明します.

最急降下法の概形を書くと,次のようになります.
1. 初期点$x = (x_0, y_0)$, 誤差のパラメータ$eps$を入力として与える.
2. $\nabla f(x) < eps$ならば,5へ.
3. ステップ幅 $t > 0$を決める.
4. $ x \gets x + t (-\nabla f(x)) $として,2へ.
5. return $x$

最適解 $ ( x_{\text{opt}}, y_{\text{opt}} ) $ では勾配が $\mathbf{0}$ になるのでそれを終了条件とすればただしく動くはずです.

モジュールのインポート

今回は,Haskellでは標準モジュールのみで解くので,何もインポートしません.

main.hs
-- 何もいらない(便利なモジュールを知らない)

Pythonでは,numpyのみを使います.

main.py
import numpy as np

型宣言

つぎに,ベクトルを表現するのに[Double]などを用いるとわかりにくくなるので,名前を変えます.

main.hs
type Scalar = Double -- Double型をScalarと呼ぶことにする
type Vector = [Double] -- Double型のリストをVectorと呼ぶことにする

目的関数の定義

今回は,$f(x, y) = (x-1)^2 + (y-2)^2$を最小化することを考えます.
$f(x, y)$の$(x, y)$における微分は次のようになります.

\nabla f(x, y) = 
\begin{bmatrix}
    2(x-1) \\
    2(y-2)
\end{bmatrix}

この$f, \nabla f$をHaskell, Pythonでそれぞれ書きます.

main.hs
f :: Vector -> Scalar
f (x:y:_) = (x-1.0)^2 + (y-1.0)^2

f' :: Vector -> Vector
f' (x:y:_) = [a, b]
    where a = 2*(x - 1.0)
          b = 2*(y - 2.0)
main.py
def f(x):
    return (x[0] - 1.0)**2 + (x[1] - 1.0)**2

def df(x):
    a = 2*(x[0] - 1.0)
    b = 2*(x[1] - 2.0)
    return np.array([a, b])

基本的な関数の定義

僕はHaskellにおけるpythonのnumpyのようなモジュールを知らないので,ベクトルの内積,和,スカラー倍等の計算をする関数を定義します(便利なモジュールがあれば教えてください).

main.hs
-- ベクトルのスカラー倍
mul :: Scalar -> Vector -> Vector
mul c v = map (* c) v

-- ベクトルの和
add :: Vector -> Vector -> Vector
add xs ys = zipWith (+) xs ys

-- ベクトルの内積
dot :: Vector -> Vector -> Scalar
dot xs ys = foldr (+) 0 $ zipWith (*) xs ys

-- ベクトル v を負にする i.e. (-1) * v を計算する
negative :: Vector -> Vector
negative = map (* (-1))

-- ベクトルの2-ノルムを計算する
norm :: Vector -> Scalar
norm xs = sqrt $ dot xs xs

-- Double型の最大値(maxBoundがなかった?)
doubleMax :: Double
doubleMax = 1e10

pythonでは,numpyがあるのでベクトルの内積,和,スカラー倍等は定義しなくても使えます.
2-ノルムを計算する関数だけ作ります.

main.py
def norm(x):
    n = np.sum(x**2)
    return np.sqrt(n)

 main関数の完成像

今回の記事で書くプログラムの完成像は,Haskellだと次のような感じになります.

main.hs
main :: IO ()
main = do
    let x  = [5.0, 5.0] -- 初期点
        c1 = 0.4        -- パラメータ(後で説明)
        c2 = 0.8        -- パラメータ(後で説明)
        (v_opt, x_opt) = solve x c1 c2 -- 最適解を求める.
    putStr "Optimal Solution : "
    print x_opt
    putStr "Optimal Value    : "
    print v_opt

Pythonだと,次のような感じを目標にプログラムを書いていきます.

main.py
if __name__ == "__main__":
    x0  = np.array([5.0, 5.0]) # 初期点
    c1  = 0.4                  # パラメータ(後で説明)
    c2  = 0.8                  # パラメータ(後で説明)
    eps = 1e-9                 # 終了条件用のパラメータ
    sd = SteepestDescent(x0=x0, c1=c1, c2=c2, eps=eps)
    sd.solve()
    print("Optimal Solution : {}".format(sd.opt))
    print("Optimal Value    : {}".format(f(sd.opt)))

じゃあ,後は最適解を計算する関数(Pythonだとメソッド)solveを作れば終わりです.

その前にPythonでは

クラスを使ってプログラムを書いていきます.
SteepestDescentクラスのコンストラクタを次のように作ります.

main.py
class SteepestDescent(object):
    """
    STEEPEST DESCENT METHOD
    """
    def __init__(self, x0, c1=0.4, c2=0.8, eps=1e-6):
        """
        PARAMETERS:
            x    : 時刻tに選んだ点
            opt  : 最適解を格納する変数
            c1   : パラメータ(armijoの条件)
            c2   : パラメータ(wolfの条件)
            eps  : 誤差パラメータ
        """
        self.x    = x0
        self.opt  = None
        self.c1   = c1
        self.c2   = c2
        self.eps  = eps

終了条件

終了条件は,$\nabla f(x, y) < eps$なので,Haskellでは次のように書きます.

main.hs
eps :: Double
eps = 1e-9

isFiniteState :: Vector -> Bool
isFiniteState v =
    let n = (norm . f') v
        in if n < eps then True else False

Pythonでは,次のように書きます.

main.py
def isFinite(self):
    g = df(self.x)
    n = norm(g)
    if n < self.eps:
        return True
    else:
        return False

次に,solve関数(メソッド)の概形を書きます.

main.hs
solve :: Vector -> Scalar -> Scalar -> (Scalar, Vector)
solve x c1 c2
    | isFiniteState x = (f x, x)                  -- 終了条件を満たしたら終了
    | otherwise       =                           -- 満たしていなかったら,
        let d = negative (f' x)                   -- -f'(x) を計算して,
            t = getStepSize x 1 c1 c2 0 doubleMax -- ステップ幅を求めて,
            y = add x (mul t d)                   -- x = x - t f'(x)として
            in solve y c1 c2                      -- 解き直す

Pythonだと,こう書きました.

main.py
    def solve(self):
        while not self.isFinite():   # 終了条件を満たしていない間
            d = -df(self.x)          # - f'(x)を計算して,
            t = self.getStepSize()   # ステップ幅を求めて,
            self.x = self.x + t * d  # x = x - t f'(x)として解き直す
        self.opt = self.x

Haskell, Pythonのどちらにおいても,あとはステップ幅を求めるだけです.

ステップ幅をもとめる

最初の方で言ったように,ステップ幅は大きすぎても小さすぎてもいけません(再掲).
Picture1.png

理想としては,$ \mathbf{d} = - \nabla f(\mathbf{x}) $としたときに

\min_t  f( \mathbf{x} + t \mathbf{d})

を満たすtが最も良いステップ幅です.これをexact line searchと言うそうです.
Picture2.png
しかし,exact line searchをしようとすると別の最適化問題を解かなければならないので,
今回はinexact line searchでステップ幅を求めようと思います.
inexact line search とは,大きすぎず小さすぎないステップ幅$ t > 0 $をもとめる方法です.
つまり,次のようなプログラムを書きたいわけです.

  1. ステップ幅$t>0$を渡す
  2. if $t$が大きすぎる$ \to $ $t$を小さくして,1へ
  3. else if $t$が小さすぎる$ \to $ $t$を大きくして,1へ
  4. else return t

ステップ幅が大きすぎないか,という問題はArmijoの条件というものを使って判定できます.
Armijoの条件:
$ \phi(t) = f( \mathbf{x} + t \mathbf{d})$としたときに,$c_1 \in (0, 1)$に対して

\phi (t) \leq \phi(0) + c_1 t \phi'(0)

を満たす$t > 0$を選ぶ.

また,ステップ幅が小さすぎないか,という問題は,Wolfの条件というものを使って判定できます.
Wolfの条件:
$0 < c_1 < c_2 < 1$に対して,

\phi (t)  \leq \phi(0) + c_1 t \phi'(0)\\
\phi' (t) \geq c_2 \phi' (0)

を満たす$t > 0$を選ぶ.
(Wolfの条件の1つ目は,Armijoの条件のことです)

この辺は数学書に書いてあると思うので,なんでこれを満たせば良いか,等が気になる人は調べてみてください.
どういう式を満たせばステップ幅が適当か,という式がわかったので,あとはプログラムを書くだけです.

追記:コメントアウトしてあるところは僕が書いたコードで,
コメントアウトされていないところは,@kmrt7010さんにおしえていただいたコードです!
ありがとうございます!

main.hs
-- コメントアウトされているところは
-- 僕が書いた綺麗ではないコードです

{-
armijo :: Vector -> Scalar -> Scalar -> Bool
armijo v t c =
    let z  = f v
        d  = negative (f' v)
        dx = add v (mul t d)
        y  = f dx
        p0 = dot (f' v) d
        in if y <= z + t*c*p0 then True
           else False
-}

armijo :: Vector -> Scalar -> Scalar -> Bool
armijo v t c = y <= z + t*c*p0
  where z  = f v
        d  = negative (f' v)
        dx = add v (mul t d)
        y  = f dx
        p0 = dot (f' v) d

{-
wolf :: Vector -> Scalar -> Scalar -> Bool
wolf v t c =
    let d  = negative (f' v)
        dx = add v (mul t d)
        p0 = dot (f' v) d
        pt = dot (f' dx) d
        in if pt >= c*p0 then True
           else False
-}

wolf :: Vector -> Scalar -> Scalar -> Bool
wolf v t c = pt >= c*p0
    where d  = negative(f' v)
          dx = add v (mul t d)
          p0 = dot (f' v) d
          pt = dot (f' dx) d

{-
getStepSize :: Vector -> Scalar -> Scalar -> Scalar -> Scalar -> Scalar -> Scalar
getStepSize x t c1 c2 a b
    | not (armijo x t c1) = if t < doubleMax then getStepSize x ((a + t)/2.0) c1 c2 a t
                            else getStepSize x (2*a) c1 c2 a t
    | not (wolf x t c2) = if b < doubleMax then getStepSize x ((t + b)/2.0) c1 c2 t b
                          else getStepSize x (2*t) c1 c2 t b
    | otherwise = t
-}

getStepSize :: Vector -> Scalar -> Scalar -> Scalar -> Scalar -> Scalar -> Scalar
getStepSize x c1 c2 a t b
   | not (armijo x t c1) && t < doubleMax = getStepSize x c1 c2 a ((a+t)/2) t
   | not (armijo x t c1)                  = getStepSize x c1 c2 a (2*a)     t
   | not (wolf x t c2)   && t < doubleMax = getStepSize x c1 c2 t ((t+b)/2) b
   | not (wolf x t c2)                    = getStepSize x c1 c2 t (2*t)     b
   | otherwise                            = t

この辺の書き方が少し汚い気がするので,良い書き方があれば教えて欲しいです.
一方,Pythonでは次のようにかけます.

main.py
    def armijo(self, t, c):
        z  = f(self.x)
        d  = -df(self.x)
        dx = self.x + t*d
        y  = f(dx) # left hand side of armijo's condition
        p0 = np.dot(df(self.x), d) # armijo's cond. phi'(t)
        if y <= z + t*c*p0:
            return True
        else:
            return False

    def wolf(self, t, c):
        d  = -df(self.x)
        dx = self.x + t*d
        p0 = np.dot(df(self.x), d) # phi'(0)
        pt = np.dot(df(dx), d) # phi'(t)
        if pt >= c*p0:
            return True
        else:
            return False

    def getStepSize(self):
        # initialize parameters
        alpha = 0
        beta  = sys.maxsize
        t     = 1
        c1    = self.c1
        c2    = self.c2
        while True:
            if not self.armijo(t, c1):
                beta = t
            elif not self.wolf(t, c2):
                alpha = t
            else:
                break
            if beta < sys.maxsize:
                t = (alpha + beta) / 2.0
            else:
                t = 2*alpha
        return t

はい,これで欲しい関数などは全部求められました.
あとは実行するだけです.

実行結果

$f(x, y) = (x-1)^2 +(y-2)^2$に対して実行してみました.
Haskellでの実行結果のみ載せてます.

result.txt
Optimal Solution : [1.0, 2.0]
Optimal Value    : 0.0

正しく計算できてますね!

まとめとか

ほんとはNewton法を書きたかったけど,Haskellで固有値を求める関数がどれかわからなかったので断念しました(Pythonだとnp.linalg.eigvals(A)で行列Aの固有値を求められた).
こういうのがあるよ!っていうのを知ってる方がいらっしゃいましたら,是非教えてください.
この前@lotzさんにNumeric.LinearAlgebraに固有値計算の関数があることを教えていただきました.調べてみると,コレスキー分解もできるみたいなので,Newton法も書けそうです!
近いうちに,研究や心に余裕ができたら書かせていただきます.

@lotzさん,@kmrt7010さん,アドバイスありがとうございます!