甲斐性なしのブログ

うどんくらいしか食べる気がしない

二乗誤差を評価値とする焼きなましの差分計算をEigenで行う(Atcoder Heuristic Contest 030)

はじめに

しばらくブログを書いてなかったのでAHC参加ネタを投稿。

先日Atcoder Heuristic Contest 030(AHC030)に参加し最終56位という結果だったが、提出した基本方針に二乗誤差を最小化する焼きなましが含まれていた。この方法では二乗誤差の計算を何度も行うことになるためその計算時間が肝になるのだが、線形代数ライブラリのEigenを使うと下手に差分計算するより愚直に全計算した方が速かったりする*1。とはいえ、少し工夫して差分計算にもEigenを活用すれば、状況によっては愚直計算よりも速くなる。

ということで、本記事では実際にAHC030で使ったEigenによる二乗誤差の差分計算の方法を説明する。このテクニックはAHC030に限らず二乗誤差を評価値として焼きなましする場合にも応用できる(かも?)。

AHC030の問題については以下を参照。 atcoder.jp

AHC030解法基本方針概要

解法基本方針についてはtwitterのリンクおよびそのスレッドを参照。

二乗誤差最小化によるポリオミノの配置最適化

まずは評価値の計算を行列演算の形で表す。占いをn回行ったとしてi回目の占いで得られる値をy_iy_ii番目の要素に持つn要素のベクトルを\mathbf{y}i行目にi回目の占いにおいてどのマスを選んだか?という情報を持つn \times N^2の行列を\mathbf{X}(各列がマスに対応し、選ばれたマスなら1, そうでなければ0)とする*2。そして、\mathbf{w}を現状態のポリオミノの配置において、各マス何個のポリオミノが重なっているかを表すN^2要素のベクトルとし、これが最適化対象の変数である。これらを使って最適化問題として定式化すると、

\displaystyle 
\begin{aligned}
&{\rm arg}\min_{\mathbf{p}} \parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2 \\
&{\rm s.t. \ \ }   \mathbf{w}はポリオミノの配置の組み合わせであり得る値
\end{aligned}
\tag{1}

となる。これを1つのポリオミノを選択して置く場所を変える焼きなましによって最適化を図る。そのため二乗誤差\parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2の計算を高速化できれば短い時間で焼きなましのループ回数が稼げ、より正解を出せる確率が上がる。

二乗誤差計算の高速化

Eigenの活用

C++においては行列計算を高速に実行できるライブラリEigenを活用すれば、愚直計算でもかなり速く計算できる。コードもこんな感じで1行で書ける。

double cost = (y - X * w).squaredNorm();

EigenではSIMD命令の活用、ループ変形等によりメモリアクセスの効率化などが行われているため演算回数だけ見れば遅くなりそうな処理も高速にできる場合がある。逆に差分計算はメモリアクセスが飛び飛びになりがちで、SIMD化もしにくいため処理時間上不利になることがある。

二乗誤差差分計算の行列表現

差分計算を行列演算で行うことを考える。ポリオミノの配置変更による\mathbf{w}の変化を\Delta \mathbf{w}とすると、二乗誤差の変化量Dは以下のようになる。

\displaystyle
\begin{align}
D &= \parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2 - \parallel \mathbf{y} - \mathbf{X} \left(\mathbf{w} + \Delta \mathbf{w} \right) \parallel_2^2 \\
&=2  \left(\mathbf{y}^T \mathbf{X} \Delta \mathbf{w} - \mathbf{w}^T \mathbf{X}^T \mathbf{X} \Delta \mathbf{w} \right) - \Delta \mathbf{w}^T \mathbf{X}^T \mathbf{X} \Delta \mathbf{w}
\end{align}
\tag{2}

この式の中で \mathbf{X}^T \mathbf{X} は焼きなましの中で変化しないので、事前に計算しておくことができる。 加えて、焼きなましの1ループでは1つのポリオミノしか動かさないので、\mathbf{w}の中で変化する要素はごく一部、つまり\Delta \mathbf{w}の大半の要素は0になる。このような多くの要素が0のベクトルを扱うためEigenにはSparseVectorが用意されている。\Delta \mathbf{w}をSparseVector型として扱うことで、式(2)の計算がだいぶ速くなる。コードはこんな感じ。

Eigen::SparseVector<double> sdw(N * N);
vector<Eigen::Triplet<double>> trivec;
rep(i, N * N){
    if(abs(dw(i)) > 1e-5){
        sdw.insert(i) = dw(i);
    }
}
VectorXd XXdw = XX * sdw;
VectorXd Xdw = X * sdw;
double D = 2.0 * (y.dot(Xdw) - old_w.dot(XXdw)) - dw.dot(XXdw);

比較

上記2方法で評価値計算をする焼きなましがatcoderのジャッジサーバー上で400[ms]内に何回回るかをnを変えて比較する。その際のテストデータはseed=1とする(N = 13, M = 6, \epsilon = 0.04)。結果は以下の通り。

n 全計算 差分計算
20 814301 342042
40 598892 338786
60 312957 320429
80 311225 319219
100 175905 303552
120 226957 297810
140 143947 275129
160 176061 250388
180 115353 249101
200 146547 229857

このように\mathbf{X}の行数(=占い回数)nが小さい間は愚直に全計算した方が回数を稼げるが、n \geq 100あたりで逆転し差分計算した方が速い*3。従って、実際提出したコードでは占いした回数nに応じて愚直全計算と差分計算を切り替える処理とした。

ちなみに

今回の手法を使ったseed=0のビジュアライザは以下の通り。

*1:焼きなましでは1回の更新の中で状態が一部しか変化しないことがしばしばある。そのため基本的にはその差分のみを使って評価値を計算できれば評価値計算の時間を短縮につながる、ヒューリスティックコンテストにおいてはこの差分計算は典型的なテクニックの1つとなっている。

*2:実際の提出したコードでは、問題文にある占いで得られる平均値\muの式に基づき\mathbf{v}\mathbf{X}を補正している。

*3:愚直全計算で顕著だがnが8の倍数の時にLoop回数が盛り返す。これはおそらくEigenの内部でAVX512命令を使っているのが関係していると思われる。

IIRフィルタの設計問題を焼きなまし法で解いてみる

これは数理最適化 Advent Calendar 2022 の 23 日目の記事です。

はじめに

ディジタルフィルタとは、離散時間信号に対し所望の周波数帯域の成分だけを通過させて他は阻止する機能のことを言う。中でもIIR (Infinite Impulse Response) フィルタはフィードバック機構を用いることにより無限長の畳み込み計算を行うフィルタのことを言う。 IIRフィルタは特定の制約を満たしつつ所望の周波数特性が得られるよう適切に設計する必要がある。フィルタ設計方法として様々な方法が提案されており、その中には最適化手法を利用したものも数多くある。

本記事は主に参考文献[1]の1章の内容を基に、最適化によるIIRフィルタ設計を実装してその感触を掴んでみたという内容になっている。 文献ではPSO(Particle Swarm Optimization)で最適化問題を解いているが、今回は手始めに(個人的に)実装ハードルが低い焼きなまし法による設計を試している。

ディジタルフィルタの基本

まずはディジタルフィルタについての基本事項を述べる。ここでは最低限のみ触れるので詳しくは参考文献[2][3]等を参照のこと。

フィルタの設計問題

IIRフィルタの処理は、離散時間の入力信号をx [  i ] 、出力信号をy [  i ] \ \  (i = 0, \cdots, \infty)とすると、以下のような差分方程式で表現できる。ただし、i \lt 0の場合はx[i] =0, y[i ] = 0とする。

\displaystyle y [i ] = -\sum_{m = 1}^{M} b_m y [i -  m ] + \sum_{n = 0}^{N} a_n x [i - n ] \tag{1}

このb_m \in \mathbb{R}, (m = 1, \cdots, M)a_n \in \mathbb{R}, (n = 0, \cdots, N)がフィルタ係数で、所望の周波数特性が得られるようにこの係数を決めるのがフィルタの設計にあたる。

式(1)の右辺第1項は過去の出力値を利用するフィードバック機構になっており、第2項が過去の入力値を利用するフィードフォワード機構になっている。このフィードバック機構があることによってフィルタの設計が不適切な場合、出力信号が無限に発散してしまう可能性がある。このことを不安定といい、逆に発散する心配のないフィルタは安定という*1。従って安定かつ所望の周波数特性が得られるようにIIRフィルタを設計する必要がある。

Z変換

連続信号を周期Tでサンプリングした離散信号はデルタ関数を用いて以下のように表すことができる。

\displaystyle x_s(t) = \sum_{n = 0}^{\infty} x(nT) \delta(t - nT) \tag{2}

このサンプリングされた信号をラプラス変換すると、

 \displaystyle
\begin{align}
\mathcal{L}[x_s(t)]  &= \int_{0}^{\infty} \sum_{n = 0}^{\infty} x(nT) \delta(t - nT) \exp(-st) dt \\
&= \sum_{n = 0}^{\infty} x(nT)  \exp(-nsT)
\end{align}
\tag{3}

となる。ここでz = \exp(sT)としたのがZ変換であり、次のような式になる*2

 \displaystyle X(z) = \sum_{i = 0}^{\infty} x [ i ]  z^{-i} \tag{4}

ラプラス変換ではs = j\omegaを代入することでフーリエ変換を得ていたが、Z変換も同様にz = e^{j\omega}を代入すれば離散時間フーリエ変換を得ることができる(\omegaは角周波数で周波数をfとすると\omega=2\pi fである)。そのため、この操作で周波数領域における特性を知ることができる。ただし、離散時間フーリエ変換ができるのは\sum_{i = 0}^{\infty}  | x [ i ]| \lt \inftyという条件を満たす場合のみであるため、これを満たさない信号に対してz = e^{j\omega}を代入する操作を行っても意味はない。

伝達関数

フィルタの周波数特性を見るには伝達関数が便利である。式(1)の両辺をZ変換にすると、

\displaystyle Y (z) = -\sum_{m = 1}^{M} b_m Y(z) z^{-m} + \sum_{n = 0}^{N} a_n X (z) z^{-n}\tag{5}

となる。ここでH(z) = Y(z)/X(z)伝達関数と呼ばれ、式(5)を整理してIIRフィルタの伝達関数を求めると以下のようになる。

\displaystyle H (z) =\frac{Y(z)}{X(z)} = \frac{\sum_{n = 0}^{N} a_n z^{-n}}{1 + \sum_{m = 1}^{M} b_m  z^{-m}} \tag{6}

この伝達関数に対して、 z = e^{j \omega}を代入すると伝達関数を離散時間フーリエ変換した結果が得られる。これを以下のようにH(\omega)と定義する。

\displaystyle H (\omega) =  \frac{\sum_{n = 0}^{N} a_n e^{-jn \omega}}{1 + \sum_{m = 1}^{M} b_m  e^{-jm\omega}} \tag{7}

H(\omega)はフィルタの周波数特性であり、このフィルタによって各周波数成分がどれくらい減衰するかという振幅特性| H(\omega) |と、どれくらい位相が遅れるかという位相特性\angle H(\omega)がわかる*3。振幅特性はどの周波数成分を通過させてどの周波数成分を阻止するかを左右するのでフィルタにおいて当然重要な特性だが、位相特性も重要で通過させたい帯域に関してはなるべく線形であることが望ましい。位相特性が線形でないと、出力信号の波形が歪んでしまうためである。

極と安定性

伝達関数の分母=0を満たす点を極という。フィルタが安定であるためには全ての極の絶対値が1未満である必要がある。そのため、分母のフィルタ係数b_m1 + \sum_{m = 1}^{M} b_m  z^{-m} = 0の解が全て絶対値1未満となるように定めなければならない。このことは最適化問題における制約条件となる。

焼きなまし法の適用

定式化

まず所望の周波数特性をH_D(\omega)とする。例えば、以下のような周波数特性を設定する。

\displaystyle
H_D(\omega) = 
\begin{cases}
 e^{-j \omega \tau_d} & \omega \in \Omega_p \\
0 & \omega \in \Omega_c
\end{cases}
\tag{8}

ここで、\tau_dは所望の群遅延、\Omega_pは通過させたい角周波数の集合、\Omega_cは阻止したい角周波数の集合である。 この理想的な周波数特性に近づくよう式(6)、(7)の係数を決めたいので、係数を並べたベクトル\mathbf{x} = [a_0, \cdots, a_N, b_1, \cdots, b_M ] を変数として、以下のような最適化問題を解くことになる。

\displaystyle
\begin{align}
\min_{\mathbf{x}}& \ \  J(\mathbf{x})\\
\rm{s.t.}& \ \  | p_i(\mathbf{x}) | < 1 \ \ (i = 1, \cdots, M) 
\end{align}
\tag{9}

ここでJ(\cdot)は理想周波数特性との誤差関数、p_i(\mathbf{x})はフィルタ係数が\mathbf{x}の時の極である。誤差関数は色々考えられるがよく使われるのは以下の2つ*4

  • 最大絶対誤差
\displaystyle
J(\mathbf{x}) =\max_{\omega \in \Omega} |H_D(\omega) - H(\omega; \mathbf{x})|
\tag{10}
  • 平均2乗誤差
\displaystyle
J(\mathbf{x}) = \frac{1}{|\Omega|}\sum_{\omega \in \Omega} |H_D(\omega) - H(\omega; \mathbf{x})|^2
\tag{11}

\Omegaは誤差計算に使用する周波数点の集合である。

さて、式(9)には極に関する制約条件があるがこのままでは扱いにくいので、制約条件違反をペナルティ項として組み込んだ新たな目的関数 F(\mathbf{x})を作り、こちらを焼きなましで解くことにする。

\displaystyle
\begin{align}
\min_{\mathbf{x}}& \ \ F(\mathbf{x}) =  \min_{\mathbf{x}} \ \ (J(\mathbf{x}) + c \phi(\mathbf{x})) \\

\phi(\mathbf{x}) &= 
\begin{cases}
 \max_{i} |p_i(\mathbf{x})|^2 & \max_i | p_i(\mathbf{x}) | \geq 1 \\
0 & \max_i | p_i(\mathbf{x}) | < 1
\end{cases}

\end{align}
\tag{12}

cはペナルティの強さを決めるハイパーパラメータになっている。

極の計算方法

式(12)のままだと更新の度に分母=0のM次方程式を解くなどして安定性を判断する必要がある。これでは効率が悪いため少し工夫する。まず簡単のためフィルタ係数の次数N, Mは偶数とする*5。そうすると式(6)は次のように変形できる。


\begin{align}
\displaystyle 
H (z) &= \frac{\sum_{n = 0}^{N} a_n z^{-n}}{1 + \sum_{m = 1}^{M} b_m  z^{-m}}  \\
&= a0 \frac{\prod_{n = 1}^{N/2} (1 - q_n z^{-1})(1 - \bar{q}_n z^{-1})}{\prod_{m = 1}^{M/2} (1 - p_m z^{-1})(1 - \bar{p}_m z^{-1})} \\
&= a0 \frac{\prod_{n = 1}^{N/2}  (1 - 2\rm{Re} (q_n) z^{-1} + |q_n|^{2} z^{-2}) } {\prod_{m = 1}^{M/2} (1 - 2\rm{Re}(p_m) z^{-1} + |p_m|^{2} z^{-2}) }
\end{align}

\tag{13}

ここで、a_0 \in \mathbb{R}, p_m, \bar{p}_m, q_n,  \bar{q}_n \in \mathbb{C}であり、\bar{p}_m, \bar{q}_np_m, q_n複素共役である。こう変形するとシステムの極はp_m, \bar{p}_mなので直接的にわかる。従って、\mathbf{x} = [a_0, p_1, \cdots, p_{M/2}, q_1, \cdots, q_{N/2} ] *6のように直接極を変数として扱う形で式(12)の最適化問題を解けば更新毎の安定性の判別は簡単にできる。

最適化アルゴリズム

ここでは実際に実装したアルゴリズムの流れを述べる。焼きなまし法自体の解説は割愛するので、詳しくは参考文献[5]等を参照。

まず、理想周波数特性の関数H_D(\cdot)、誤差計算に用いる周波数点集合\Omega、フィルタ次数M, N、ペナルティの強さc、繰り返し回数L、更新量の標準偏差\sigma、初期温度T、冷却率\alpha、使用する誤差関数J(\cdot)アルゴリズムのパラメータとして入力する。そして以下の処理を行う。

  1. \mathbf{x} = [a_0, p_1, \cdots, p_{M/2}, q_1, \cdots, q_{N/2} ] の各要素を正規分布\mathcal{N}(0, 1)に従う乱数で初期化する。
  2. \mathbf{x}'  = \mathbf{x} + \Delta \mathbf{x} を求める。ここで \Delta \mathbf{x}は各要素が正規分布 \mathcal{N}(0, \sigma^2)に従う乱数。
  3.  F(\mathbf{x}') \lt F(\mathbf{x} )であれば、 \mathbf{x} \leftarrow \mathbf{x}'に更新する(F(\cdot)は式(12)の目的関数)。 F(\mathbf{x}') \geq F(\mathbf{x})の場合は、\exp ( \frac{F(\mathbf{x}) -  F(\mathbf{x}')}{T})の確率で \mathbf{x} \leftarrow \mathbf{x}'に更新する。
  4. 温度を T \leftarrow \alpha Tとして更新する。
  5. 繰り返し回数がL回未満なら2.に戻る。L回に達したら、探索した中で最もF(\mathbf{x})が小さかった\mathbf{x}を出力して終了。

実装と実験

ここまで述べたことを実装して実験してみた。なお、ここに示す実験用コードはgithubにも上げている。

github.com

焼きなまし法の実装と実験

焼きなまし処理のコード

まずは焼きなまし法でフィルタ係数を求める処理の実装を以下に示す。

import pickle
import numpy as np
import matplotlib.pyplot as plt

def compute_mean_squared_error(H_, ref_h, omega_list):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H_(omg)

    e = np.real((ref_h - h) * np.conjugate(ref_h - h))
    return np.mean(e)

def compute_abs_max_error(H_, ref_h, omega_list):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H_(omg)
    e = np.abs(ref_h - h)
    return np.max(e)

# フィルタの周波数特性を計算する関数
def frequency_characteristic_func(a0, p_array, q_array, omega):
    denomi = 1.0
    for p in  p_array:
        denomi *= (1 - p * np.exp(-1j * omega)) *\
            (1 - np.conjugate(p) * np.exp(-1j * omega))

    numer = 1.0
    for z in q_array:
        numer *= (1 - z * np.exp(-1j * omega)) *\
            (1 - np.conjugate(z) * np.exp(-1j * omega))
    return a0 * (numer/denomi)

def obj_func(p_array, q_array, a0, ref_h, omega_list, c, error_func):
    H = lambda omega: frequency_characteristic_func(a0, p_array, q_array, omega)
    max_p2 = np.max(np.real(p_array * np.conjugate(p_array)))
    e = error_func(H, ref_h, omega_list)

    return e + (c * max_p2 if max_p2 >= 1.0 else 0.0)

# 焼きなまし法のメイン処理
def annealing(HD, omega_list, M, N, c, L, sigma, T, alpha, error_func,
              init_a0=None, init_p_array=None, init_q_array=None):

    # 理想周波数特性の関数から各周波数点に対する周波数特性を計算する
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = HD(omg)

    # 初期化
    if init_a0 is None:
        a0 = 1.0
    else:
        a0 = init_a0
    if init_p_array is None:
        p_array = np.random.randn(M//2) + 1j * np.random.randn(M//2)
    else:
        p_array = np.copy(init_p_array)
    if init_q_array is None:
        q_array = np.random.randn(N//2) + 1j * np.random.randn(N//2)
    else:
        q_array = np.copy(init_q_array)

    # 目的関数F
    F = lambda p_array_, q_array_, a0_: \
        obj_func(p_array_, q_array_, a0_, h, omega_list, c, error_func)

    now_cost = F(p_array, q_array, a0)

    # 最終的に出力するフィルタ係数を初期化
    best_p_array = np.copy(p_array)
    best_q_array = np.copy(q_array)
    best_a0 = a0
    best_cost = now_cost

    for i in range(L):
        # ロールバック用に保持
        save_a0 = a0
        save_p_array = np.copy(p_array)
        save_q_array = np.copy(q_array)

        a0 += np.random.randn() * sigma
        p_array += (np.random.randn(M//2) + 1j * np.random.randn(M//2)) * sigma
        q_array += (np.random.randn(N//2) + 1j * np.random.randn(N//2)) * sigma

        tmp_cost = F(p_array, q_array, a0)
        if np.exp((now_cost - tmp_cost)/T) > np.random.rand():
            now_cost = tmp_cost
            print("now_cost=", now_cost, "T=", T, "i=", i)
            if now_cost < best_cost:
                best_a0 = a0
                best_p_array = np.copy(p_array)
                best_q_array = np.copy(q_array)
                best_cost = now_cost
                print("best_cost = ", best_cost)
        else:
            # ロールバック
            a0 = save_a0
            p_array = np.copy(save_p_array)
            q_array = np.copy(save_q_array)
        T *= alpha

    return best_a0, best_p_array, best_q_array, best_cost

基本的には先述のアルゴリズムに書いた通りの実装になっている。annealingという関数が実際の焼きなましの処理を行う関数である。F(\cdot)の計算などを工夫すればもっと効率化できるかもしれない。

実験

それでは実際に理想的な周波数特性を設定し、この焼きなまし処理を用いてフィルタ係数を求めてみる。 理想的な周波数特性は式(8)の形式において、通過域の集合を\Omega_p = \{\omega |0.4 \pi \leq \omega \leq 0.6 \pi \} 、阻止域の集合を\Omega_c = \{\omega | 0 \leq \omega \leq 0.3 \pi, 0.7 \pi \leq \omega \leq \pi \} 、群遅延を\tau_d = 12と設定して実験した。なお、\Omega_p, \Omega_cを見てわかる通り、急激な変化を避けるため通過域と阻止域の間に特に周波数特性を指定しない遷移域も設けている。この理想周波数特性のグラフは以下のようになる。

理想周波数特性

焼きなましのパラメータについては、

M = 8, N  = 8, c = 5, L = 10000, \sigma = 0.01, T=10, \alpha = 0.998

とした。誤差関数については両方の方法を試している。

また、問題の性質上、得られる結果が初期値や乱数に依存して大きく変わるので、この実装では初期値を変えながら100回最適化処理を行い、最もF(\mathbf{x})が小さくなったものを採用している。

このような実験を行うコードは以下のようになっている。

def visualize_characteristic(H, omega_list, save_file_name):
    h = np.zeros(len(omega_list), dtype=np.complex)
    for i, omg in enumerate(omega_list):
        h[i] = H(omg)

    fig = plt.figure()
    amp_h = np.abs(h)
    ax1 = fig.add_subplot(2, 1, 1)
    ax1.plot(omega_list, amp_h, marker=".")
    ax1.set_ylim(-0.1, 1.1)
    ax1.set_ylabel("amplitude")

    angle_h = np.angle(h)
    ax2 = fig.add_subplot(2, 1, 2)
    ax2.plot(omega_list, angle_h, marker=".")
    ax2.set_ylim(-np.pi - 0.2, np.pi + 0.2)
    ax2.set_xlabel("omega")
    ax2.set_ylabel("phase")
    fig.savefig(save_file_name)
    plt.close()

if __name__ == '__main__':
    np.random.seed(0)

    # 理想的な周波数特性の関数
    HD = lambda omg: np.exp(- 1j * 12 * omg) \
        if omg >= 2 * 0.2 * np.pi and omg <= 2 * 0.3 * np.pi  else 0.0

    # 誤差関数に利用する誤差計算に使用する周波数点の集合
    omega_list1 = np.linspace(0, 2 * 0.15 * np.pi, 200)
    omega_list2 = np.linspace(2 * 0.2 * np.pi, 2 * 0.3 * np.pi, 200)
    omega_list3 = np.linspace(2 * 0.35 * np.pi, np.pi, 200)
    omega_list = np.concatenate([omega_list1, omega_list2, omega_list3])

   # 理想的な周波数特性のグラフを作成
    visualize_characteristic(HD, omega_list, "ideal.png")

    # 100回焼きなましをして最もよかったものを採用
    best_cost = 1e9
    for i in range(100):
        a0_, p_array_, q_array_, cost =\
            annealing(HD, omega_list, 8, 8, 5,
                      10000, 0.01, 10, 0.998,
                      compute_mean_squared_error) # 実験する誤差関数を指定
        if cost < best_cost:
            a0 = a0_
            p_array = np.copy(p_array_)
            q_array = np.copy(q_array_)
            best_cost = cost

    # 求めたフィルタの周波数特性の関数
    H = lambda omega: frequency_characteristic_func(a0, p_array, q_array, omega)

    # 求めたフィルタの周波数特性のグラフを作成
    visualize_characteristic(H, np.arange(0, np.pi, 0.01), "solution.png")

    # 制約条件を満たしていることを確認
    for p in p_array:
        print("|p|=", np.abs(p))

    print("final cost=", best_cost)

    # 求めたフィルタ係数を保存
    save_dict = {"a0":a0, "p_array":p_array, "q_array":q_array}
    with open('filter_coef.pickle', 'wb') as f:
        pickle.dump(save_dict, f)

これを実行して得たIIRフィルタの周波数特性は以下のようになった。

  • 誤差関数が最大絶対誤差の場合

求めたIIRフィルタの周波数特性(最大絶対誤差)

  • 誤差関数が平均2乗誤差の場合

求めたIIRフィルタの周波数特性(平均2乗誤差)

目視の印象だが通過域については最大絶対誤差で求めた方が平均2乗誤差に比べて理想の周波数特性に近いように見える。一方で、阻止域については最大絶対誤差だと振幅特性が0付近まで落ち切っておらず平均2乗誤差の方が優位に見える。

フィルタ処理の実装と実験

フィルタ処理のコード

得られたフィルタ係数を用いて実際の離散信号に対してフィルタをかけてみる。式(13)を展開して式(6)の形式にした上で式(1)の計算でフィルタをかけても良いが、ここでは式(13)の形式のまま2次のフィードバック(もしくはフィードフォワード)処理を縦続に接続することでフィルタ処理を実現している。つまり、2次フィルタの出力をそのまま次の2次フィルタに入力するような形にしている。これについては実装を見た方がわかりやすいだろうと思うので早速実装を貼る。

class IIRFilter:
    def __init__(self, a0_, p_array_, q_array_):
        self.a0 = a0_
        self.p_array = np.copy(p_array_)
        self.q_array = np.copy(q_array_)

    def _feedback_filter(self, x, p):
        b1 = -np.real(p) * 2
        b2 = np.real(p * np.conjugate(p))
        y = np.zeros(x.shape[0])
        for i in range(x.shape[0]):
            s = x[i]
            if i - 1 >= 0:
                s += -b1 * y[i - 1]
            if i - 2 >= 0:
                s += -b2 * y[i - 2]
            y[i] = s
        return y

    def _feedforward_filter(self, x, q):
        a1 = -np.real(q) * 2
        a2 = np.real(q * np.conjugate(q))
        y = np.zeros(x.shape[0])
        for i in range(x.shape[0]):
            s = x[i]
            if i - 1 >= 0:
                s += a1 * x[i - 1]
            if i - 2 >= 0:
                s += a2 * x[i - 2]
            y[i] = s
        return y

    def filtering(self, x):
        y = np.copy(x)
        for q in self.q_array:
            y = self._feedforward_filter(y, q)
        for p in self.p_array:
            y = self._feedback_filter(y, p)
        y *= self.a0
        return y
実験

この実装を使っていくつかの信号にフィルタをかけてみる。この実験用コードについては上述のgithubを参照。

  • \cos (0.1 \pi t)

まずは低周波数帯の遮断域にある周波数の信号をフィルタリングしてみる。以下のように減衰された信号が得られるが、誤差関数が最大絶対誤差の場合は平均2乗誤差の場合に比べて減衰が甘い。この結果は上述の振幅特性通りの結果と言える。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

  • \sin (0.4 \pi t)

次に通過域ギリギリにある周波数の信号を試す。誤差関数が最大絶対誤差の場合はt=0付近以外ほとんど減衰のない出力が得られているが、平均2乗誤差では少し減衰した結果になっている。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

  • \sin (1.0 \pi t)

高周波数帯の遮断域にある周波数の信号。低周波数帯の時と同様に、誤差関数が最大絶対誤差の場合は減衰しきれていない結果になっている。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

  • \cos (0.1 \pi t) + \sin (0.4 \pi t) + \sin (1.0 \pi t)

最後に混合した信号で試す。通過域の成分のみ抽出できているように見える。波形としてきれいなのは、遮断域をきちんと遮断できている平均2乗誤差の方か。

最大絶対誤差によるフィルタリング結果

平均2乗誤差によるフィルタリング結果

所感

上でも述べたが良い解が得られるかどうかは初期値と乱数にかなり依存し、あまり好ましくない局所最適解が得られるケースも多々あった。特に誤差関数として最大絶対誤差を選択した場合では、振幅特性がほとんど0になるような解が得られるケースも多かった。そのため今回は初期値を変えながら何度も焼きなましを行う方法で対処したが、そもそも焼きなましで解くには非線形性が強すぎる問題だった可能性がある。 そう考えると、PSOのように複数の個体を有する最適化手法を使うのは妥当な選択なのかもしれない。この辺は次の記事辺りでも比較してみたいところ。

また、フィルタの次数M, Nを増やすほど最適化が難しくなる。変数が増えるので当然と言えば当然だが、今回の試行ではM=8, N=8(9変数)程度がまともな解が得られる限界だったのは少し残念だった。

加えて、通過域と阻止域の境界となる不連続な領域は遷移域として誤差関数の計算に含めないことも適切な最適化を得るのに思いのほか寄与した。不連続に変化する点は線形フィルタの周波数特性では得ることができず、より最適化が難しい問題になるのであろうと予想する。

最後に

焼きなまし法自体は比較的簡単に様々な最適化問題に適用できるため、今回のように最適化問題の感触を掴むのには良いかもしれない。ただし、そこからより良い解をより速く得たい場合は工夫や調整が必要だろう*7。今回ももう少し工夫すればより良いフィルタが得られた可能性もある。

IIRフィルタの設計には既存のアナログフィルタから変換する等最適化を利用しない方法もある。この方法で得られたフィルタ係数を初期値として最適化するのも面白いかもしれない。

参考文献

[1] 今回の記事を書くきっかけになった書籍。今回は1章の内容を参考にしているが、他の章には他にも凸関数に近似して解く方法や、他の信号処理の問題を最適化手法を用いて解く事例が載っていて面白い。

[2] Z変換など信号処理における基本事項は以下の書籍を参考にしている*8

[3] 信号処理の基本に関してはこちらも参考にした。ネットで見れる最も丁寧なディジタル信号処理解説の1つだと思う。

www.ic.is.tohoku.ac.jp

[4] 電子情報通信学会が公開してるディジタルフィルタの解説。こちらの2章も今回の記事を書くにあたって参考にした。

www.ieice-hbkb.org

[5] お世話になっている最適化の書籍。この書籍には焼きなまし法についての解説もある。

*1:第2項のみのフィルタはFIR(Finite Impulse Response)フィルタと呼ばれる。FIRフィルタはフィードバック機構がないため常に安定だが、一般にIIRフィルタの方が少ない次数で良い周波数特性が得られる。

*2:式(2)、(3)では便宜上xを時刻の関数で表していたが、式(4)では添え字の形式で表現している。つまり、t = 0, T, 2T, \cdotsがそれぞれi = 0, 1, 2, \cdotsに対応している。

*3:このような操作で周波数特性が得られるのはフィルタが安定の場合に限る。フィルタが不安定の場合は伝達関数の絶対値総和が発散し、そもそも離散時間フーリエ変換が成り立たたない。

*4: H(\omega; \mathbf{x})は式(7)のH(\omega)と同じものだが、\mathbf{x}の変数であることを明示的に示すためこのように表現している。

*5:奇数の次数にもできるが少し面倒なので、ここでは偶数にしている。

*6: \mathbf{x} の中に実数の要素と複素数の要素が混在していてあまりよろしくないが、面倒なのでこういう表記をしている。

*7:競プロのマラソンマッチ(ヒューリスティックコンテスト)と呼ばれる形式でも焼きなましが有効になる問題が多々ある。しかし、それでもハイスコアを狙うにはただ単に焼きなましを適用するだけでなく、初期値の工夫、定式化の工夫、高速化等様々な工夫が必要になる。

*8:10年以上前に購入して個人的には割とお世話になっている書籍だし内容も分かりやすいと思うが、どうやら今は中古でしか手に入らないらしい。

カーネル関数は何故『半』正定値でよいのか

はじめに

タイトルの通りカーネル関数の話。内積のことを考えるとカーネル関数から作られるグラム行列は正定値でないとダメな気がしたがそんなことはなくて半正定値であればよい、それは何故か?という話。

いくつかの書籍に書いてある話なのでわざわざ記事にしなくてもよいかもしれないが、自分が抱いた疑問やその解消に至るまでの経緯などを含めて備忘録として残しておく。

なお、カーネル関数については以前書いた以下の記事でも触れているので参考までに。

yamagensakam.hatenablog.com

カーネル関数と再生性

本題に入る前に前提となるカーネル関数の定義や対象とする空間および再生性という性質について簡単に述べておく。

まず集合\mathcal{X}上の2変数関数であって、対称性と半正定値性を有する関数kカーネル関数と呼ぶ。 これは、任意 x_1 , \cdots, x_n \in \mathcal{X}に対しij列目がk(x_i, x_j)のグラム行列\mathbf{K} を作ると、それが半正定値対称行列になることを意味する。

また、カーネル関数kについて、第2引数をx \in \mathcal{X}で固定した1変数関数をk(\cdot, x) = k_x(\cdot)  = k_xと表記する。

そして、この1変数関数を用いて内積空間\mathcal{V}を以下のように定義する。

\displaystyle 
\mathcal{V} = \left \{ \sum_{i=1}^{n} a_i k_{x_i} \  \middle| \ x_i \in \mathcal{X}, \ a_i \in \mathbb{R}, \ n \in \mathbb{N} \right\}
\tag{1}

この空間の内積f = \sum_{i = 1}^{n} a_i k_{x_i}, \ g  =  \sum_{j = 1}^{m} b_j k_{y_j} \in \mathcal{V}として、

\displaystyle 
\langle f, g \rangle =  \sum_{i=1}^{n} \sum_{j=1}^{m} a_i b_j k(y_j, x_i)
\tag{2}

と定義する。後述の零関数の例にもある通りf \in \mathcal{V}の表現方法は複数存在する可能性がある(つまり、a_ix_iが異なっていても同じfになることがある)。それでも式(2)で定義する内積はwell-definedであり、f, gの表し方によらず1つに定まる。

このように定義した内積空間では、任意のf \in \mathcal{V}, \ x \in \mathcal{X}に対し以下のような再生性という性質を有する。

\displaystyle 
f(x) =  \langle f, k_x \rangle
\tag{3}

これは関数fxを代入する操作が内積の形で表せることを意味し、この性質から様々な定理が導かれる。

ちなみに、このままでは完備性を有しておらずヒルベルト空間とは呼べない(と思われる*1)。この\mathcal{V}に対し完備化という操作を施し完備性を有した空間は再生核ヒルベルト空間になるが、今回の内容には絡まないのでこれ以上は触れない。

疑問点

このように定めた内積空間\mathcal{V}において

\displaystyle 
\langle f, f \rangle = \sum_{i = 1}^{n} \sum_{j = 1}^{n} a_i  a_j k(x_i, x_j) = \mathbf{a}^T \mathbf{K} \mathbf{a}
\tag{4}

を考える(\mathbf{a}a_ii番目の要素に持つベクトル)。カーネル関数が半正定値ということはグラム行列\mathbf{K}が半正定値で \rm{rank} (\mathbf{K}) \lt nにもなり得る。ということは、 \rm{rank} (\mathbf{K}) \lt nであれば\langle f, f \rangle  = 0となる\mathbf{0}以外の\mathbf{a}が存在することになる。

ところが、内積の定義を考えると\langle f, f \rangle = 0 \Leftrightarrow f = 0のはずで、それにもかかわらず\mathbf{a}\mathbf{0}でなくとも式(4)を0にできるのはなんか奇妙、というのが今回の論点。

疑問の解消

そもそもの話

考えてみればf = 0とはどういうことだろう?f = \sum_{i = 1}^{n} a_i k_{x_i}a_i = 0\ ( i = 1, \cdots, n) というなら分かるが、それ以外に意味することはあるのだろうか?

これはベクトル空間に立ち返るとfは零ベクトルであることを意味する。零ベクトルなので任意のg \in \mathcal{V}に対して、f + g = g + f = gとなればよい。 関数なので表現方法が異なっていても入出力関係が変わらなければよい。すなわち、任意のx\in \mathcal{X}に対してf(x) + g(x) = g(x) + f(x) = g(x)であればfは零ベクトルと言える*2。これを達成するfは何を入力しても0を出力する関数、つまり零関数である。fが零関数であれば、a_iが全て0である必要もない。

\langle f, f \rangle = 0 \Leftrightarrow f = 0の証明

f=0 \Rightarrow \langle f, f \rangle = 0は自明なので \langle f, f \rangle = 0 \Rightarrow f = 0のみ示す。 まず、f, g \in \mathcal{V}について以下のコーシー・シュワルツの不等式が成り立つことを示す。

\displaystyle 
| \langle f, g \rangle |^2 \leq  \langle f, f \rangle \langle g, g \rangle
\tag{5}

これは\langle tf + g, tf + g \rangle  \ \  (t \in \mathbb{R}) を考えると、

\displaystyle 
\begin{align}
\langle tf + g, tf + g \rangle  = t^2 \langle f, f \rangle + 2t \langle f, g \rangle + \langle g, g \rangle \geq 0

\end{align}

であり、これが任意の実数tに対して成立するためには判別式が0以下、すなわち4\langle f, g \rangle^2  - 4 \langle f, f \rangle  \langle g, g \rangle  \leq 0となる必要があるため、これを展開して式(5)が導かれる。

また、任意のx \in \mathcal{X}に対し|f(x)|^2 = 0であればf = 0である。|f(x)|^2を式(3)の再生性と式(5)を用いて展開すると、

\displaystyle 
| f(x) |^2 = | \langle f, k_x \rangle |^2 \leq  \langle f, f \rangle \langle k_x , k_x  \rangle
\tag{6}

となる。そして、\langle f, f\rangle = 0なので 0 \leq |f(x)|^2 \leq 0となり f = 0が導かれる。

具体例で検証

このことと式(4)より\mathbf{a}^T \mathbf{K} \mathbf{a} = 0となるa_i, \  x_i  \ \    (i =1, \cdots, n) の組み合わせではf = \sum_{i = 1}^{n} a_i k_{x_i}は零関数であり、任意のy \in \mathcal{X}についてf(y) = \sum_{i = 1}^{n} a_i k(y, x_i) = 0となる。

これをカーネル関数として2次の多項式カーネル k(\mathbf{y}, \mathbf{x}) = (1 + \mathbf{x}^T \mathbf{y})^2 \ \  (\mathbf{x}, \mathbf{y} \in \mathbb{R}^d) を採用した場合を例にとって考えてみると、\mathbf{a}^T \mathbf{K} \mathbf{a} = 0であれば任意の\mathbf{y}に対してf(\mathbf{y}) = \sum_{i = 1}^{n} a_i (1 +\mathbf{x}_i^T \mathbf{y})^2 = 0が成り立つ。

ここまで見てもまだ直感と反する気がするので、一応このことをスクリプトを書いて確かめてみる。多項式カーネルなのでd \ll nなら \rm{rank} (\mathbf{K}) \lt nとなる。そして、 \mathbf{K} の零固有値に対応する固有ベクトル\mathbf{a}とすれば\mathbf{a}^T \mathbf{K} \mathbf{a} = 0となる。この\mathbf{a}を用いて作ったfは果たして適当な\mathbf{y}に対してf(\mathbf{y}) = 0となるのか、ということを確かめるスクリプトが以下。

import numpy as np

def func(a, X, y, k):
    res = 0
    for i in range(a.shape[0]):
        res += a[i] * k(y, X[i, :])
    return res

if __name__ == '__main__':
    np.random.seed(3407)
    d = 6
    n = 30
    X = np.random.randn(n, d)

    # 2次多項式カーネル
    k = lambda x, y: (1 + x @ y) ** 2

    K = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            K[i, j] = k(X[i, :], X[j, :])

    print ("rank(K) =", np.linalg.matrix_rank(K))

    # Kの零固有値に対応する固有ベクトルを抽出
    # このベクトルaはa @ K @ a = 0を満たす
    l, U = np.linalg.eigh(K)
    a = U[:, np.abs(l) < 1e-12]

    # 入力を適当に生成
    y = -50 + 100 * np.random.rand(d)

    for i in range(a.shape[1]):
        f = lambda y : func(a[:, i], X, y, k)
        print("f" + str(i) + "(x) =", f(y))

結果としては、

rank(K) = 28
f0(x) = -1.4438228390645236e-11
f1(x) = -1.4743761767022079e-12

演算誤差もあるので完全に0とはならないが、ほぼ0になったといえる。また乱数のseedを変えて何度か試しても同じくらいの値になるので、とりあえず納得。

おわりに

最後に参考文献。冒頭に述べた以前書いた記事でも同様の参考文献を挙げたが今回もお世話なったので改めて記載。

*1:これは自分の中でまだ理解が怪しい。式(1)の定義だと有限個の和なので完備になったりしないか?n \rightarrow \inftyとして解析するときに完備性の議論が必要なる?

*2:関数同士の和は(f + g)(x) = f(x) + g(x)としていいのかという話もある。おそらく一般的にも関数の和はこのように定義するのだと思うが、少なくとも今回の場合は、再生性より(f + g)(x) = \langle f + g, k_x \rangle = \langle f, k_x \rangle + \langle g, k_x \rangle = f(x) + g(x)として得られる。

HTTFで超球面上の最適化を試してみた話

はじめに

AtCoder上で開催されたフューチャー社主催のヒューリスティックコンテストHACK TO THE FUTURE2022(HTTF)に参加した*1。問題等詳細は以下を参照。

atcoder.jp

この問題は大きく分けて各作業者のスキル推定と作業スケジューリングという2つの要素があるが、その中でもスキル推定の方に以前から気になってたリーマン多様体上の最適化なるものが使えそうだったので試してみたというのが今回の話。実のところ自分も多様体は何ぞやとか数学的なことはよく分かっていないが、今回行った単位超球面上の最急降下法であればそこまで難しくないし、以下の資料を参考に比較的簡単に実装できる。

http://www.orsj.or.jp/archive2/or60-9/or60_9_549.pdf

結果としては割とうまく推定できるが普通の最急降下法でやった場合とスコア的にはさほど差がなく、最終的にはスコアが高く出た普通の最急降下法を採用した(誤差範囲の気もする)。

以下、具体的行ったことを書いていく。

スキル推定問題の定式化

問題全体は上記のリンク先を見ていただくとして、ここではスキル推定問題について定式化をしておく。スキルは作業者間で独立なので、ここでは1人の作業者のスキルを推定する問題を定式化する。 作業者はK種類のスキルを有しており、そのk番目のスキルをs_k \geq 0とする*2。また、タスクの要求スキルも同様にK種類あり、タスクiの要求スキルをd_{ik} \geq 0とする。そして、実際にその作業者がタスクiを実行したときにかかった時間はt_iであり、この時間は \sum_{k = 1}^{K}\max(0, d_{ik} - s_k) = 0ならt_i = 1 \sum_{k = 1}^{K}\max(0, d_{ik} - s_k) > 0ならこの値に-3 \sim 3の乱数を加えたもので得られる(この場合も0以下になるならt_i = 1とする)。n個タスクが作業者に割り振られ、\left\{(\mathbf{d}_1, t_1), \cdots, (\mathbf{d}_n, t_n) \right\}が得られているので、この集合から作業者のスキル\mathbf{s}を推定したい。

ということで、このスキル推定問題を2乗誤差最小化問題として以下のように定式化する。

\displaystyle
{\rm minimize} \ \  E(\mathbf{s}) = {\rm minimize} \ \  \frac{1}{2n} \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - |s_k|) \right\}^2
\tag{1}

ここで制約付き最適化になるとややこしいので、s_k \geq 0という制約はつけず負値を取ることを許して、式の中で絶対値を取るようにしている。

通常の最急降下法でスキル推定

最終的に採用したのはこっちの方。式(1)の中にmaxやら絶対値やらが入っているので一部の人から怒られそうだが、気にせず\mathbf{s}微分して勾配を求める*3

\displaystyle
\nabla E(\mathbf{s}) =   \frac{1}{n} \mathbf{h} \odot \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - |s_k|) \right\} \odot \mathbf{g}_i
\tag{2}

ただし、\mathbf{h} \mathbf{g}_i K次元のベクトルで、そのk番目の要素をそれぞれh_k, g_{ik}とすると

\displaystyle
h_k = \left\{
\begin{array}{ll}
1 & (s_k > 0) \\
0 & (s_k = 0) \\
-1 & (s_k < 0)
\end{array}
\right., 
\ \ \ 
g_{ik} = \left\{
\begin{array}{ll}
1 & (d_{ik} - |s_k| > 0) \\
0 & (d_{ik} - |s_k| \leq 0)
\end{array}
\right.

である。また、\odotアダマール積を表す。スキル\mathbf{s}はこの勾配を用いて、

\displaystyle
\mathbf{s}^{(l + 1)} = \mathbf{s}^{(l)} - \alpha \ \nabla E(\mathbf{s}^{(l)})
\tag{3}

という更新式を適当な回数まわして推定した。本来なら収束判定を勾配のノルムとかで行うべきだが、手を抜いてループ回数を固定している。また、ステップ幅\alphaについても本来なら線形探索などで適切な値を求めるべきだが、これまた手を抜いて適当な値に固定している。

超球面上の最急降下法でスキル推定

ここからが本題。すぐ後に示すように式(1)を少し変形するとノルム1制約の最適化問題が出てくるので、そこに超球面上の最急降下法を使った。

スキルの成分分解

問題文をよく見ると真のスキルの生成方法が書かれている。これを見ると、\mathbf{s} = \frac{{\rm randdouble}(20, 60) | \mathbf{s}' |}{\parallel \mathbf{s}' \parallel_2} \ \ (\mathbf{s}' \sim \mathcal{N}(\mathbf{0}, \mathbf{I}))という式で生成されている。従って、ノルムが1のベクトル\mathbf{w}と実数q \ \ (20 \leq q \leq 60)を用いて、\mathbf{s} = q|\mathbf{w}|と書くことができる。 これを使うと式(1)は

\displaystyle
\begin{align}
{\rm minimize} \ \  E_{d}(\mathbf{s}) = {\rm minimize} &\ \  \frac{1}{2n} \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - q|w_k|) \right\}^2\\
{\rm s.t.} &\ \  \sum_{k=1}^K w_k^2 = 1 \\
& 20 \leq q \leq 60
\end{align}
\tag{4}

というq\mathbf{w}に関する最適化に変換される。 この式(4)の最適化のため、\mathbf{w}を固定してqを最適化する処理と、qを固定して\mathbf{w}を最適化する処理を交互に繰り返す方法を取った。

qの最適化

\mathbf{w}を固定すると、式(4)はqについて凸な関数になる。また、qは1次元の実数。ということで三分探索で最適解を求めた。ただ、タイムリミットの関係上三分探索のループは5回にしてるので、かなり粗い探索になっている。

\mathbf{w}の最適化

\mathbf{w}はノルムが1という制約条件があり、これは\mathbf{w}が単位超球面上に存在する制約だといえる。 このような制約条件つき最適化問題に対しては様々なアルゴリズムが提案されているが*4、今回は冒頭で挙げた資料の通り超球面上での制約条件なしの最適化問題とみなす手法を実装した。

まずは一応式(4)の\mathbf{w}に関する通常の勾配を書いておくと、

\displaystyle
\begin{align}
\nabla E_{d}(\mathbf{w}) =  \frac{q}{n} \mathbf{h}' \odot \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - q|w_k|) \right\} \odot \mathbf{g}_i'
\end{align}
\tag{5}

ただし、

\displaystyle
h_k' = \left\{
\begin{array}{ll}
1 & (w_k > 0) \\
0 & (w_k = 0) \\
-1 & (w_k < 0)
\end{array}
\right., 
\ \ \ 
g_{ik}' = \left\{
\begin{array}{ll}
1 & (d_{ik} - q|w_k| > 0) \\
0 & (d_{ik} - q|w_k| \leq 0)
\end{array}
\right.

となる。この勾配方向に式(3)のように動かしても、超球面のことは全く考慮されていないので明後日の方向に行ってしまう可能性がある。ということで、

  • なるべく超球面から離れないよう現在点における接平面に沿った方向に動かす
  • それでも超球面からはみ出るので超球面上に戻す(レトラクション)

ということを繰り返すのが当該のアルゴリズム。以下もう少し具体的に書く。

動かす方向

今超球面上の点\mathbf{w}^{(l)}にあって、そこから動かす方向を決める。\mathbf{w}^{(l)}における接平面の法線ベクトルは\mathbf{w}^{(l)}である。この法線ベクトルが張る空間の直交補空間は接平面に平行な部分空間(接ベクトル空間)であり、\nabla E_{d}(\mathbf{w}) のその部分空間へ直交射影したベクトルを\mathbf{d}とすると、

\displaystyle
\mathbf{d} = \nabla E_d(\mathbf{w}^{(l)}) - \langle \mathbf{w}^{(l)}, \nabla E_d(\mathbf{w}^{(l)}) \rangle \mathbf{w}^{(l)}
\tag{6}

となる。これを用いて、

\displaystyle
\mathbf{w'}^{(l)} = \mathbf{w}^{(l)} - \alpha \mathbf{d}
\tag{7}

にまずは動く。こうすれば、とりあえずは接平面に沿った移動になるので、超球面から離れないという観点で言うと少なくとも-\nabla E_d(\mathbf{w}^{(l)})方向に移動するよりはマシになるはず、というのが直感的な理解。なお、\alphaはステップ幅でこちらも本来は(曲線)探索するべきだが実装では手を抜いて固定値にした。

超球面へのレトラクション

式(7)で得た\mathbf{w'}^{(l)} は、それでもやはり超球面からははみ出ている。従ってこれを超球面上に戻す必要がある。この操作をレトラクションと呼ぶ。これは単純に\mathbf{w'}^{(l)} のノルムで割ればよい。よって、

\displaystyle
\mathbf{w}^{(l + 1)} = \frac{\mathbf{w'}^{(l)}}{\parallel \mathbf{w'}^{(l)} \parallel_2}
\tag{8}

として\mathbf{w}を更新できる。

結果

seed=0の時のビジュアライザ動画を貼りたかったけどはてなブログに動画がアップロードできなかったので、最後の状態を画面キャプチャした画像を貼る。

うーん、どっちがいいとも悪いとも言いにくい。結局のところそこまで変わらないので、処理時間的にも短い通常の最急降下法を最終的に採用したわけだが。気持ちとしては後者の方がテクくて好きだったので何とかして採用したかった。。。

大半の人はスキル推定に焼きなましやMCMCを使っているので、実際のところそちらの方がいいのかもしれない。

あと、これを書いてて思ったが、振られたタスクの要求スキルが小さいと各スキルは過小評価されるので、最初のうちは満遍なくいろいろなタスクを振ってスキル推定の精度を向上させた方がよかったかもしれない。

ちなみに

作業スケジューリングについては、重み付き二分マッチング問題を誰かがタスクを終える度に実行して作業の割り当てを行った。この際、辺の重みは予測される作業時間に加えて、なるべくそのタスクを終えることで次に実行できるタスクが増えるように調整している。 スキル推定はそこまで重要じゃない気もするので、作業スケジューリングの方にもっと時間を割いて考えるべきだったかもしれない。

所謂ヒューリスティックコンテスト形式のコンテストはAtCoder Heuristic Contestが始まってからまともに参加するようになったが、今回のように連続最適化の手法も使えるコンテストがたまにあって中々楽しい。

*1:終結果はスコア5,664,033で129位。

*2:スキルとして取りうる値は整数だが、ここでは実数値として考え、後から四捨五入する方法を取った。

*3:一応、最初は絶対値ではなく2乗で0以上を保証するなど微分可能な式を立ててやってたが結果的には今の定式化の方がよかった。今や深層学習なんかでReLU等をガンガン微分してるのでまぁ。。。

*4:以前書いた内点法とかyamagensakam.hatenablog.com

超平面と点の距離の求め方を少し抽象的に書いてみる

はじめに

以下の書籍を読み関数解析のお気持ちを理解しつつあるので、今回は備忘録がてら超平面と点の距離を抽象ヒルベルト空間の観点で記述してみる。

なお話を簡単にするために、超平面は必ず原点を通るものとする。

n次元実ベクトル空間\mathbb{R}^{n}における超平面と点の距離

まずは通常の線形代数の範疇で超平面と点の距離を記述する。超平面上の点を\mathbf{x} \in \mathbb{R}^{n}とすると、超平面は原点を通るとしているのでその方程式は適当な法線ベクトル\mathbf{v} \in \mathbb{R}^{n}を用いて、

\displaystyle
\mathbf{v}^T \mathbf{x} = 0
\tag{1}

と表される。この超平面と任意の\mathbf{x}_0 \in \mathbb{R}^nの距離だが、これは\mathbf{x}_0を超平面上に射影した点\bar{\mathbf{x}}_0\mathbf{x}_0の距離にあたる。超平面上の部分空間と\mathbf{v}が張る空間は直交しているため\bar{\mathbf{x}}_0\mathbf{x}_0の関係は\alpha \in \mathbb{R}を用いて

\displaystyle
\mathbf{x}_0 = \bar{\mathbf{x}}_0 + \alpha \mathbf{v}
\tag{2}

が成り立つ。ということで超平面と\mathbf{x}_0の距離d

\displaystyle
d = \parallel \mathbf{x}_0 - \bar{\mathbf{x}}_0 \parallel_2 = \parallel  \bar{\mathbf{x}}_0 + \alpha \mathbf{v} -  \bar{\mathbf{x}}_0 \parallel_2 =  |\alpha|  \parallel \mathbf{v} \parallel_2
\tag{3}

となる。\alphaだが、これは式(2)の両辺に\mathbf{v}^Tをかけて、\bar{\mathbf{x}}_0 は超平面上の点なので\mathbf{v}^T\bar{\mathbf{x}}_0 = 0であることに注意すると、

\displaystyle
\alpha = \frac{\mathbf{v}^T \mathbf{x}_0}{\mathbf{v}^T\mathbf{v}} = \frac{\mathbf{v}^T \mathbf{x}_0}{\parallel \mathbf{v} \parallel_2^2}
\tag{4}

となるので、後はこれを式(3)に代入して

\displaystyle
d =\frac{|\mathbf{v}^T \mathbf{x}_0|}{\parallel \mathbf{v} \parallel_2} = \frac{|\mathbf{v}^T \mathbf{x}_0|}{\sqrt{\mathbf{v}^T \mathbf{v}}} 
\tag{5}

を得る。

ヒルベルト空間の定義と諸定理

次にこれをヒルベルト空間上の議論へ抽象化してみる。が、そもそもヒルベルト空間とは?その超平面とは?距離とは?という話なので、本題に入る前に今回の議論に必要な定義および定理を紹介する。

ヒルベルト空間の定義

まずヒルベルト空間\mathcal{H}の定義を述べる。

  • ベクトル空間であること
  • 任意の元\mathbf{x}, \mathbf{y} \in \mathcal{H}に対し、内積\langle \mathbf{x}, \mathbf{y} \rangle_{\mathcal{H}}が定義されていること
    • この内積内積の定義さえ満たしていればよい
    • この内積によってノルムが \parallel \mathbf{x} \parallel_{\mathcal{H}} = \sqrt{\langle \mathbf{x}, \mathbf{x} \rangle_{\mathcal{H}}}として定まる
    • さらにこのノルムより2点\mathbf{x}, \mathbf{y}の距離関数がd(\mathbf{x}, \mathbf{y} ) =\parallel\mathbf{x} - \mathbf{y} \parallel_{\mathcal{H}} として定まる
  • 上述の距離関数に関して完備であること
    • すなわち、任意のコーシー列(\parallel \mathbf{x}_m - \mathbf{x}_n \parallel_{\mathcal{H}} \rightarrow 0 \ (m, n \rightarrow \infty))が\mathcal{H}の元に収束すること

なお、完備でない内積空間も完備化という操作によって元の追加と内積の拡張を行うことでヒルベルト空間に変換することができる。

射影定理

ヒルベルト空間には完備性があるおかげて射影定理が成り立ち、これにより通常の線形代数と同じような議論を行うことができる。射影定理とは以下のような定理である。

\mathcal{M}\mathcal{H}の閉部分空間(任意の\alpha, \beta \in \mathbb{R}, \ \  \mathbf{x}_1, \mathbf{x}_2 \in \mathcal{M}に対して\alpha \mathbf{x}_1 + \beta \mathbf{x}_2 \in \mathcal{M}が成り立つ閉集合)、\mathcal{M}^{\perp}\mathcal{M}と直交するベクトル全体とする。このとき、\mathbf{x} \in \mathcal{H}に対し、\mathbf{x} = \mathbf{y} + \mathbf{z} \ \ (\mathbf{y} \in \mathcal{M}, \mathbf{z} \in \mathcal{M}^{\perp}) となる分解が一意に存在する。また、この時\mathbf{y} \mathcal{M}上で最も\mathbf{x}との距離が最小になる唯一の点である*1

この\mathbf{y}\mathbf{x}\mathcal{M}上への直交射影と呼びその作用素P_{\mathcal{M}}と表す。同様に \mathcal{M}^{\perp}上への直交射影作用素P_{\mathcal{M}^{\perp}}と表す。従って、\mathbf{x} = P_{\mathcal{M}}\mathbf{x}+ P_{\mathcal{M}^{\perp}}\mathbf{x} と表すことができる。また、このような分解を直交分解と呼ぶ。

証明は書籍等に譲るが、完備性が射影定理にどう絡んでくるかだけ見ておくと、d = \inf_{\mathbf{y} \in \mathcal{M}} \parallel \mathbf{x} -   \mathbf{y} \parallel_{\mathcal{H}}となる \mathcal{H}の元\mathbf{y} が存在し、d = \min_{\mathbf{y} \in \mathcal{M}} \parallel \mathbf{x} -   \mathbf{y} \parallel_{\mathcal{H}}となることを保証するために必要になる。具体的には、任意のn \in \mathbb{N}に対し、 d^{2} \leq \parallel \mathbf{x} -   \mathbf{y}_n \parallel_{\mathcal{H}}^2 \lt d^{2} + \frac{1}{n^{2}} となる\mathbf{y}_n \in \mathcal{M}が存在し、また点列(\mathbf{y}_n)_{n=1}^\inftyがコーシー列となることから\parallel \mathbf{y}_n -  \mathbf{y}_{*} \parallel = 0 \ \  (n \rightarrow \infty)である\mathbf{y}_{*}  \in \mathcal{M}の存在を保証するために完備性が必要となる。

線形汎関数とリースの表現定理

ヒルベルト空間\mathcal{H}上に定義された写像\varphi: \mathcal{H} \rightarrow \mathbb{R}\varphi(\alpha \mathbf{x} + \beta \mathbf{y}) = \alpha \varphi(\mathbf{x}) + \beta \varphi( \mathbf{y}) \ \ (\alpha, \beta \in \mathbb{R}, \ \  \mathbf{x}, \mathbf{y} \in \mathcal{H})を満たすとき\varphi\mathcal{H}上の線形汎関数と呼ぶ。また、\parallel \varphi \parallel = \sup_{ \parallel \mathbf{x} \parallel_{\mathcal{H}} \leq 1} |\varphi(\mathbf{x})|と定め、この量が有限である場合\varphi有界であるという。加えて、\varphi有界であることと\varphiが連続であることは同値である。

リースの表現定理はこの有界な線形汎関数\varphiについて以下のことを主張する定理である。

\displaystyle 
\varphi(\mathbf{x}) = \langle \mathbf{x}, \mathbf{v} \rangle_{\mathcal{H}}
\tag{6}

となる\mathbf{v} \in \mathcal{H}が一意に存在する。

これまた詳細な証明は書籍に譲り、ここでは存在性についてのみ述べる。存在性の証明には以下の補題が必要になる。

  1. 零空間\mathcal{M} = \ker \varphi=\left\{\mathbf{x} \in \mathcal{H} : \varphi(\mathbf{x}) =0 \right\} は閉部分空間
  2. \dim \mathcal{M}^{\perp} = 1

1.に関しては部分空間であることは\varphiの線形性から明らかで、閉集合であることは\varphiの連続性と\left\{0\right\}\mathbb{R}閉集合であることから逆像\mathcal{M} = \varphi^{-1}(\left\{0\right\})によって示される。

2.は\dim \mathcal{M}^{\perp}  \gt 1である場合は \mathbf{y} \in  \mathcal{M}^{\perp} に対して直交な \mathbf{z} \in  \mathcal{M}^{\perp} を取れるが、\varphi(\mathbf{y}) \neq 0, \varphi(\mathbf{z}) \neq 0であるため適当な定数を選んで\varphi(\mathbf{y}) = \lambda \varphi(\mathbf{z})が成り立つ。線形性より\varphi(\mathbf{y}- \lambda \mathbf{z}) = 0となるため、\mathbf{y}- \lambda \mathbf{z} \in \mathcal{M}である。従って、\langle \mathbf{y}, \mathbf{y}- \lambda \mathbf{z} \rangle_{\mathcal{H}} = 0であり、 \mathbf{z} \mathbf{y}と直交になるように選んでいることから\langle \mathbf{y}, \mathbf{y}- \lambda \mathbf{z} \rangle_{\mathcal{H}} = \langle \mathbf{y}, \mathbf{y} \rangle_{\mathcal{H}} - \langle  \mathbf{y}, \lambda \mathbf{z} \rangle_{\mathcal{H}}  =\parallel \mathbf{y} \parallel_{\mathcal{H}}^2 = 0が得られるため\mathbf{y} = \mathbf{0}である。よって、\mathcal{M}^{\perp}の中で\mathbf{z} \in \mathcal{M}^{\perp}と直交するのは\mathbf{0}のみということは\dim \mathcal{M}^{\perp} = 1だといえる。

さて、この2によって\mathcal{M}^{\perp}上の任意の点は\mathbf{z} \in \mathcal{M}^{\perp} \ \  (\parallel \mathbf{z} \parallel_{\mathcal{H}}=1)を用いて、\alpha \mathbf{z} \ (\alpha \in \mathbb{R})と表現できる。これと1によって任意の\mathbf{x} \in \mathcal{H}\mathbf{x} = P_{\mathcal{M}} \mathbf{x} + \alpha \mathbf{z}と直交分解でき、加えて、

\displaystyle 
\begin{align}
\varphi(\mathbf{x}) &= \varphi(P_{\mathcal{M}} \mathbf{x} + \alpha \mathbf{z}) \\
&= \varphi(P_{\mathcal{M}} \mathbf{x}) + \alpha \varphi(\mathbf{z})\\
&= \alpha \varphi(\mathbf{z}) \\
&= \alpha \varphi(\mathbf{z}) \left(\langle \frac{1}{\alpha} P_{\mathcal{M}} \mathbf{x}, \mathbf{z}  \rangle_{\mathcal{H}} + \langle \mathbf{z}, \mathbf{z} \rangle_{\mathcal{H}} \right) \\
&= \alpha \varphi(\mathbf{z}) \langle \frac{1}{\alpha} P_{\mathcal{M}} \mathbf{x} +  \mathbf{z}, \mathbf{z} \rangle_{\mathcal{H}}\\
&= \langle P_{\mathcal{M}} \mathbf{x} + \alpha  \mathbf{z}, \varphi(\mathbf{z}) \mathbf{z}  \rangle_{\mathcal{H}} \\
&= \langle \mathbf{x} , \varphi(\mathbf{z}) \mathbf{z}   \rangle_{\mathcal{H}} \\
\end{align} 
\tag{7}

と展開できる。そして、有界性より\varphi(\mathbf{z})が有限値になることが保証されているので、式(7)の最後の式において\mathbf{v} =\varphi(\mathbf{z}) \mathbf{z} とすれば式(6)が得られリースの表現定理における存在性が示される。

ヒルベルト空間における超平面と点の距離

ようやく本題だが、ここまでくると最初に述べた実ベクトル空間上の議論とあまり変わらない。以下、前章の結果を用いてヒルベルト空間上の超平面と点の距離の導出を行い、その後に具体的な元と内積が定義されたヒルベルト空間に対し導出した結果を適用した例を示す。

導出

これまでの議論によりヒルベルト空間における超平面とは\mathcal{H}上の有界線形汎関数\varphiにより、

\displaystyle
\varphi(\mathbf{x}) = 0 \ \ (\mathbf{x} \in \mathcal{H})
\tag{8}

となる\mathbf{x}の集合のことだといえる。つまりそれは零空間\mathcal{M} = \ker \varphi であり、それと任意の\mathbf{x}_0 \in \mathcal{H}との距離は

\displaystyle
d = \parallel \mathbf{x}_0 - P_{\mathcal{M}} \mathbf{x}_0\parallel_\mathcal{H}=\parallel P_{\mathcal{M}} \mathbf{x}_0 + P_{\mathcal{M}^{\perp}} \mathbf{x}_0 - P_{\mathcal{M}} \mathbf{x}_0\parallel_\mathcal{H}= \parallel P_{\mathcal{M}^{\perp}} \mathbf{x}_0\parallel_\mathcal{H}
\tag{9}

だといえる。また、前章の議論より\varphi(\mathbf{x}) = \langle \mathbf{x}, \mathbf{v} \rangle_{\mathcal{H}} と書け、P_{\mathcal{M}^{\perp}} \mathbf{x}_0= \alpha \mathbf{v}となる。よって、

\displaystyle
d = \parallel P_{\mathcal{M}^{\perp}} \mathbf{x}_0\parallel_{\mathcal{H}} = | \alpha | \parallel \mathbf{v} \parallel_{\mathcal{H}}
\tag{10}

と式(3)とノルム以外同じ形になる。さらに、

\displaystyle
\varphi(\mathbf{x}_0) = \varphi(P_{\mathcal{M}} \mathbf{x}_0 + \alpha \mathbf{v}) = \varphi(P_{\mathcal{M}} \mathbf{x}_0 )+ \alpha \varphi(\mathbf{v})

\tag{11}

であるが、\mathcal{M}\varphiの零空間であることを思い出すと \varphi(P_{\mathcal{M}} \mathbf{x}_0 ) = 0なので、

\displaystyle
\alpha = \frac{\varphi(\mathbf{x}_0) }{\varphi(\mathbf{v})} =  \frac{\varphi(\mathbf{x}_0) }{\parallel \mathbf{v} \parallel_{\mathcal H}^2}

\tag{12}

となる。これを式(10)に代入することで、

\displaystyle
d =  \frac{|\varphi(\mathbf{x}_0) | }{ \parallel \mathbf{v} \parallel_{\mathcal{H}} }  = \frac{|\langle \mathbf{x}_0, \mathbf{v}\rangle_{\mathcal{H}} |}{\parallel \mathbf{v} \parallel_{\mathcal{H}} } =  \frac{|\langle \mathbf{x}_0, \mathbf{v}\rangle_{\mathcal{H}} |}{\sqrt{\langle \mathbf{v}, \mathbf{v}  \rangle_{\mathcal{H}}} }
\tag{13}

が得られる。

具体例

式(13)で超平面と点の距離を求める例として、最初に述べたn次元実ベクトル空間の例とカーネル関数と呼ばれる関数から導かれるヒルベルト空間の例を挙げる。

n次元実ベクトル空間再考

式(13)を使って、n次元実ベクトルにおける超平面と点の距離を求める。まず、\mathbf{x}, \mathbf{y} \in \mathbb{R}^n内積\langle \mathbf{x}, \mathbf{y} \rangle_{\mathbb{R}^n} = \mathbf{x}^T \mathbf{y}であり、超平面は\varphi(\mathbf{x}) = \langle \mathbf{x}, \mathbf{v} \rangle_{\mathbb{R}^n} = \mathbf{x}^T \mathbf{v}  = \mathbf{v}^T \mathbf{x} = 0なので、そのまま式(13)に当てはめると式(5)が得られる。

カーネル関数によって定義されるヒルベルト空間

まず、2変数x_1, x_2 \in \mathcal{X}上の半正定値対称関数k(x_1, x_2)カーネル関数と呼ぶ。また、カーネル関数の2変数のうち一方を固定した1変数関数をk(\cdot, x) = k_{x}(\cdot) = k_{x}と表記する。今、この1変数関数の有限個の線形和f = \sum_{i} {a_i} k_{x_i} \ (a_i \in \mathbb{R}) が要素である集合\mathcal{V}を考える。k_{x} \mathcal{V}の要素である。そして、この集合\mathcal{V}に対して内積を次のように定義する。

f = \sum_{i} {a_i} k_{x_i}, \ \  g =  \sum_{j} {b_j} k_{y_j} \in \mathcal{V}に対し、

\displaystyle
\langle f, g \rangle_{\mathcal{V}} = \sum_{i}   \sum_{j} {a_i} {b_j} k(y_j, x_i)

\tag{14}

カーネル関数の半正定値対称性から内積の定義を満たす*2。 このようにして定義した内積空間\mathcal{V}を完備化することによりヒルベルト空間\mathcal{H}_kを得る。これでこの節で考える空間を構築できた。なお、\mathcal{V} \subset \mathcal{H}_kであり、f, g \in \mathcal{V}であれば\langle f, g \rangle_{\mathcal{V}} = \langle f, g \rangle_{\mathcal{H}_k}である。

それでは、\varphi(f) = \langle f,   \sum_{i} {c_i} k_{z_i} \rangle_{\mathcal{H}_k} = 0\ \ (c_i \in \mathbb{R}, z_i \in \mathcal{X})となる超平面とf = k_{x} の距離を式(13)を使って計算してみる。そのまま式(13)に代入することで、

\displaystyle
d =  \frac{|\langle k_x, \sum_{i} {c_i} k_{z_i} \rangle_{\mathcal{H}_k} |}{\parallel \sum_{i} {c_i} k_{z_i} \parallel_{\mathcal{H}_k} } =  \frac{ |\sum_{i}    {c_i}  k( z_i, x) |}{\sqrt{\sum_i \sum_j c_i c_j k(z_i, z_j)}}

\tag{15}

が得られる。

もう少しこの空間について詳しく述べると、このように定めた\mathcal{H}_kは再生核ヒルベルト空間と呼ばれ、\mathcal{X}上の関数f \in \mathcal{H}_kx \in \mathcal{X}を代入する操作がf(x) = \langle f, k_{x} \rangle_{\mathcal{H}_k}内積の形で表される「再生性」という性質がある。この性質によってリプレゼンター定理などカーネル法における重要な定理が導かれる。また、再生性は完備化によって得られる元についても成り立つので任意のf \in \mathcal{H}_kと上述の超平面の距離は

\displaystyle
d =  \frac{|\langle f, \sum_{i} {c_i} k_{z_i} \rangle_{\mathcal{H}_k} |}{\parallel \sum_{i} {c_i} k_{z_i} \parallel_{\mathcal{H}_k} } = \frac{|\sum_{i} {c_i}\langle f,  k_{z_i} \rangle_{\mathcal{H}_k} |}{\parallel \sum_{i} {c_i} k_{z_i} \parallel_{\mathcal{H}_k} } =  \frac{ |\sum_{i}    {c_i}  f(z_i) |}{\sqrt{\sum_i \sum_j c_i c_j k(z_i, z_j)}}

\tag{16}

になる。

最後に

今回、抽象的なヒルベルト空間における超平面と点の距離導くため、必要な定義と諸定理について記述した。もちろん今回挙げた例以外にもヒルベルト空間上であれば同じ議論が成り立つ。例えば、閉区間[a, b]上の連続関数の集合C[a, b]についても内積を定義し、それを完備化すればヒルベルト空間になるため今回の議論が成り立つ*3

本記事では超平面と点の距離をトピックとして取り上げたが、それに限らず一度抽象的な世界で議論してその中で成り立つ法則を見つけると、具象の世界ではそれを適用するだけでいいのでかなり有難い。関数解析はその有難みを特に感じられる分野な気がする。

参考文献

*1:この最良近似性は閉部分空間よりも広い集合である閉凸集合に対して成り立つ。

*2:『半』正定値だと内積の定義のうち\langle f, f \rangle_{\mathcal{V}} = 0 \Leftrightarrow f = 0 が成り立たない気もしたが、これは後述の再生性とコーシー・シュワルツの不等式を使って成り立つことが示せる(詳しくは書籍を参照)。

*3:これを具体例の1つとして挙げようかとも思ったがルベーグ積分の知識が必要になりそうで今回は力が及ばなかった。

数理最適化をしっかり学ぶために主双対内点法とそれを使ったSVMを実装

はじめに

シリーズ第3弾で今回は制約付き非線形計画問題を解く主双対内点法を実装した*1。前回の準ニュートン法の記事はこちら。

yamagensakam.hatenablog.com

制約付き非線形計画問題を式で表すと、以下のような最適化問題になる。

\displaystyle 
\begin{align}
\min_{{\bf x}   \in {\mathbb R}^n} \  & f({\bf x}) \\
{\rm s.t. } \  & g_i({\bf x}) \leq 0 \ (i = 1, \cdots, l) \\
 & h_i({\bf x}) = 0 \ (i = 1, \cdots, m) \\
\end{align} \tag{1}

なお今回取り上げる主双対内点法は、すぐ下に記すKKT条件が最適性の必要条件となることが保証されなければ成立しない。これが保証されるために制約条件が満たすべき条件が存在する。そのような制約条件に関する条件を制約想定と言い、1次独立条件やSlater条件などの種類がある。制約想定に関してはこちらが参考になる*2

qiita.com

手法

内点法は最適性の1次の必要条件であるKKT条件を満たす点を実行可能領域の内部を通りながら探す手法である。まずKKT条件について簡単に触れた後、主双対内点法の考え方、アルゴリズムを述べる。

ラグランジュ関数とKKT条件

式(1)の最適化問題に対して、以下のような関数を考える。

\displaystyle 
L({\bf x}, {\bf u}, {\bf v}) = f({\bf x}) + \sum_{i = 1}^{l} u_i g_i({\bf x}) + \sum_{i = 1}^{m} v_i h_i({\bf x}) \tag{2}

これはラグランジュ関数と呼ばれる関数で、式(1)の最適化問題が上述の制約想定を満たす場合、この問題における1次の最適性の必要条件はラグランジュ関数を使って、

\displaystyle
\begin{align}
& \nabla_{{\bf x}} L({\bf x}^{*}, {\bf u}^{*},  {\bf v}^{*})  = \nabla f({\bf x}^{*}) + \sum_{i = 1}^{l} u_i^{*} \nabla g_i({\bf x}^{*}) + \sum_{i = 1}^{m} v_i^{*}  \nabla h_i({\bf x}^{*}) = {\bf 0} \\
& g_i({\bf x}^{*}) \leq 0  , \ \ \   u_i^{*} \geq 0, \ \ \   u_i^{*} g_i({\bf x}^{*}) = 0 \ \ \ (i = 1, \cdots, l) \\
& h_i({\bf x}^{*}) = 0  \ \ \  (i = 1, \cdots, m) \\
\end{align}
\tag{3}

と表すことができる。これはKKT条件と呼ばれ、主双対内点法ではこのKKT条件を満たす({\bf x}^{*}, {\bf u}^{*}, {\bf v}^{*})(に十分に近い点)を探索することになる。

スラック変数とバリア関数の導入

ここで式(1)にスラック変数とバリア関数というものを導入した以下のような最適化問題を考える。

\displaystyle 
\begin{align}
\min_{{\bf x}, {\bf s}} \  & f({\bf x}) - \rho \sum_{i = 1}^{l} \log s_i\\
{\rm s.t. } \  & g_i({\bf x}) + s_i = 0,  \ \ \  s_i \geq 0 \ \ \  (i = 1, \cdots, l) \\
 & h_i({\bf x}) = 0  \ \ \   (i = 1, \cdots, m) \\
\end{align} \tag{4}

s_iがスラック変数でこの変数の導入により不等式制約は等式制約に置き換わる。また、 -  \sum_{i = 1}^{l} \log s_iがバリア関数で実行可能領域の内部では有効な値を取り、\rho > 0であれば制約の境界であるs_i = 0に近づくにつれ無限に発散していく。従って、バリア関数は不等式制約を満たさないことに対するペナルティとして働く。\rho \rightarrow 0で式(4)は元の最適化問題式(1)と等価になる。そのため、主双対内点法では式(4)の停留点方向に移動しながら徐々に\rho0に近づけていくという操作を繰り返すことで式(1)の停留点を見つける。

主双対内点法アルゴリズム

主双対内点法では現時点が({\bf x}^{(k)}, {\bf s}^{(k)}, {\bf u}^{(k)}, {\bf v}^{(k)})だとして、({\bf x}^{(k + 1)}, {\bf s}^{(k + 1)}, {\bf u}^{(k + 1)}, {\bf v}^{(k + 1)}) = ({\bf x}^{(k)} + \alpha {\bf \Delta x}, {\bf s}^{(k)} + \alpha {\bf \Delta s}, {\bf u}^{(k)} + \alpha {\bf \Delta u}, {\bf v}^{(k)} + \alpha {\bf \Delta v})と更新する。以下、更新方向、ステップ幅、\rhoをどのように更新するか述べる。

更新方向の計算

式(4)の停留点{\bar{\bf x}}, {\bar {\bf s}}ではKKT条件より、

\displaystyle
\begin{align}
&  \nabla f({\bar{\bf x}}) + \sum_{i = 1}^{l} \bar{u}_i \nabla g_i({\bar{\bf x}}) + \sum_{i = 1}^{m} {\bar v}_i  \nabla h_i({\bar{\bf x}}) = {\bf 0} \\
& {\bar u}_i {\bar s}_i= \rho, \ \ \   g_i({\bar{\bf x}}) + {\bar s}_i = 0  ,  \ \ \  {\bar u}_i > 0, \ \ \  {\bar s}_i > 0   \ \ \  (i = 1, \cdots, l) \\
& h_i({\bar{\bf x}}) = 0 \ \ \  (i = 1, \cdots, m) \\
\end{align}
\tag{5}

が満たされる。そのため更新方向は式(5)を近似的に満たすように決める。つまり、

\displaystyle
\begin{align}
&  \nabla f({\bf x}^{(k)} + {\bf \Delta x}) + \sum_{i = 1}^{l} {u}_i^{(k)} \nabla g_i({\bf x}^{(k)} + {\bf \Delta x}) + \sum_{i = 1}^{m} v_i^{(k)}  \nabla h_i({\bf x}^{(k)} + {\bf \Delta x}) = {\bf 0} \\
& (u_i^{(k)} + {\Delta u}_i) (s_i^{(k)} + {\Delta s}_i) = \rho_k, \ \ \   g_i({\bf x}^{(k)} + {\bf \Delta x}) + (s_i^{(k)} + {\Delta s}_i)  = 0 \\
& (u_i^{(k)} + {\Delta u}_i)  > 0, \ \ \   (s_i^{(k)} + {\Delta s}_i)   > 0   \ \ \ (i = 1, \cdots, l) \\
& h_i({\bf x}^{(k)} + {\bf \Delta x}) = 0 \ \ \ (i = 1, \cdots, m) \\
\end{align}
\tag{6}

を満たすように({\bf \Delta x}, {\bf \Delta s}, {\bf \Delta u}, {\bf \Delta v})を決めたいのだが、式(6)の方程式を解くのは困難なので1次近似して解く。1次近似なので、

\displaystyle
\begin{align}
\nabla f({\bf x}^{(k)} + {\bf \Delta x}) & \simeq \nabla f({\bf x}^{(k)} ) +\nabla^2 f({\bf x}^{(k)})  {\bf \Delta x} \\
\nabla g_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq \nabla g_i({\bf x}^{(k)} ) +\nabla^2 g_i({\bf x}^{(k)})  {\bf \Delta x} \\
\nabla h_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq \nabla h_i({\bf x}^{(k)} ) +\nabla^2 h_i({\bf x}^{(k)})  {\bf \Delta x} \\
g_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq g_i({\bf x}^{(k)} ) +(\nabla g_i({\bf x}^{(k)}))^T  {\bf \Delta x} \\
h_i({\bf x}^{(k)} + {\bf \Delta x}) & \simeq h_i({\bf x}^{(k)} ) +(\nabla h_i({\bf x}^{(k)}))^T  {\bf \Delta x} 
\end{align}
\tag{7}

を式(6)の等式に代入し、2次項\Delta u_i \Delta s_iを無視すると、以下の1次連立方程式が得られる。

\displaystyle
\left(
    \begin{array}{rrrr}
      \nabla_{{\bf x}{\bf x}}^2 L & {\bf O} & {\bf J}_g^T & {\bf J}_h^T \\
       {\bf O} & {\bf D}_u & {\bf D}_s &  {\bf O} \\
       {\bf J}_g &  {\bf I}  &  {\bf O}  &  {\bf O}  \\
       {\bf J}_h  &  {\bf O}  &  {\bf O}  &  {\bf O}
    \end{array}
\right)
\left(
    \begin{array}{r}
      {\bf \Delta x} \\
      {\bf \Delta s} \\ 
      {\bf \Delta u} \\ 
      {\bf \Delta v} \\ 
    \end{array}
\right)
=
\left(
    \begin{array}{r}
      -{\nabla}_{\bf x} L \\
      {\bf \rho_k {\bf 1} - {\bf s} \odot {\bf u}} \\ 
      -({\bf g} + {\bf s}) \\ 
      -{\bf h} \\ 
    \end{array}
\right)
\tag{8}

ただし、

\displaystyle
\begin{align}
&\nabla_{{\bf x}{\bf x}}^2 L =  \nabla^2 f({\bf x}^{(k)}) + \sum_{i = 1}^{l} {u}_i^{(k)} \nabla^2 g_i({\bf x}^{(k)}) + \sum_{i = 1}^{m} {v}_i^{(k)}  \nabla^2 h_i({\bf x}^{(k)}), \\
&{\bf J}_g =\left( \nabla  g_1({\bf x}^{(k)}), \cdots,  \nabla g_l({\bf x}^{(k)}) \right)^T, \\
&{\bf J}_h =\left( \nabla  h_1({\bf x}^{(k)}), \cdots,  \nabla  h_m({\bf x}^{(k)}) \right)^T, \\
&{\bf D}_u = {\rm diag} ({\bf u^{(k)}}), \\
&{\bf D}_s = {\rm diag} ({\bf s^{(k)}}), \\
&{\nabla}_{\bf x} L =\nabla f({\bf x}^{(k)}) + \sum_{i = 1}^{l} {u}_i^{(k)} \nabla g_i({\bf x}^{(k)}) + \sum_{i = 1}^{m} {v}_i^{(k)}  \nabla h_i({\bf x}^{(k)}) ,\\
&{\bf g} =\left( g_1({\bf x}^{(k)}), \cdots,  g_l({\bf x}^{(k)}) \right)^T, \\
&{\bf h} =\left( h_1({\bf x}^{(k)}), \cdots,  h_m({\bf x}^{(k)}) \right)^T ,\\
&{\bf 1} = \left( 1, \cdots, 1 \right)^T
\end{align}
\tag{9}

であり、\odotアダマール積を表す。この式(8)を解くことで更新方向が得られる。

ステップ幅の決め方

次にステップ幅\alphaだが、まず{\bf s}^{(k)} + \alpha {\bf \Delta s} > 0,  {\bf u}^{(k)} + \alpha {\bf \Delta u} > 0となる必要があるため、これを満たす 0 \lt \alpha \leq 1を見つける*3。更にその条件を満たす中でもなるべく適切な\alphaを選びたいので、以下の目的関数と制約のペナルティを組み合わせたメリット関数を定義し、

\displaystyle
\phi_{(\eta, \rho)}({\bf x}, {\bf s}) = f({\bf x}) - \rho \sum_{i = 1}^{l} \log s_i+ \eta \sum_{i = 1}^{l} |g_i({\bf x}) + s_i| + \eta \sum_{i = 1}^{m} |h_i({\bf x})| \tag{10}

\phi_{(\eta, \rho_k)}({\bf x}^{(k)} + \alpha {\bf \Delta x} ,  {\bf s}^{(k)} + \alpha {\bf \Delta s} ) を最小にする\alphaを直線探索により見つける*4\eta > 0はペナルティの重みを表すパラメータで、今回の実装では固定値をあらかじめ与えている。

\rhoの更新

最後に\rhoの更新だが、これは{\bf x}^{(k + 1)}, {\bf s}^{(k + 1)},  {\bf u}^{(k + 1)},  {\bf v}^{(k + 1)} が式(6)を満たすのならs_i^{(k + 1)} u_i^{(k + 1)} = \rho_kなので*5 \rho_k = \frac{({{\bf s}^{(k + 1)}})^{T}  {\bf u}^{(k + 1)}}{l}となっているのが理想的である。このことから、適切なパラメータ0 \lt \delta \lt 1をあらかじめ与えて、

\displaystyle
\rho_{k + 1} = \delta \frac{({{\bf s}^{(k + 1)}})^{T}  {\bf u}^{(k + 1)}}{l}
\tag{11}

として更新する。

アルゴリズムまとめ

主双対内点法アルゴリズムは以下の通り。

  1. 初期点{\bf x}^{(0)}, {\bf s}^{(0)},  {\bf u}^{(0)},  {\bf v}^{(0)} とバリア関数のパラメータの初期値\rho_0を定める。
  2. \rho_kが十分に小さければ終了。
  3. 式(8)を解いて探索方向 ( {\bf \Delta x}, {\bf \Delta s}, {\bf \Delta u}, {\bf \Delta v})を決める。
  4. ステップ幅\alphaを決める。
  5. ({\bf x}^{(k + 1)}, {\bf s}^{(k + 1)}, {\bf u}^{(k + 1)}, {\bf v}^{(k + 1)}) = ({\bf x}^{(k)} + \alpha {\bf \Delta x}, {\bf s}^{(k)} + \alpha {\bf \Delta s}, {\bf u}^{(k)} + \alpha {\bf \Delta u}, {\bf v}^{(k)} + \alpha {\bf \Delta v})として更新する。
  6. 式(11)で\rho_{k+1}を求める。
  7. 2に戻る。

その他補足

実装する上で困ったこと、まだよく分かってないこと、その他細かい話など。

初期点について

内点法は基本的に実行可能領域の内部を通って解を見つけるアルゴリズムなので初期点も内部に置かないといけないと思ったのだが、実は今回述べたアルゴリズムだと初期点を実行可能領域外の点にしても試した範囲では実行可能領域の停留点に収束した。

今回の方法が実行可能領域外を初期点にしても収束するという理論的保証があるかはよく分かっていないが、雰囲気的には\rhoが十分に大きければ式(4)の最適解は実行可能領域内の境界から遠く離れた点であり、式(8)を満たす方向に点を動かすということはその式(4)の最適解に近づけることを意味するので、たとえ初期点が実行可能領域外だったとしても内部に入る方向に動き最終的には内部の停留点に収束するという感じだろう*6。ということで実装では、適当な値を初期値としている。

非凸な最適化問題の場合

ラグランジュ関数のヘッセ行列\nabla_{{\bf x}{\bf x}}^2 L(式(9)参照)が正定値でない場合、式(8)を満たす{\bf \Delta x} が下降方向とは限らない。一方、ステップ幅を探索するときは、なるべくメリット関数が小さくなる\alphaを探すので、{\bf \Delta x} が上昇方向であれば限りなく0に近い\alphaが得られなかなか収束しないということが起こりうる。そのため、今回の実装では上述の通り\alphaが小さくなりすぎないよう探索の範囲を絞っている。また、そのほかの対策方法として、ヘッセ行列の代わりに準ニュートン法同様にヘッセ行列を近似する正値対称行列を用いる場合もある模様。

実装

以下が実装した主双対内点法pythonコード。

import numpy as np
import copy

# 以下の最適化問題を定義するクラス
# min f(x) : objective
# s.t. g_i(x) <= 0 (i = 1..l) : inequality_constraints_list
#      h_i(x) = 0  (i = 1..m) : equality_constraints_list
# objective, inequality_constraints_list[i], equality_constraints_list[i]は、
# それぞれ以下のメンバを持ってる必要あり
#      func(x)   : return f(x)
#      grad(x)   : return ∇f(x)
#      hessian(x): return ∇^2 f(x)
#      dim       : 変数の次元
class problem():
    def __init__(self, objective, 
                inequality_constraints_list,
                equality_constraints_list):

        for constraint in inequality_constraints_list:
            if objective.dim != constraint.dim:
                raise Exception('dimension mismatch')
        
        for constraint in equality_constraints_list:
            if objective.dim != constraint.dim:
                raise Exception('dimension mismatch')
        
        self.dim = objective.dim
        self.objective = objective
        self.inequality_constraints_list = inequality_constraints_list
        self.equality_constraints_list = equality_constraints_list

def merit_function(problem, rho, eta, x, s):
    res = problem.objective.func(x)
    res -= rho * np.sum(np.log(s))
    g_sum = 0
    for i, constraint in enumerate(problem.inequality_constraints_list):
        g_sum += np.abs(constraint.func(x) + s[i])
    res += eta * g_sum
    h_sum = 0
    for i, constraint in enumerate(problem.equality_constraints_list):
        h_sum += np.abs(constraint.func(x))
    res += eta * h_sum 
    return res 

def backtracking_line_search(problem, x, s, u, dx, ds, du, rho, eta, beta):
    phi = lambda x, s:merit_function(problem, rho, eta, x, s)

    alpha1 = 1.0
    alpha2 = 1.0

    # αの初期値はs + α ds > 0 and u + α du > 0を満たす最大値
    # ただし、α <= 1
    if(np.min(ds) < 0):
        alpha1 = np.min(-s[ds < 0]/ds[ds < 0])
    if(np.min(du) < 0):
        alpha2 = np.min(-u[du < 0]/du[du < 0])
    alpha = min(1.0, alpha1 * beta, alpha2 * beta)

    min_merit = np.inf
    res_alpha = alpha
    init_alpha = alpha

    # 探索範囲はαの初期値の1/10
    while(alpha > init_alpha * 0.1):
        new_x = x + alpha * dx
        new_s = s + alpha * ds

        m = phi(new_x, new_s)
        if (m < min_merit):
            min_merit = m
            res_alpha = alpha
        alpha *= beta
    return res_alpha

# problem : 最適化問題
# eps     : 収束判定基準
# eta     : メリット関数パラメータ
# beta    : 直線探索の減衰パラメータ
# t       : ρ更新時のパラメータ
def interior_point(problem, eps = 1e-9, eta =0.1, beta = 0.9, t = 0.5):
    
    l = len(problem.inequality_constraints_list) #不等式制約数
    m = len(problem.equality_constraints_list) #等式制約数
    n = problem.dim #変数の次元

    rho = 1
    # 初期値xはとりあえず制約を満たさなくてもよい
    x = np.ones(n)
    # s の初期値は s > 0
    s = np.ones(l)
    # s_i u_i = ρを満たすようにuを設定
    u = rho / s

    v = np.zeros(m)
    
    while(1):
        if rho < eps:
            break
        H = problem.objective.hessian(x)
        d = problem.objective.grad(x)
        Jg = np.zeros((l, n))
        g = np.zeros(l)
        # 不等式制約の関数がらみの計算
        for i, constraint in enumerate(problem.inequality_constraints_list):
            H += constraint.hessian(x) * u[i]
            Jg[i, :] = constraint.grad(x)
            d += Jg[i, :] * u[i]
            g[i] = constraint.func(x)
        
        Jh = np.zeros((m, n))
        h = np.zeros(m)
        # 等式制約の関数がらみの計算
        for i, constraint in enumerate(problem.equality_constraints_list):
            H += constraint.hessian(x) * v[i]
            Jh[i, :] = constraint.grad(x)
            d += Jh[i, :] * v[i]
            h[i] = constraint.func(x)

        Du = np.diag(u)
        Ds = np.diag(s)

        # make array P       
        # |H   O   Jg.T  Jh.T|
        # |O   Du  Ds    O   |
        # |Jg  I   O     O   |
        # |Jh  O   O     O   | 
        P = np.zeros((n + 2 * l + m, n + 2 * l + m))
        P[:n, :n] = H
        P[:n, n + l: n + 2 * l] = Jg.T
        P[:n, n + 2 * l:] = Jh.T
        P[n:n + l, n:n + l] = Du
        P[n:n + l, n + l:n + 2 * l] = Ds
        P[n + l:n + 2 * l, :n] = Jg
        P[n + l:n + 2 * l, n:n + l] = np.eye(l)
        P[n + 2 * l:, :n] = Jh

        # make vector r
        r = np.concatenate([-d, rho - s * u, -(g + s), -h])

        delta = np.linalg.solve(P, r)

        dx = delta[:n]
        ds = delta[n:n + l]
        du = delta[n + l:n + 2 * l]
        dv = delta[n + 2 * l:]

        alpha = backtracking_line_search(problem, x, s, u,
                                            dx, ds, du, rho, eta, beta)

        # x,s,u,vの更新
        x += alpha * dx
        s += alpha * ds
        u += alpha * du
        v += alpha * dv

        # ρの更新
        rho =  t * u @ s / l

    return x, s, u, v

式(1)の最適化問題problemというクラスに記述し、そのインスタンスを関数interior_pointに渡すと停留点を探して出力する。problemobjectiveに目的関数を表すクラスのインスタンスをセットし、inequality_constraints_listには不等式制約の関数(式(1)のg_i({\bf x}) )を表すクラスのインスタンスのリスト、inequality_constraints_listには等式制約の関数(式(1)のh_i({\bf x}) )を表すクラスのインスタンスのリストをそれぞれセットする。各クラスはf({\bf x}), g_i({\bf x}), h_i({\bf x})の値を返すメソッドfunc\nabla f({\bf x}), \nabla g_i({\bf x}),\nabla  h_i({\bf x})を返すメソッドgrad\nabla^{2} f({\bf x}), \nabla^{2} g_i({\bf x}),\nabla^{2}  h_i({\bf x})を返すメソッドhessianを実装してる必要がある。

なお、式(9)の連立方程式np.linalg.solveで解いているが、零要素が多いためもう少し効率的に解く方法があるかもしれない。また、この実装だと不等式制約が存在しない問題は解けない(等式制約が存在しない問題は解ける)。

ついでにSVM最適化問題内点法で解く

SVMは分類問題を解く教師あり学習の一種であり、その「学習」は最適化問題を解くことに対応する。通常はSMOなどの専用の効率が良い最適化アルゴリズムなどで解くが、今回はこれを主双対内点法で解いてみる。

SVMで解く最適化問題

SVMの学習で解くことになる最適化問題を簡単に説明する。SVMの学習では教師データ集合\left\{({\bf x}_i, y_i) \right\}_{i=1}^N, ({\bf x}_i \in \mathbb{R}^n, y_i \in \left\{-1, 1 \right\})が得られてるとして、それをうまいこと分離できる関数fとバイアス項bを求める。具体的にはヒンジ損失関数に正則化項を加えた以下の最小化問題を解く。

 \displaystyle
\min_{f, b} \left\{ \frac{1}{2} \parallel f \parallel_{\mathcal H}^2+ C \sum_{i = 1}^{N} \max(1 - y_i (f({\bf x}_i) + b), 0) \right\}
\tag{12}

この問題は関数のfに関する最適化問題であり、加えて第2項(ヒンジ損失関数)は微分できないので、このまま解くのは難しい。そこでリプレゼンター定理などを使ってゴニョゴニョしたりヒンジ損失関数の制約条件への変換を行ったりすると、以下の凸2次計画問題に変形できる。

\displaystyle
\begin{align}
\min_{{\bf a}, b, {\bf \xi}} \  & \frac{1}{2}{\bf a}^T {\bf K} {\bf a} + C \sum_{i=1}^{N} \xi_i  \\
{\rm s.t.} \ & -  y_i \left( \sum_{j = 1}^N a_j k({\bf x}_i,{\bf x}_j)+ b \right) + 1 - \xi_i \leq 0, \ \xi_i  \geq 0 \ (i = 1, \cdots, N)
\end{align}
\tag{13}

ここでk({\bf x}, {\bf x}')カーネル関数と呼ばれ適切な正定値対称関数を選ぶ(例えば、 k({\bf x}, {\bf x}') = \exp(-\gamma \parallel {\bf x} - {\bf x}' \parallel^{2}) など)。また{\bf K} \in{\mathbb R}^{N \times N}ij列目にk(x_i, x_j)を持つグラム行列である。 これらの理屈などについては、以下の資料や参考文献として挙げた『カーネル多変量解析』という書籍を参考にするとよい*7

https://www.ism.ac.jp/~fukumizu/ISM_lecture_2004/Lecture2004_kernel_method.pdf

この最適化問題式(13)の解{\bf a}^{*}b^{*}を用いて、新たに与えられた入力{\hat{\bf x}}

\displaystyle
{\rm sign} \left( \sum_{i = 1}^N a_i^{*} k({\bf x}_i, {\hat {\bf x}})+ b^{*} \right)
\tag{14}

としてクラスを分類できる。

SVMの実装

以下が実装したSVMとそれを用いて検証用データの学習と分類を行うpythonコード。上述の主双対内点法の関数を利用している。

import numpy as np
import interior_point_method as ipm
import matplotlib.pyplot as plt

# 目的関数を表すクラス
class objective_function:
    def __init__(self, K, C, var_dim, constraint_num):
        self.K = K
        self.C = C
        self.a_dim = var_dim
        self.z_dim = constraint_num # ξはzで表現
        self.dim = var_dim + constraint_num + 1 # aの次元 + ξの次元 + bの次元

    # a_z_bはaとξとbを1つに並べたベクトル(array)
    def func(self, a_z_b):
        return 0.5 * a_z_b[:self.a_dim] @ self.K @ a_z_b[:self.a_dim] + \
                self.C * np.sum(a_z_b[self.a_dim:self.a_dim + self.z_dim])
           
    def grad(self, a_z_b):
        return np.concatenate([self.K @ a_z_b[:self.a_dim],
                               self.C * np.ones(self.z_dim),
                               [0]])

    def hessian(self, a_z_b):
        H = np.zeros((self.dim, self.dim))
        H[:self.a_dim, :self.a_dim] = self.K
        return H

# 教師データによる不等式制約クラスを生成するクラス
class train_data_constraint_creator:
    def __init__(self, K, y):
        self.K = K
        self.y = y

    # 教師データによる不等式制約を表すクラス
    class train_data_constraint:
        def __init__(self, K, y, idx):
            self.idx = idx
            self.K = K
            self.y = y
            self.a_dim = K.shape[0]
            self.z_dim = y.shape[0]
            self.dim = self.a_dim + self.z_dim + 1
        
        def func(self, a_z_b):
            return -a_z_b[self.a_dim + self.idx] + 1 - \
                    self.y[self.idx] * \
                    (a_z_b[:self.a_dim] @ self.K[self.idx, :] + a_z_b[-1])
             
        def grad(self, a_z_b):
            return np.concatenate([-self.K[self.idx, :] * self.y[self.idx],
                                    np.ones(self.z_dim) * -1,
                                    [-self.y[self.idx]]])

        def hessian(self, a_z_b):
            return np.zeros((self.dim, self.dim))

    def create(self, idx):
        return self.train_data_constraint(self.K, self.y, idx)

# 非負制約クラスを生成するクラス
class nonnegative_constraint_creator:
    def __init__(self, K, y):
        self.constraint_num = y.shape[0]
        self.var_dim = K.shape[0]
        self.constraint_num = y.shape[0]
        self.dim = self.var_dim + self.constraint_num + 1

    # 非負制約を表すクラス
    class function_nonnegative_constraint:
        def __init__(self, var_dim, constraint_num, idx):
            self.idx = idx
            self.a_dim = var_dim
            self.z_dim = constraint_num
            self.dim = self.a_dim + self.z_dim + 1
        
        def func(self, a_z_b):
            return -a_z_b[self.a_dim + self.idx]
        
        def grad(self, a_z_b):
            res = np.zeros(self.dim)
            res[self.a_dim + self.idx] = -1
            return res

        def hessian(self, a_z_b):
            return np.zeros((self.dim, self.dim))
             
    def create(self, idx):
        return self.function_nonnegative_constraint(self.var_dim,
                                                self.constraint_num, idx)

class svm:
    def __init__(self, C, kernel_func):
        self.C = C
        self.K = None
        self.X = None
        self.kernel_func = kernel_func

    def fit(self, X, y):
        self.X = X
        self.K = np.zeros((X.shape[0], X.shape[0]))

        # グラム行列の計算
        for i in range(X.shape[0]):
            for j in range(i, X.shape[0]):
                self.K[i, j] = self.kernel_func(X[i, :], X[j, :])
                self.K[j, i] = self.K[i, j]

        # 各種不等式制約インスタンス生成用
        tdcc = train_data_constraint_creator(self.K, y)
        ncc = nonnegative_constraint_creator(self.K, y)
        # 最適化問題のインスタンスを生成
        Pro = ipm.problem(objective_function(self.K, self.C,
                                            self.K.shape[1], y.shape[0]), 
                        [tdcc.create(i) for i in range(y.shape[0])] + 
                        [ncc.create(i) for i in range(y.shape[0])],
                        [])
        
        a_z_b_opt, s_opt, u_opt, v_opt = ipm.interior_point(Pro, eps = 1e-15)
        self.a = a_z_b_opt[:self.K.shape[0]]
        self.b = a_z_b_opt[-1]
        self.u = u_opt
        self.support_vectors = np.where(self.u[:X.shape[0]] > 1e-8)

    def predict(self, X_):
        if X_ is None:
             raise Exception('not learned')
            
        K_ = np.zeros((self.K.shape[0], X_.shape[0]))
        for i in range(self.K.shape[0]):
            for j in range(X_.shape[0]):
                K_[i, j] = self.kernel_func(self.X[i], X_[j])
        return np.sign(self.a @ K_ + self.b)


# テストデータ生成用
def make_data(xmin, xmax, n):
    np.random.seed(1111)
    cnt = 0
    t = 0
    N = n * 9
    X_train = np.zeros((N, 2))
    y_train = np.zeros(N)
    center1 =  np.linspace(-(-xmin[0] + xmax[0])/3.5, 
                            (-xmin[0] + xmax[0])/3.5, 3)
    center2 =  np.linspace(-(-xmin[1] + xmax[1])/3.5, 
                            (-xmin[1] + xmax[1])/3.5, 3)
    for i in range(3):
        for j in range(3):
            y_train[t:t + n] = 1 if cnt % 2 == 0 else -1
            X_train[t:t + n, :] = 1.5 * np.random.randn(n, 2) + \
                                np.array([center1[i], center2[j]])
            t += n
            cnt += 1

    #テスト入力点生成
    T1 = np.linspace(xmin[0], xmax[0], 60)
    T2 = np.linspace(xmin[1], xmax[1], 60)
    X_test = np.zeros((T1.shape[0] * T2.shape[0], 2))
    i = 0
    for t1 in T1:
        for t2 in T2:
            X_test[i, :] = np.array([t1, t2])
            i += 1
    return X_train, y_train, X_test

def main():
    #教師データとテストデータの生成
    xmin = np.array([-9, -9])
    xmax = np.array([9, 9])
    n = 20
    X_train, y_train, X_test = make_data(xmin, xmax, n)

    # カーネル関数(ガウスカーネル)
    gamma = 0.5
    kernel_func = lambda x_i, x_j : np.exp(-gamma * (x_i - x_j) @ (x_i - x_j))

    sv = svm(0.1, kernel_func)
    # SVM学習
    sv.fit(X_train, y_train)
    # テストデータの予測
    y_pred = sv.predict(X_test)
    y_pred[y_pred == -1] = 0

    # 以下結果のプロット
    size = 100
    sc = plt.scatter(X_test[:, 0], X_test[:, 1],
                        size, y_pred, marker = ',', cmap='bwr')

    plt.scatter(X_train[y_train == 1, 0], X_train[y_train == 1, 1] ,
                edgecolors="black", marker="o" , color='r',
                label= "train data (y = 1)")
    plt.scatter(X_train[y_train == -1, 0], X_train[y_train == -1, 1] ,
                edgecolors="black", marker="o" , color='b',
                label= "train data (y = -1)")
    plt.axis([xmin[0], xmax[0], xmin[1], xmax[1]])
    plt.legend(loc='upper left')
    plt.show()

if __name__ == '__main__':
    main()

先述のproblemクラスとinterior_point関数を使ってSVMを実装している。基本的にはこれまで述べてきたことをそのまま書いてるだけなので特筆すべきことはないと思う。これを実行すると以下の結果が得られる。

f:id:YamagenSakam:20210124133035p:plain

ガウスカーネルを使っているので非線形な分離できている。 ちなみに、グラム行列{\bf K}が正則でない場合(例えば、\gammaを小さくしすぎた場合や全く同じ学習サンプルが2つ以上ある場合などにあり得る)、主双対内点法のコードにおけるPも正則でなくなり計算が破綻する。このあたりの弱点はSMO法などで克服される。

ここに挙げた主双対内点法のコードやSVMのコードに加えて、線形のSVMやその他の最適化問題で確認したコードもあるので、それらは全て以下のgithubリポジトリに上げている。

github.com

最後に

今回は制約付き非線形計画問題を解く主双対内点法を実装した。制約付き非線形計画問題を解く手法としては他にも逐次2次凸計画法や乗数法(拡張ラグランジュ関数法)などがある。ちなみに乗数法については、以前その改良の1つである交互方向乗数法(ADMM)を実装した記事を書いている。

yamagensakam.hatenablog.com

ということで、非線形計画問題については一旦ここまでにして次から整数計画問題に移る予定。

参考文献

*1:主双対内点法は数ある内点法の一種だが、単に内点法と書いて主双対内点法のことを指してる文献もある。

*2:制約想定以外にも主双対内点法について色々書かれているので、今書いてる自分の記事は不要じゃないかとも思ったが...まぁいいか。

*3:{\bf s}^{(k)} > 0, {\bf u}^{(k)} > 0なので、条件を満たす\alphaは必ず存在する。

*4:実装では条件を満たす中でもなるべく大きい\alphaを初期値としてバックトラッキングを行っている。ただし、探索範囲の下限を小さくすると、極端に小さい\alphaが得られることがあるので探索範囲の下限は初期値の1/10程度までにとどめている。

*5:1次近似して式(6)を解いているので一般には満たされない。

*6:実行可能とは限らない初期点を使うことのできる内点法を『インフィージブル内点法』というらしい。今回がその一手法に当たるかは不明。

*7:かなり説明を端折ったのでもう少し補足すると、SVMの学習をソフトマージンのマージン最大化問題として定式化→それは式(12)の正則化付きヒンジ損失最小化と等価→これにリプレゼンター定理を適用して有限次元の最適化問題に帰着→これをソフトマージンの問題に戻すと有限次元の2次計画問題である式(13)が出てくるという流れ。詳しくは貼った資料のP39~P45参照。

数理最適化をしっかり学ぶために準ニュートン法を実装

はじめに

以前、以下の単体法の記事を書いた。

yamagensakam.hatenablog.com

今回はこのシリーズ第2弾で、以下の制約なし非線形計画問題を準ニュートン法で解くプログラムを実装する。

\displaystyle \min_{{\bf x} \in \mathbb{R}^n} f({\bf x}) \tag{1}

手法

制約なし非線形計画問題の解き方

非線形計画問題を解く際よく行われるのは、最適性の1次の必要条件である\nabla f({\bf x}_{*}) = {\bf 0}を満たす {\bf x}_{*}(に十分に近い点)を以下の反復計算で見つける方法である。

\displaystyle {\bf x}^{(k + 1)} = {\bf x}^{(k)} + \alpha_{k} {\bf d}({\bf x}^{(k)}) \tag{2}

ここで{\bf d}({\bf x}^{(k)})は探索方向で、f({\bf x}^{(k)} +  \alpha_{k}  {\bf d}({\bf x}^{(k)}))が減少する方向に探索方向を選ぶ。探索方向の選び方には様々なバリエーションがあり、その1つが今回実装した準ニュートン法。なお\alpha_{k}はステップ幅で、これについても目的関数が減少するよう適切に選ぶ必要があるが、こちらの選び方の説明は今回省略する*1

ニュートン法とその問題点

準ニュートン法の前にニュートン法について書く。f({\bf x})のヘッセ行列を\nabla^{2} f({\bf x})とするとニュートン法では探索方向を

\displaystyle {\bf d}({\bf x}^{(k)}) = - {\nabla^{2} f({\bf x}^{(k)})}^{-1} \nabla f({\bf x}^{(k)}) \tag{3}

とする。これは、目的関数を{\bf x}^{(k)}まわりでの2次近似した関数

\displaystyle f({\bf x}^{(k)} + {\bf d}) \sim f({\bf x}^{(k)})  + \nabla f({\bf x}^{(k)})^T {\bf d} + \frac{1}{2} {\bf d}^T \nabla^2 f({\bf x}^{(k)}) {\bf d} \tag{4}

{\bf d}に関する停留点に当たり、ヘッセ行列が正定値なら2次近似の最小解となる。従って、ニュートン法は目的関数を2次近似してその停留点を求めるという操作を収束するまで繰り返すアルゴリズムだといえる。ニュートン法は2次収束するため収束性はよい。

一方問題点としては、目的関数が狭義凸関数でない場合ヘッセ行列が正定値にならないケースがある*2。そのため、{\bf d}が下降方向になることは保証されない。そもそもヘッセ行列が逆行列を持たない可能性もある。更にはヘッセ行列を計算すること自体が難しい場合もある。

準ニュートン法アルゴリズム

ニュートン法の上記の問題を解消するため、ヘッセ行列そのものではなくヘッセ行列を近似する正定値対称行列{\bf B}_kを使って探索方向を

\displaystyle {\bf d}({\bf x}^{(k)}) = - {{\bf B}_{k}}^{-1} \nabla f({\bf x}^{(k)}) \tag{5}

とするのが準ニュートン法{\bf B}_kが正定値なので探索方向は常に下降方向であり逆行列も存在する。

ヘッセ行列の近似

直前の反復でヘッセ行列を近似した正定値対称行列{\bf B}_kが得られているとして、これを用いて次の近似行列{\bf B}_{k+1}を求めることを考える。まず、\nabla f({\bf x}^{(k)}) の1次近似として、

\displaystyle \nabla f({\bf x}^{(k)}) = \nabla f({\bf x}^{(k + 1)} + ({\bf x}^{(k)}  - {\bf x}^{(k + 1)} )) \sim   \nabla f({\bf x}^{(k + 1)}) + \nabla^2 f({\bf x}^{(k + 1)}) ({\bf x}^{(k)}  - {\bf x}^{(k + 1)} ) \tag{6}

が得られる。式(6)の関係よりヘッセ行列\nabla^{2} f({\bf x}^{(k + 1)})の近似行列{\bf B}_{k+1}としては以下の条件を満たすのが適切だといえる。

\displaystyle {\bf B}_{k + 1}  ({\bf x}^{(k + 1)}  - {\bf x}^{(k)}  )  =\nabla f({\bf x}^{(k + 1)}) - \nabla f({\bf x}^{(k)})\tag{7}

この式(7)はセカント条件と呼ばれる。これに加えて {\bf B}_{k + 1}  は正定値対称行列でなければならない。これらの条件を満たす行列はいくつか提案されているが、中でもよく知られているのはBFGS公式と呼ばれるもので、これは直前の反復で得られた正定値対称行列 {\bf B}_{k} および{\bf s}^{(k)} ={\bf x}^{(k + 1)}  - {\bf x}^{(k)} , \ {\bf y}^{(k)} = \nabla f({\bf x}^{(k + 1)}) - \nabla f({\bf x}^{(k)}) を用いて

\displaystyle {\bf B}_{k + 1} = {\bf B}_{k} - \frac{{\bf B}_{k} {\bf s}^{(k)} ({\bf B}_{k} {\bf s}^{(k)})^T}{({\bf s}^{(k)})^T  {\bf B}_{k} {\bf s}^{(k)}} + \frac{{\bf y}^{(k)} ({\bf y}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} }\tag{8}

と計算できる。詳細は省くが式(8)の両辺に右から{\bf s}^{(k)}をかければ {\bf B}_{k + 1} {\bf s}^{(k)} = {\bf y}^{(k)}となりセカント条件を満たすことを示せる。また、ステップ幅\alphaを適切に決めれば正定値性は満たされる*3

逆行列の更新

式(5)にあるように探索方向を求めるために {\bf B}_{k} 逆行列が必要になるが、次元が大きくなると更新毎に式(8)で {\bf B}_{k+1} を計算してその逆行列を計算する、なんてやってると計算量的に厳しい。そこで式(8)の右辺第2項、第3項が共にランク1の行列であることを考えると、Sherman-Morrisonの公式*4を2回適用することで逆行列そのものの更新式を導くことができる。具体的な式展開などは以下の記事が参考になる。

qiita.com

結果としては、

\displaystyle ({\bf B}_{k + 1})^{-1} = \left({\bf I} - \frac{{\bf s}^{(k)} ({\bf y}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} } \right)({\bf B}_{k})^{-1} \left({\bf I} - \frac{{\bf y}^{(k)} ({\bf s}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} } \right) + \frac{{\bf s}^{(k)} ({\bf s}^{(k)})^T }{({\bf s}^{(k)})^T {\bf y}^{(k)} }  \tag{9}

という更新式が得られる。

BFGS公式について

式(8)でやや唐突にBFGS公式を書いたが、ここではBFGS公式の考え方を示す。まずBFGS公式含めて提案されている多くの更新式は、正定値対称性とセカント条件を満たす行列のうち、直前の行列{\bf B}_{k}にもっとも近い行列を{\bf B}_{k+1}としている。これはヘッセ行列は急激に変化しないという事前知識を基づく。これを式で書くと適当な行列間の非類似度を測る関数D(\cdot, \cdot)を導入して、

\displaystyle
\begin{aligned}
\min_{{\bf B}} \ &  D({\bf B}, {\bf B}_{k}) \\
{\rm s.t.} \ & {\bf B} {\bf s}^{(k)} = {\bf y}^{(k)}  \\
\ & {\bf B} \succ {\bf O} \\
\ & {\bf B} = {\bf B}^T
\end{aligned}
\tag{10}

という最適化問題の解を{\bf B}_{k+1}としている。

BFGS公式は非類似度としてKLダイバージェンスを採用して式(10)を解くことに相当する。KLダイバージェンスは確率分布間の非類似度を測る尺度であるが、ここでは平均は共通で共分散行列として{\bf B}{\bf B}_{k}を用いた多次元正規分布同士のKLダイバージェンスD({\bf B}, {\bf B}_{k})とする。多次元正規分布のKLダイバージェンスについては以下が参考になる。

qiita.com

これを踏まえてD({\bf B}, {\bf B}_{k})を式で表すと、

\displaystyle 
\begin{aligned}
D({\bf B}, {\bf B}_{k}) &= D_{KL} \left[{\mathcal N}({\bf \mu}, {\bf B}) \parallel {\mathcal N}({\bf \mu}, {\bf B}_{k}) \right]  \\
&= {\rm tr}({\bf B} {\bf B}_k^{-1}) - \log \det ({\bf B} {\bf B}_k^{-1}) - n 
\end{aligned}
\tag{11}

となり、この式(11)を式(10)に当てはめて最適化問題を解くと式(8)が得られる*5

ちなみに、KLダイバージェンスは対称性がないため、式(11)の{\bf B}{\bf B}_{k}を入れ替えたKLダイバージェンスD({\bf B}_{k}, {\bf B}) を使うと別の更新式が得られ、こちらはDFP公式と呼ばれている。実はKLダイバージェンスに基づく非類似度にはD({\bf B}_{k}, {\bf B}) = D({\bf B}^{-1}, {\bf B}_{k}^{-1}) という性質があり、式(10)の{\bf B}^{-1}に関する最適化問題を解くことと等価になるため、その最適解である式(9)の{\bf B}_{k}^{-1}{\bf B}_{k}に置き換えたものがDFP公式の更新式となる(従って、その逆行列の更新式は式(8)の{\bf B}_{k}{\bf B}_{k}^{-1}に置き換えたものとなる)。

アルゴリズムの流れ

BFGS公式に基づく準ニュートン法アルゴリズムは以下の通り。

  1. 近似行列の初期値({\bf B}_0)^{-1} *6と初期点{\bf x}^{(0)}を定める。
  2.  \parallel \nabla f(x^{(k)}) \parallelが十分に小さければ終了。
  3. 式(5)で探索方向を求める。
  4. 直線探索で適切なステップ幅\alpha_kを求める。
  5. 式(2)でx^{(k+1)}を求める。
  6. 式(9)で({\bf B}_{k + 1})^{-1}を求める。
  7. 2.に戻る。

実装

以下が実装した準ニュートン法pythonコード。

import numpy as np
import copy

# アルミホ条件を満たすαを線形探索
def backtracking_line_search(objective, x, d, tau1, beta):
    alpha = 1
    g = lambda alpha : objective.func(x + alpha * d)
    g_dash = lambda alpha : objective.grad(x + alpha * d) @ d    
    while(1):
        if ((g(alpha) <= g(0) + tau1 * g_dash(0) * alpha)):
                break
        alpha *= beta
    return alpha

# 目的関数のクラスobjectiveは以下のメンバを持っている必要あり
# objective.func(x) : return f(x)
# objective.grad(x) : return ∇f(x)
# objective.dim : 変数の次元
def quasi_newton(objective, eps = 1e-9, tau1 = 0.5, beta = 0.9):
    H = np.eye(objective.dim) # = inv(B)
    x = np.random.randn(objective.dim)
    
    while(1):
        grad = objective.grad(x)
        if(np.linalg.norm(grad) <= eps):
            break
        d = - H @ grad
        alpha = backtracking_line_search(objective, x, d, tau1, beta)
        x_new = x + alpha * d
        s = (x_new - x)[:, None]
        y = (objective.grad(x_new) - grad)[:, None]
        I = np.eye(H.shape[0])
        sy = (s.T @ y)
        H = (I - (s @ y.T) / sy)  @ H @ (I - (y @ s.T) / sy) + (s @ s.T / sy)
        x = x_new
    return x, objective.func(x)

引数objectiveに目的関数を表すクラスのインスタンスを入力するが、これはf({\bf x})の値を返すメソッドfunc\nabla f({\bf x})を返すメソッドgradを実装してる必要がある。準ニュートン法の処理自体は素直に書いており特別なことはしていない。なお、脚注*3にあるように({\bf B}_k)^{-1}の正定値性を保証するならウルフ条件を満たすように直線探索するべきだが、alphaの初期値次第ではウルフ条件を満たすalphaが見つからず無限ループに陥るなど扱いにくかったのでアルミホ条件のバックトラッキングで妥協した*7

前回同様、上記のソースコードと検証用コードを以下のgithubリポジトリに上げている。

github.com

最後に

今回は制約なしの非線形計画問題を解く準ニュートン法を実装したが、次は制約ありの非線形計画問題を解く内点法を実装する予定。実は内点法の実装自体はもう出来てるので後は記事を書くだけ。。。

参考文献

*1:実装ではアルミホ条件を満たすように、バックトラック法で直線探索している。

*2:逆に狭義凸関数なら正定値行列、凸関数なら半正定値行列になる。

*3: 参考文献である『しっかり学ぶ数理最適化』ではg(\alpha) = f({\bf x}^{(k)} + \alpha {\bf d}({\bf x}^{(k)}) ) g'(\alpha) = 0を達成する\alphaを選択したときの正定値性が示されているが、別の参考文献『機械学習のための連続最適化』ではウルフ条件を満たした\alphaでも正定値になることが示されている。

*4:\displaystyle({\bf A} + {\bf u}{\bf v}^T)^{-1} = {\bf A}^{-1} - \frac{{\bf A}^{-1} {\bf u}{\bf v}^T {\bf A}^{-1}}{1 + {\bf v}^T{\bf A}^{-1} {\bf u}}

*5:この最適化問題は解析的に解ける。

*6:正定値対称行列でなければならないので{\bf I}などを初期値とする。

*7:パウエルの修正BFGS公式を使って更新すればアルミホ条件でも問題ないかも。