甲斐性のない男が機械学習とか最適化とかを綴るブログ

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

ベイズ推論による教師あり分類 ~ガウス混合分布編~

はじめに

以前、ベイズ推論による機械学習入門という書籍を読み、以下のような記事を公開しました。

yamagensakam.hatenablog.com

この書籍には、多次元ガウス混合モデルによるクラスタリング手法が詳しく述べられていますが、多次元ガウス混合モデルによる教師あり分類の手法は少し触れる程度にとどまっています(クラスタリングで述べられている導出過程を流用すれば、比較的簡単に導出できるためだと思います)。
ということで、今回は多次元ガウス混合モデルによる教師あり分類の導出、実装をやってみました。

ガウス混合モデル

モデル自体は教師ありだろうが、教師なしだろうが同じモデルを考えることができます。まず、全て独立でサンプリングされるD次元のデータが{\bf X}  = \left\{ {\bf x}_1 , \cdots, {\bf x}_N \right\}N個あるとし、そのそれぞれに対応するラベルデータを{\bf S} = \left\{ {\bf s}_1 , \cdots, {\bf s}_N \right\} とします。ここで、{\bf s}1 of K表現と呼ばれ、クラス数がKの場合、所属するクラスに対応する要素のみ1で他は全て0K次元ベクトルです。この{\bf s}をサンプルするための分布として、パラメータ{\bf \pi}を持つカテゴリ分布を考えます。
また、K個のガウス分布が混合したモデルを考えるので、それぞれの分布の平均を{\bf \mu}_1, \cdots, {\bf \mu}_K、精度行列(共分散行列の逆行列)を{\bf \Lambda}_1,\cdots, {\bf \Lambda}_Kとします。なお、各パラメータを個別に扱わない場面では表記を簡単にするため、これらのガウス混合分布のパラメータを全てまとめて{\bf \Theta}と表記します。
更に、カテゴリ分布のパラメータ{\bf \pi}や各ガウス分布{\bf \mu}_k{\bf \Lambda}_kといったパラメータもベイズ推論の枠組みでは、確率分布に従うと考えます。ここでは、それぞれ共役事前分布に従うことにします。具体的には、

\displaystyle p({\bf \pi}) = {\rm dir}({\bf \pi} | {\bf \alpha}) \tag{1}
\displaystyle p({\bf \mu}_k, {\bf \Lambda}_k)  = {\rm NW} ({\bf \mu}_k, {\bf \Lambda}_k | {\bf m}, \beta, \nu, {\bf W} ) \tag{2}

という分布に従うとします。ここで、{\rm dir}(\cdot)はディリクレ分布、{\rm NW}(\cdot)ガウス・ウィシャート分布です(具体的な分布の式は書籍を参照ください)。ここで事前分布として共役事前分布を採用しておくと事後分布や予測分布が解析的に求めることができます。{\bf \alpha}{\bf m}, \beta, \nu, {\bf W} は超パラメータで今回は事前に与えるものとします。また、{\bf m}, \beta, \nu, {\bf W} については、それぞれのガウス分布で異なる値を与えることもできますが、簡単のため今回は全て共通の値とします。
さて、ここまで出てきた確率変数{\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}の同時分布を考えると、以下のようになります。
{
\displaystyle
\begin{eqnarray}
p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}) &=& p({\bf X} | {\bf S}, {\bf \Theta}) p ({\bf S}|{\bf \pi})p({\bf \Theta})p({\bf \pi} ) \\
&=& \left\{\prod_{i = 1}^{N} \prod_{k = 1}^{K} {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1})^{{\bf s}_{i, k}} {\rm Cat}({\bf s}_k | {\bf \pi})\right\} \left\{ \prod_{k = 1}^{K} p({\bf \mu}_k, {\bf \Lambda}_k) \right\} p({\bf \pi}) 
\end{eqnarray}
} \tag{3}

ここで、{\rm N}(\cdot)ガウス分布{\rm Cat}(\cdot)はカテゴリ分布を表し、p({\bf \pi})p({\bf \mu}_k, {\bf \Lambda}_k) は式(1)、式(2)で表される分布です。
このモデルのグラフィカルモデルを下図に示します。この図で{\bf X}{\bf S}の四角で囲まれているのが学習データとして与えられるデータ、{\bf x}_*{\bf s}_*がクラスを予測する対象のデータ(テストデータ)です。


f:id:YamagenSakam:20190210190923p:plain

このように図示しておくと有効分離の考え方を用いて、変数間の条件付き独立性がわかりやすくなり、後々式展開する際に便利です。有効分離については、以下の記事がわかりやすいと思います。

条件付き独立と有向分離を用いた統計モデルの妥当性チェック - StatModeling Memorandum

推論

上述式(3)のモデルから、パラメータの事後分布とそれを元に、新たなデータ{\bf x}_*が与えられた時の{\bf s}_*の確率分布である予測分布を導出します。このモデルの場合、上述のように共役事前分布を事前分布として採用しているため、事後分布および予測分布は解析的に求めることができます。ここでは、逐次的にデータが与えられるとして、オンラインで予測分布を更新していけるような式を導出します。

事後分布の導出

今回扱う問題では{\bf X}{\bf S}がデータとして与えられるので、これらが与えられている下でのパラメータの事後分布を求めます。つまり、

 \displaystyle p({\bf \Theta}, {\bf \pi}|{\bf X}, {\bf S}) = \frac{p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta})}{p({\bf X}, {\bf S})}  \varpropto p({\bf X}, {\bf S}, {\bf \pi}, {\bf \Theta}) \tag{4}

を求めることになります。なお、クラスタリングの場合、{\bf S}がデータとして与えられているわけではないので、求めるべき事後分布としてはp({\bf \Theta}, {\bf \pi}, {\bf S}|{\bf X}) となり、こちらは解析的に求めることが難しくなります(なので、ギブスサンプリングや変分推論などの近似手法を使う)。
ここで式(3)を見ると、\Theta\piは完全に分離していることがわかります。つまり、{\bf X}{\bf S}が与えられた下では、\Theta\piは独立であり、

 \displaystyle  p({\bf \Theta} |{\bf X}, {\bf S}) \varpropto  p({\bf X} | {\bf S}, {\bf \Theta})p({\bf \Theta}) \tag{5}
 \displaystyle  p({\bf \pi} |{\bf X}, {\bf S}) =p({\bf \pi} |{\bf S})  \varpropto   p ({\bf S}|{\bf \pi})p({\bf \pi} ) \tag{6}

を求めればよいことになります。
それでは、式(5)を展開してみます。式(5)の右辺について、対数をとってみると、

{
\displaystyle 
\begin{eqnarray}
\ln p({\bf X} | {\bf S}, {\bf \Theta}) p({\bf \Theta}) &=& \ln \prod_{i = 1}^{N} \prod_{k = 1}^{K} {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1})^{{\bf s}_{i, k}} p({\bf \mu}_k, {\bf \Lambda}_k)\\
 &=& \sum_{k = 1}^{K} \left\{ \sum_{i = 1}^{N} s_{i,k} \ln {\rm N}({\bf x}_i| {\bf \mu}_k, {\bf \Lambda}_k^{-1}) + {\rm NW} ({\bf \mu}_k, {\bf \Lambda}_k | {\bf m}, \beta, \nu, {\bf W} ) \right\}
\end{eqnarray}
} \tag{7}

となります。これを見ると各kに対して分離しているのがわかるので、各kに対し{\bf \mu}_k, {\bf \Lambda}_kそれぞれ独立に事後分布を求めればよいことになります。また、式(7)の中括弧の中はガウス・ウィシャート分布を事前分布として用いたガウス分布パラメータの事後分布推定の計算となります。この計算方法は件の書籍に詳しく載っているので、ここでは割愛しますが、最終的な形は以下のようになります。
\displaystyle p({\bf \mu}_k, {\bf \Lambda}_k | {\bf X}, {\bf S}) =  {\rm WA}({\bf \mu}_k, {\bf \Lambda}_k |  \hat{{\bf m}}_k, {\hat \beta}_k, {\hat \nu}_k, {\hat {\bf W}}_k) \tag{8}
ただし、

{
\displaystyle
\begin{eqnarray}
{\hat \beta}_{k} &=& \sum_{i = 1}^{N} s_{i, k} + \beta  \\
{\hat \nu}_{k} &=& \sum_{i = 1}^{N} s_{i, k} + \nu \\ 
{\hat {\bf m}}_{k} &=& \frac{\sum_{i = 1}^{N} s_{i, k} {\bf x}_{i} + \beta {\bf m}  } {{\hat \beta}_{k}} \\
{\hat {\bf W}}_{k} ^{-1} &=&   \sum_{i = 1}^{N} s_{i, k} {\bf x}_{i} {\bf x}_{i}^T + \beta {\bf m}{\bf m}^T - {\hat \beta}_{k}{\hat {\bf m}}_{k}  {\hat {\bf m}}_{k}^{T} + {\bf W}^{-1}\\
\end{eqnarray}
} \tag{9}

です。同様に式(6)はディリクレ分布を事前分布としてカテゴリ分布の事後分布そのものであるため、

\displaystyle p({\bf \pi} |{\bf S}) = {\rm dir}({\bf \pi}| {\hat {\bf \alpha}}) \tag{10}

というディリクレ分布になります。ただし、

\displaystyle  \alpha_{k} = \sum_{i = 1}^{N} s_{i, k} + \alpha \tag{11}

です。

ここで、式(9)および式(11)の各超パラメータをデータが与えられるたびに逐次更新できる形に変更すると、

{
\displaystyle
\begin{eqnarray}
{\hat \beta}_{n, k} &=&  {\hat \beta}_{n-1, k} + s_{n, k} \\
{\hat \nu}_{n, k} &=&  {\hat \nu}_{n-1, k}  + s_{i, k} \\ 
{\hat {\bf m}}_{n, k} &=&  \frac{s_{n, k} {\bf x}_{n} + {\hat \beta}_{n-1, k} {\hat {\bf m}}_{n-1, k}  } {{\hat \beta}_{n, k}}\\
{\hat {\bf W}}_{n, k}^{-1} &=& { \hat {\bf W}}_{n-1, k}^{-1} + s_{i, k}{\bf x}_{i} {\bf x}_{i}^T +{\hat \beta}_{n-1, k}{\hat {\bf m}}_{n-1, k}  {\hat {\bf m}}_{n-1, k}^{T}  - {\hat \beta}_{n, k}{\hat {\bf m}}_{n, k}  {\hat {\bf m}}_{n, k}^{T} \\
{\hat \alpha}_{n, k} &=&  {\hat \alpha}_{n-1, k}  + s_{i, k} 
\end{eqnarray}
} \tag{12}

となり、これによりn-1の時点で、計算した超パラメータを用いて、n番目のデータが来たときの超パラメータに更新することが可能です(初期値は事前分布の超パラメータ{\bf m}, \beta, \nu, {\bf W}, \alpha )。

予測分布の導出

最終的に欲しいのは、新たなデータ{\bf x}_*が与えられた時のクラス分類結果{\bf s}_*です。そのため、{\bf x}_*が与えられた下での{\bf s}_*の予測分布を求めます。予測分布は、

\displaystyle p({\bf s}_* | {\bf x}_*, {\bf X}, {\bf S}  ) = \frac{p({\bf s}_* , {\bf x}_*, {\bf X}, {\bf S}  )}{p( {\bf x}_*, {\bf X}, {\bf S} )} = \frac{p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  )p( {\bf X}|{\bf s}_*, {\bf S})p({\bf s}_*|{\bf S})p({\bf S})}{{p( {\bf x}_*, {\bf X}, {\bf S} )}}  \tag{13}

となります。ここで、{\bf s}_*に関係のない項は無視していいので、p({\bf S})p( {\bf x}_*, {\bf X}, {\bf S} )は考慮から外します。更に、p( {\bf X}|{\bf s}_*, {\bf S})について、グラフィカルモデルから有効分離性を読み解くと{\bf S}が与えられた下では、{\bf X}{\bf s}_*は独立であり({\bf X}に繋がる{\bf S}が観測され、かつ、head-to-head型を形成している{\bf s}_*, {\bf x}_*, {\bf \Theta} において{\bf x}_*が観測されていないため)、p( {\bf X}|{\bf s}_*, {\bf S})= p( {\bf X}| {\bf S}) となるため、こちらも考慮しなくてよくなります。まとめると、式(13)は

\displaystyle p({\bf s}_* | {\bf x}_*, {\bf X}, {\bf S}  ) \varpropto p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  ) p({\bf s}_*|{\bf S}) \tag{14}

と、2つの項のみ考えればよくなります。1つ目の項は、
 \displaystyle p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  ) = \int  p({\bf x}_* | {\bf s}_* ,{\bf X},{\bf S}, {\bf  \Theta} )  p({\bf  \Theta}| {\bf s}_*,{\bf X},{\bf S}) d {\bf \Theta} \tag{15}

であり、ここでもグラフィカルモデルから有効分離性を考えるとp({\bf x}_* | {\bf s}_* ,{\bf X},{\bf S}, {\bf  \Theta} )  =p({\bf x}_* | {\bf s}_* , {\bf  \Theta} ) p({\bf  \Theta}| {\bf s}_*,{\bf X},{\bf S})  = p({\bf  \Theta}| {\bf X},{\bf S})とすることができます。従って、

\displaystyle p({\bf x}_*| {\bf s}_* , {\bf X}, {\bf S}  ) =  \int  p({\bf x}_* | {\bf s}_* , {\bf  \Theta} )  p({\bf  \Theta}| {\bf X},{\bf S}) d {\bf \Theta} \tag{16}

となります。
 p({\bf  \Theta}| {\bf X},{\bf S}) は、式(5)、式(8)、式(9)に記している事後分布そのものです。今、各kに対して{s}_{*,k} = 1であるそれぞれの予測分布を求めたいので、式(16)より

\displaystyle p({\bf x}_*| {\bf s}_{*,k} = 1 , {\bf X}, {\bf S}  ) =  \int  \int p({\bf x}_* | {\bf  \mu}_k, {\bf \Lambda}_k )  p({\bf  \mu}_k, {\bf \Lambda}_k | {\bf X},{\bf S}) d {\bf  \mu}_k d {\bf  \Lambda}_k \tag{17}

を求めればよいことになります。これは、ガウス・ウィッシャート分布を事前分布として用いた場合のガウス分布の予測分布計算になります。こちらも予測分布導出は書籍に載っているので、導出過程は省略し最終形を示します。

\displaystyle p({\bf x}_*| {\bf s}_{*,k} = 1 , {\bf X}, {\bf S}  ) = {\rm St}({\bf x}_* |{\hat {\bf m}}_{n, k}, \frac{(1-D+ {\hat \nu}_{n, k}){\hat \beta}_{n, k} }{1 + {\hat \beta}_{n, k} } {\hat {\bf W}}_{n,k}, 1-D+ {\hat \nu}_{n, k}) \tag{18}

ただし、 {\rm St}(\cdot)は多次元のスチューデントのt分布で、D{\bf x}の次元です。
次に、式(14)のp({\bf s}_*|{\bf S})ですが、こちらも同様に、

\displaystyle p({\bf s}_*|{\bf S}) = \int p({\bf s}_*|{\bf S},{\bf \pi})p({\bf \pi}| {\bf S}) d {\bf \pi} =  \int p({\bf s}_*|{\bf \pi})p({\bf \pi}| {\bf S}) d {\bf \pi} \tag{19}

となります。ただし、第2式から第3式への展開はグラフィカルモデルより {\bf \pi}を観測した下での{\bf S} {\bf s}_*の条件付き独立性を利用しました。そしてこの式(19)は、ディリクレ分布を事前分布として使ったカテゴリ分布の予測分布です。ということで、この予測分布についても導出過程は割愛し最終形を示します。

\displaystyle p({\bf s}_*|{\bf S})  = {\rm Cat}({\bf s}_*|{\hat {\bf \eta}}_{n}) \tag{20}

ただし、{\hat {\bf \eta}_{n}}の各要素{\hat \eta}_{n,k}
\displaystyle {\hat \eta}_{n,k}=  \frac{{\hat \alpha}_{n, k}}{\sum_{i=1}^K  {\hat \alpha}_{n, i}}  \tag{21}
です。
ここまでの話をまとめると、n番目のデータが与えら手た時の予測分布の更新式は、式(14)、式(17)、式(20)、式(21)より

 \displaystyle p(s_{*, k} = 1 | {\bf x}_*, {\bf X}, {\bf S}  )  \varpropto  {\hat \eta}_{n,k} {\rm St}({\bf x}_* |{\hat {\bf m}}_{n, k}, \frac{(1-D+ {\hat \nu}_{n, k}){\hat \beta}_{n, k} }{1 + {\hat \beta}_{n, k} } {\hat {\bf W}}_{n,k}, 1-D+ {\hat \nu}_{n, k}) \tag{22}

となります(各超パラメータは式(12)、式(21))。そして、各クラスの所属確率の平均を得たいときは、k = 1 \cdots Kに対し、式(22)の右辺に{\bf x}_*を代入し、得られる値を合計が1になるように正規化すればよいです。

実装

ということで実装です。

import numpy as np
import scipy as sci
import math
import random
import matplotlib.pyplot as plt
import matplotlib.animation as anm

class multivariate_student_t():
    def __init__(self, mu, lam, nu):
        self.D   = mu.shape[0]
        self.mu  = mu
        self.lam = lam
        self.nu  = nu

    def pdf(self, x):
        temp1 = np.exp( math.lgamma((self.nu+self.D)/2) - math.lgamma(self.nu/2) )
        temp2 = np.sqrt(np.linalg.det(self.lam)) / (np.pi*self.nu)**(self.D/2)
        temp3 = 1 + np.dot(np.dot((x-self.mu), self.lam),  (x-self.mu))/self.nu
        temp4 = -(self.nu+self.D)/2
        return temp1*temp2*((temp3)**temp4)

class mixture_gaussian_classifier:
    def __init__(self, alpha0, m0, invW0, beta0, nu0):
        self.K = alpha0.shape[0]

        #予測分布(多次元t分布)の生成
        self.pred_distr = (lambda m, W, beta, nu:
                                      multivariate_student_t
                                      (m,
                                      (1 - m.shape[0] + nu) * beta * W/(1 + beta),
                                      (1 - m.shape[0] + nu)))

        #超パラメータの初期化
        self.alpha = alpha0
        self.eta = alpha0/np.sum(alpha0)
        self.m = [None] * self.K
        self.invW = [None] * self.K
        self.beta = np.zeros(self.K)
        self.nu = np.zeros(self.K)
        for k in range(self.K):
            self.m[k] = m0
            self.invW[k] = invW0
            self.beta[k] = beta0
            self.nu[k] = nu0



    def train(self, x, s):
        #各超パラメータの更新
        k = np.where(s == 1)[0][0]
        self.alpha[k] += 1.0
        self.eta = self.alpha/np.sum(self.alpha)
        pre_beta = self.beta[k]
        self.beta[k] += 1.0
        self.nu[k] += 1.0
        pre_m = self.m[k]
        self.m[k] = (x + pre_m * pre_beta)/self.beta[k]
        self.invW[k] = (self.invW[k] + np.dot(x.reshape(-1,1), x.reshape(-1,1).T)
                        - self.beta[k] * np.dot(self.m[k].reshape(-1,1), self.m[k].reshape(-1,1).T)
                        + pre_beta * np.dot(pre_m.reshape(-1,1), pre_m.reshape(-1,1).T))

    def test(self, x):
        s = np.zeros(self.K)
        #各クラスの予測平均を算出
        for k in range(self.K):
            s[k] = (self.eta[k] *
                    self.pred_distr(self.m[k], np.linalg.inv(self.invW[k]),
                                    self.beta[k], self.nu[k]).pdf(x))
        s /= np.sum(s)
        return s

クラスmixture_gaussian_classifierは、ここまで述べてきた学習(train)と予測平均の算出(test)を行うクラスです。また、multivariate_student_t多次元のスチューデントのt分布のクラスで、この実装については、以下の記事の内容を参考にさせてもらいました。

<ベイズ推論> 多次元ガウス分布の学習 - Helve’s Python memo

実験・結果

ということで、以下のようなコードを組んで実験してみました。クラス数はK=3で、学習データ数はそれぞれN_1=50, N_2 = 60, N_3 = 70と少し偏りを持たせています。
また、せっかく事後分布を逐次学習できるようにしたので、学習データが与えられるたびに予測平均がどのように変化するかをアニメーションにしてみました。

global used_idx

def update(i, idx, mgc, X, S):
    plt.cla
    test_range = range(0,60)
    X_new = np.zeros((2, len(test_range)**2))
    S_new = np.zeros((3, len(test_range)**2))
    i = idx[i]

    #データの学習
    mgc.train(X[:, i],S[:, i])

    used_idx.append(i)
    c = S[:, used_idx]

    plt.scatter(X[0, used_idx], X[1, used_idx],
                edgecolors="black", marker="o" , color=S[:, used_idx].T)

    count = 0
    S_res = np.zeros((len(test_range), len(test_range), 3))
    for x1 in test_range:
        for x2 in test_range:
            #各座標位置に対する予測平均を算出
            X_new[:, count] = np.array([float(x1), float(x2)])
            S_new[:, count] = mgc.test(X_new[:, count])
            S_res[int(x2 - test_range[0]), int(x1 - test_range[0]), :] = S_new[:, count]
            count += 1

    plt.imshow(S_res)


def main():
    N1 = 50
    N2 = 60
    N3 = 70
    X = np.zeros((2, N1+N2+N3))
    S = np.zeros((3, N1+N2+N3), 'int8')

    #クラス1の学習データ
    mu1 = np.array([45, 25])
    A = np.array([[3,0],[0,8]])
    sigma1 = np.dot(A.T, A)
    X[:, :N1] = np.random.multivariate_normal(mu1, sigma1, N1).T
    S[0,:N1] =1

    #クラス2の学習データ
    mu2 = np.array([35, 40])
    A = np.array([[5, 3],[3,7]])
    sigma2 = np.dot(A.T, A)
    X[:, N1:N1 + N2] = np.random.multivariate_normal(mu2, sigma2, N2).T
    S[1,N1:N1 + N2] = 1

    #クラス3の学習データ
    mu3 = np.array([30, 25])
    A = np.array([[7,-5],[-5,7]])
    sigma3 = np.dot(A.T, A)
    X[:, N1 + N2:]= np.random.multivariate_normal(mu3, sigma3, N3).T
    S[2,N1 + N2:] = 1

    #事前分布の超パラメータ
    alpha0 = np.random.randn(3)
    m0 = (mu1 + mu2 + mu3)/2
    invW0 = np.linalg.inv((sigma1 + sigma2+ sigma3)/3)
    beta0 = 1.0
    nu0 = m0.shape[0]

    #推論オブジェクトの生成
    mgc = mixture_gaussian_classifier(alpha0, m0,invW0, beta0, nu0)

    fig = plt.figure()
    idx = [i for i in range(X.shape[1])]
    random.shuffle(idx)

    global used_idx
    used_idx = []
    #結果のアニメーション表示(300msごとに関数updateを呼び出す)
    ani = anm.FuncAnimation(fig, update, fargs = (idx, mgc, X, S),
                            interval = 300, frames = N1 + N2+ N3)
    ani.save("result.gif", writer="imagemagick")

if __name__ == '__main__':
    main()

そして、こちらが結果です。

f:id:YamagenSakam:20190210212509g:plain

この図で、各色がそれぞれのクラスを表し、追加されていく各点が学習データを表しています。また、背景色が予測分布が算出した各座標点の各クラスへの所属確率を表しています。
特に左上に辺りの領域に着目すると、最初の方は青色のクラスと予測されていましたが、青と緑のクラスのデータを観測していくにつれ、緑色のクラスだと予測されるようになっていくのがわかります。これは青色のクラスが左下から右上にかけての共分散が大きいのに対し、緑色のクラスの共分散が左上から右下への共分散が大きく、その特徴をサンプルが増えるにつれ、正しく推論できるようになっていくためだと考えられます。

まとめ

今回、ガウス混合分布を用いた教師あり分類をベイズ推論の枠組みを適用して実装してみました。教師ありの分類の場合は、ガウス混合分布モデルによる予測は解析的に求めることができます。
ベイズ推論の枠組みで分類を行う場合は他にもロジスティック回帰をベイズの枠組みで定式化する方法などがあります。今後は、このベイズロジスティック回帰を含めもう少し複雑なモデルを学び実装していきたいと思います。

ベイズ推論による機械学習入門を読んだので軽くメモ

はじめに

タイトルの通りです。以下の本を一通り読んだので、再度読む際の手助けになるようメモを残します。

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

まだ理解が怪しい部分も多々ありますが、取りあえず忘れる前に。。。

全体を通してのポイント

ベイズ推論の枠組み

ベイズ推論においては、観測データも分布のパラメータ、回帰等の重み係数も潜在変数も確率変数として扱う。観測されたデータをDとし、未知の確率変数(分布のパラメータ、重み係数、潜在変数等々)を{\bf \theta}とした時、以下の2ステップを行うのがベイズ推論の枠組み。

  1. モデル\displaystyle p(D, {\bf \theta})の構築(モデリング
  2. 事後分布\displaystyle p( {\bf \theta} | D) = \frac{p(D, {\bf \theta})}{p({D})}の計算

この事後分布を求めることがベイズ推論における「学習」。ベイズ推論は全てこの枠組みの中で扱われるので、後は如何にモデリングするかと、如何に求めたい確率変数の事後分布を計算するか、という話になる。モデルが複雑な場合、事後分布は解析的に求められない場合もあるので、その時は変分推論やギブスサンプリングなどを使って近似的に事後分布を求める。

予測分布

更に、求めた事後分布を使って、未観測のデータ{\bf x}_*に対し、何らかの知見を得たい場合もある(知見の内容はタスクによって異なる。例えば、回帰の場合は、入力{\bf x}_*に対応する出力y_*を求める等)。その時は以下の予測分布を計算する。


\displaystyle p({\bf x}_* | D) = \int p( {\bf x}_* | {\bf \theta}) p( {\bf \theta} | D) d  {\bf \theta}

1章

基本的にこの章は後の章で語られる内容のアウトラインが載っている。この章で語られる内容で特に重要なのは、上記のベイズ推論のポイントと以下の関係式。

  • 周辺化 \displaystyle p(y) = \int p(x, y) dy
  • 条件付き分布 \displaystyle p(x | y) = \frac{p(x, y)}{p(y)}

また、独立の諸概念も重要。

  • 変数間の独立性p(x, y) = p(x)p(y)(これが成り立つ時、xyは独立という)
  • 条件付き独立p(x,z|y) = p(x|y)p(z|y)yを観測した下でxzは独立だが、yを観測しなければxzは独立でないことに注意)

2章

確率密度関数については書籍を参照するとして、取りあえず期待値が大事。エントロピーとかKLダイバージェンスとかも大事だけど、期待値がわからないと後々お話にならなくなる。特に、


\displaystyle \langle  f({\bf x}) \rangle_{p({\bf y})} = \int f({\bf x}) p({\bf y}) d{\bf y} =  f({\bf x})  \int p({\bf y}) d{\bf y}  =  f({\bf x })
\displaystyle \langle  f({\bf x}) \rangle_{p({\bf x}, {\bf y})} = \int \int f({\bf x}) p({\bf x}, {\bf y}) d{\bf x}d{\bf y} =  \int f({\bf x}) \int  p({\bf x}, {\bf y}) d{\bf y} d{\bf x}  = \langle f({\bf x }) \rangle_{p({\bf x})}

と、期待値を取る関数に出てこない確率変数は、期待値計算から外せることは重要。
また、同時分布の期待値とか、条件付き期待値とかは込み入ってくると混乱するので注意。迷った時は期待値の定義、


\displaystyle \langle   f({\bf x}) \rangle_{p({\bf x})}  = \int f({\bf x}) p({\bf x}) d{\bf x}

に立ち返るのが吉。

確率密度関数の対数を取った関数\ln p({\bf x})、基本的な期待値( \langle  {\bf x} \rangleとか \langle \ln {\bf x} \rangleとか)も、後々の章(特に4章、5章の変分推論辺り)で何度も使う。本当に何度も使うので、何度も見返す。

3章

2章で挙がっている確率分布の各パラメータおよび線形回帰の重み係数をベイズ推論の枠組みで推定する章。事前分布として共役事前分布を使うので、(この章で示されている問題については)事後分布が解析的に計算でき、しかも事前分布として同じ形の分布で求まるので、データを得るたびに逐次的に更新することも可能。考え方としては、

  1. 学習したいパラメータ\thetaの事前分布p(\theta)として共役事前分布を設定
  2. こうすると事後分布は事前分布と同じ形になり、p( {\bf \theta} | D) = \frac{p(D, {\bf \theta})}{p({D})}= \frac{p(D|{\theta})p(\theta)}{p({D})} \propto p(D|{\theta})p(\theta) としてp(D)は定数として扱うことが可能(後から帳尻を合わせる)。
  3. 尤度p(D | \theta)p(D | \theta) = \prod_{i=1}^N p(x_i|\theta)と観測データが独立だと仮定(後の章で観測データ間にも独立ではないモデルもあり)

また、計算のテクニックとして、以下のようなものが使われている。

  • 事後分布が事前分布と同じ形になることから、\ln p(D|{\theta})p(\theta)をうまく展開し、2章で計算済みの事後分布の対数をとった式のパラメータと対応を取ることで、事後分布のパラメータを計算
  • 予測分布を求める際、積分の計算が難しい時は、ベイズの定理に立ち返りp({\bf x}_*)=\frac{p({\bf x}_* | \theta) p(\theta)}{p(\theta | {\bf x}_* )}となることから、上記の事後分布計算と同じく対数をとって分布の式とパラメータの対応を取ればよい
  • 予測分布計算時、いきなり事後分布のパラメータで予測分布求めようとすると、ややこしくなるケースもあり、その場合はまず事前分布を用いて予測分布を計算し、後から対応するパラメータを事後分布のものに置き換えればよい

4章

生成モデル

観測データの生成過程を仮定して、確率分布p(D, {\bf  \theta})モデリングしたものが生成モデル。例えば、観測されたデータに多峰性が見れられる場合、これを単一の分布のみで表現しようとしても、うまく表現しきれない。このような多峰性のあるデータは、

  1. K個のクラスタがあり、観測データ{\bf x}はそのいずれかのクラスタ(例えば、\theta_kを持つk番目のクラスタ)から生成されたと仮定p(\bf{x} | \theta_k)
  2. {\bf x}が割り当たっているクラスタ{\bf s}  \in \left\{0, 1 \right\}^K \ \ {\rm s.t.} \ \ \sum_{i=1}^K s_i = 1(割り当たっているクラスタの要素のみ1、それ以外0)で表現すると、{\bf x}の確率分布はp({\bf x}| {\bf s} ,\Theta) = \prod_{k=1}^K p(\bf{x} | \theta_k)^{s_k}
  3. このクラスタの割り当て{\bf s}も割り当て比率パラメータ\piに支配されるカテゴリ分布の確率変数と仮定p({\bf s}| \pi)
  4. クラスタのパラメータの事前分布p(\Theta)を仮定(ここでも共役事前分布使うと幸せになれる)
  5. クラスタ割り当てのカテゴリ分布の事前分布p(\pi)を仮定(こちらもカテゴリ分布の共役事前分布であるディリクレ分布を使うと幸せ)
  6. 以上のことから、データ生成のモデルをp({\bf X}, {\bf S}, \Theta, \pi) = p({\bf X}|{\bf S} , \Theta) p({\bf S}| \pi)  p(\Theta) p(\pi) = \prod_{i = 1}^N p({\bf x}_i| {\bf s}_i ,\Theta) p({\bf s}_i| \pi)  p(\Theta)p(\pi)と仮定

等のようにモデリングする。
このモデルから、観測データに対するクラスタ割り当ての事後分布\displaystyle p({\bf S}|{\bf X}) = \int \int \frac{p({\bf X}, {\bf S}, \Theta, \pi) }{p({\bf X})} d\Theta d\pi を計算したいが、モデルが複雑すぎてp({\bf X})が解析的に求まらない。そのため、近似的に事後分布の推定を行う。

ギブスサンプリング

上のモデルを例にすると、

  1. \Theta, \piが既に観測されているものとして(2でサンプリングした\Theta, \pi を使う)、{\bf S} \sim p({\bf S}| {\bf X}, \Theta, \pi)とサンプリング
  2. {\bf S}が既に観測されているものとして(1でサンプリングした{\bf S}を使う)、\Theta, \pi \sim p(\Theta, \pi | {\bf X}, {\bf S} )とサンプリング
  3. この操作を繰り返す

この繰り返し回数が十分に多い場合、真の分布に近づくことが知られているため、繰り返し続けると所望のデータをサンプリングできるようになる。もちろん、サンプリングする対象となる分布p({\bf S}| {\bf X}, \Theta, \pi)p(\Theta, \pi | {\bf X}, {\bf S} )は十分にサンプリングしやすい分布である必要があるので、式を展開して2章で挙がっているような基本的な分布を導かないといけない(モデルによってはそれができないこともあるので、共役事前分布を使うと幸せになるというのはここに繋がる)。

変分推論


\displaystyle p({\bf S}, \Theta, \pi | {\bf X}) \simeq q({\bf S}) q(\Theta, \pi)

等のように、各確率変数間の独立性を仮定した分布で近似するのが変分推論。真の分布とのKLダイバージェンスが最小となるq({\bf S}) q(\Theta, \pi)を求めればよく、それは


\displaystyle \ln q({\bf S})= \langle \ln p({\bf X}, {\bf S}, \Theta, \pi ) \rangle_{q(\Theta, \pi)}
\displaystyle \ln q(\Theta, \pi) = \langle \ln p({\bf X}, {\bf S}, \Theta, \pi ) \rangle_{q({\bf S})}

として計算できる(あくまで、上記モデルに対する近似式。もう少し一般的な表現は書籍参照)。こちらも一方の結果を用いで、もう一方の計算を行う形になっているため、更新を十分な回数繰り返す必要がある。また、実際には、下記のような式展開して期待値を解析的に求める。

  1. 変分推論を基に、各近似分布q(\cdot)の形を求める。2章に述べられている期待値の性質(上記の無関係な確率変数が期待値計算から外せる性質含む)をうまく使い、q(\cdot)が何分布の関数かがわかるまで式を展開する。
  2. 全ての近似分布の形がわかったら、それぞれの期待値を計算する。これまた2章で計算済みの各分布の基本的な期待値を流用する。

ここでもギブスサンプリングと同様、期待値を解析的に計算するために、2章で挙がっているような基本的な分布を導く必要があるので、やはり共役事前分布を使った方がハッピー。

崩壊型ギブスサンプリング

上記のギブスサンプリングでは、パラメータ\Theta, \pi についてもサンプリングで求めたが、


\displaystyle p({\bf X}, {\bf S}) = \int \int p({\bf X}, {\bf S}, \Theta, \pi) d\Theta d\pi

として周辺化すれば、p({\bf S}| {\bf X})からサンプリングするだけでよくなる。ただし、周辺化した結果が十分にサンプリングしやすい確率分布になることが条件。このモデルの例では、周辺化した結果{\bf s}_1 \cdots, {\bf s}_Nが互いに依存し合うようになり(元々は{\bf X},\Theta, \piを与えられた下での条件付き独立)、このままではN個の潜在変数を同時にサンプリングしなければならず計算量が莫大になる。この問題をp({\bf s}_n | {\bf X}, {\bf S}_{ \backslash n})と、サンプリングしたい{\bf s}_n以外の潜在変数{\bf S}_{ \backslash n}は全て観測済みとして、個別にサンプリングするアイデアで解決。

5章

各種応用に対して、モデルの構築と変分推論(一部ギブスサンプリング)を導出している。この中のいくつかは実際に実装してみようと思うので、ここでは考え方のみメモ。

線形次元削減

観測されたデータ{\bf y}が低次元の潜在変数{\bf x}を元に、 {\bf W}^T {\bf x} + {\bf \mu}という線形変換で生成されると考える。また、実際にはノイズが乗っているはずなので、ガウス分布N({\bf y}| {\bf W}^T {\bf x} + {\bf \mu}, \sigma^2 {\bf I})で生成されると仮定。更に他のパラメータについても事前分布がガウス分布に従うと仮定(\sigma^2は既知のものとして話を進めているが、ガンマ事前分布など使って推定してもよい)。後は変分推論の期待値計算をすれば低次元の潜在変数{\bf x}や変換行列 {\bf W}を推測することができる。

非負値行列因子分解

観測された行列{\bf X}が全ての要素非負であるとし、{\bf X} \simeq {\bf W} {\bf H}という非負値の行列に近似分解することを考える。ここで、S_{dmn} \simeq W_{dm}H_{mn}となる潜在変数{\bf S}を考えることで、効率的な変分推論が導き出せる。また、非負値なので{\bf S}ポアソン分布、{\bf W}, {\bf H}の事前分布はガンマ分布を考える。

隠れマルコフモデル

これまで観測データ間はたがいに独立して得られるものとして扱っていたが、時系列データ等ではその仮定は成り立たない。隠れマルコフモデルでは、

  1. 観測されるデータ{\bf x}_nは現在の状態{\bf s}_nに基づいてデータが観測される(状態{\bf s}_nクラスタ割り当ての潜在変数として持つ混合分布)。
  2. データを観測したら、遷移確率に基づいて次の状態に遷移する({\bf s}_n \rightarrow {\bf s}_{n + 1})。

というモデルを構築する。変分推論は状態の潜在変数{\bf S}の分布を\prod_i^N q({\bf s}_i)と分離して近似する完全分解変分推論とq({\bf S})と分離せずに近似する構造化変分推論がある(その中間であるミニバッチを使った近似もあり)。前者は特別なことをせずに、変分推論の期待値を導き出せるが、後者はうまいこと式変形してフォワードバックワードアルゴリズムと呼ばれる動的計画法が使える形に持っていかないと、計算量が増大する。

トピックモデル

複数の文書を観測データとして、単語の出現頻度を元に文書を解析。各文書中には複数の潜在トピックが存在すると仮定。文書d中のn番目に出てくる単語{\bf w}_{dn}には、どのトピックから発生したかというトピック割り当て潜在変数{\bf z}_{dn}が紐づき、各単語は対応するトピックkのカテゴリ分布に従うと仮定。また、トピック割り当て自体もカテゴリ分布に従うと仮定し、これらのカテゴリ分布のパラメータはディリクレ事前分布に従うモデルを構築。後は頑張って変分推論の期待値計算。

テンソル分解

協調フィルタリングに時間軸を導入し、その時のトレンドを踏まえた推薦を行うようにしたモデル。時間方向の軸には前回の状態を平均とするガウス分布を仮定(線形動的システム)。ここでは、完全分解変分推論を考えているので、推論計算自体はこれまでと同様。

ロジスティック回帰、ニューラルネットワーク

ロジスティック回帰、ニューラルネットワークベイズ推論版。これまでの近似(平均場近似)とは違い、ガウス分布で事後分布を近似。変分パラメータを用いた勾配法で真の分布とのKLダイバージェンス最小化を行う。この勾配法を使うため、再パラメータ化トリックを使うが、これが良く分からない。

おわりに

兎にも角にもわかりやすいです。式展開がわからなくなって何度も行ったり来たりしましたが、PRML10章で卒倒した僕でも割りと短期間で最後まで行けたので、ベイズ推論の入り口として、適切な書籍の1つだと思います。ただ機械学習自体の入門としては少ししんどいので、機械学習全体を俯瞰したのち、その1つの枠組みとしてベイズ推論があると捉えられるようになった上で、読むべきかなと思います。

2019年の目標

はじめに

あけましておめでとうございます。今年は『有言実行』というものに挑戦してみようと思い、まず今年の目標を立てそれを本記事に晒します*1。後半雑になっていきますが、目標とそれを達成するためのアプローチを書いたつもりです。

競プロ

  • atcoder青になる
    • 動的計画法に弱いので典型問題を解く
    • 300点問題を解くのが遅いので、速度重視で300穴埋め
    • 400~600穴埋め
    • 主要言語をpythonからC++に乗り換え
    • セグメント木、ネットワークフロー系を理解

数学・機械学習

その他

  • 部屋をきれいに保つ
    • 週1で掃除機をかける
    • 定期的にゴミ出しをする(特にペットボトルを溜め込まない)
    • ゴミ袋をケチらない
  • 体重5kg減
    • 1日の歩数10000歩以上(現状維持)
    • 週1でジム(ランニング・筋トレ中心から水泳中心に変更?)
    • カロリーを意識(朝:400kcal以内、昼:600kcal以内、夜:800kcal以内)
    • 春からジョギング再開
  • DTM始める
    • とりあえずやってみる
  • ブログの書きかけ記事消化(タイトルだけ書いてるものを合わせると15記事ほど)
    • 気が向いたやつから
  • 行き詰ったら、すぐtwitter、slackを見る癖を直す
    • どうすれば・・・

興味あるけどやれてなくて、たぶん今年もやらないだろうなぁと思うやつ

  • ゲームAIコンテスト
    • たぶんこの中だと1番やる可能性高い
  • ラソン(競プロ)
    • これもやる可能性が高い方
  • ラソン(物理)
    • 数年前からハーフに出たいなーとか思いつつ、年々体力が減少して今や10kmすら完走が危うい・・・
  • 量子情報
    • 社内で勉強会とかがあればやるかも
  • 現代制御
  • 離散最適化、離散凸解析
  • 強化学習
  • スマホアプリ開発
  • 婚活
  • 楽器(ベース辺り)
    • ベースは数年前に衝動買いしたから機材はある
  • 格闘技
    • 運動音痴おじさんでも始められるヤツがあれば・・・
  • 資産運用
    • 証券口座は開設したからなんかやるかー

そういえば・・・

確かこのブログでは言っていなかったと思いますが、僕も昨年8月に電機メーカーを退職し、9月よりIT企業にて速いアルゴリズム考えたり、速いプログラム書いたりするお仕事に従事しています。

*1:仕事の目標は10月に職場で立てたので、プライベートの目標のみ

競プロとかに使うアルゴリズム実装メモ(二分探索、2次元累積和、しゃくとり法)

はじめに

アルゴリズムメモ第3段です。今回は二分探索法、2次元累積和、しゃくとり法と様々な問題に使える汎用的なアルゴリズムを書いていきます。
今回は勉強のため、アルゴリズムの本質的な部分を記述した抽象クラスと実際の問題を解く具象クラス(関数オブジェクト)を分け、テンプレートメソッドパターンもしくはストラテジパターンを用いてコーディングしました。

二分探索

用途

  • 広義単調増加関数または広義単調減少関数f(x)が条件を満たし、かつ、条件成立と不成立の境界となる点xを探索(条件成立と不成立の境界が1つのみ)
  • ソート済み配列から条件が成立する値が存在するかの判定

アルゴリズムのポイント

  1. 変数lを探索範囲の最小値、変数rを最大値で初期化
  2. mr,lの中点(例えば、m = (l + r)/2で求める)として、f(m)が条件を満たすか否かを判定
    • このとき、見つけたい点(条件成立と不成立の境界)がmより小さいか、大きいかはf(x)が広義単調増加(減少)関数なのでわかる
  3. 見つけたい点がm以上ならl = m、そうでなければr = mとする
  4. この2~3の操作をr - l \leq  \epsilon  を満たすまで繰り返す(\epsilonは収束判定基準の定数で、問題に応じて適切な値を定める)
  5. 収束したら、その時点のl,rから最終的な出力値を決める

抽象クラスの実装

上記アルゴリズムを実現する抽象クラスを作成し、見つけたい点が現時点以上かの判定、条件判定、中点計算、戻り値計算を純粋仮想関数とし、上記アルゴリズムをテンプレートメソッドパターンで実装。

#include <vector>

template <typename T>
class abstract_binary_search
{
private:
    virtual bool is_in_right(T x) = 0; //求めているのがxより大きいか
    virtual bool is_converge(T L, T R) = 0; //収束判定
    virtual T ret_val(T L, T R) = 0; //見つからなかった場合に返す値
    virtual T get_middle(T L, T R) = 0; //中央値を返す

protected:
    T search_rec(T L, T R);

public:
    virtual T search() = 0;
    virtual bool is_match(T x) = 0;  //それが求めている値か

    virtual ~abstract_binary_search();
};

template <typename T>
T abstract_binary_search<T>::search_rec(T L, T R) {
    if (is_converge(L, R) == true) {
        return ret_val(L, R);
    }

    T mid = get_middle(L, R);

    if (is_in_right(mid) == true) {
        return search_rec(mid, R);
    }

    return search_rec(L, mid);

}


template <typename T>
abstract_binary_search<T>::~abstract_binary_search()
{
}

具体的な問題を解く実装

Lower bownd

素数nの数列が与えられて、指定された値q以上の値が入っているindexを求める問題。もし数列内にq以上の値がなければnを返す。

#include <iostream>

class lower_bound :
    public abstract_binary_search<int>
{
private:
    std::vector<int> data_list;//検索する数列
    int N; //数列のサイズ
    int search_data;//検索対象データ
    virtual bool is_in_right(int x); //求めているのがxより大きいか
    virtual int ret_val(int L, int R); //見つからなかった場合に返す値
    virtual bool is_converge(int L, int R); //収束判定
    virtual int get_middle(int L, int R);
public:
    lower_bound(const std::vector<int>& data_list);
    void set_search_data(int a_search_data);
    virtual int search();
    virtual bool is_match(int x);  //それが求めている値か
    ~lower_bound();
};

bool lower_bound::is_in_right(int x) {
    return data_list[x] <= search_data;
}


int lower_bound::ret_val(int L, int R) {
    if (search_data <= data_list[L])
        return L;
    else
        return L + 1;
}

int lower_bound::get_middle(int L, int R) {
    return (L + R) / 2;
}

lower_bound::lower_bound(const std::vector<int>& data_list):
    data_list(data_list)
{
    N = data_list.size();
}

int lower_bound::search() {
    return search_rec(0, N);
}

bool lower_bound::is_match(int x) {
    if (x >= N) return false;
    return data_list[x] == search_data;
}

bool lower_bound::is_converge(int L, int R) {
    return R - L <= 1;
}

void lower_bound::set_search_data(int a_search_data) {
    search_data = a_search_data;
}

lower_bound::~lower_bound()
{
}

int main()
{

    int n, q;
    std::cin >> n;
    std::vector<int> S(n);
    
    for (int i = 0; i < n; i++) {
        int  s;
        std::cin >> s;
        S[i] = s;
    }
    std::cin >> q;

    int count = 0;
    
    lower_bound lbs(S);
    for (int i = 0; i < q; i++) {
        int  t;
        std::cin >> t;
        lbs.set_search_data(t);
        
        int res = lbs.search();
        std::cout << res << std::endl;
        if (lbs.is_match(res)) 
            count++;
    }
    std::cout << count << std :: endl;
    return 0;
}
POJ 1064:Cable master

蟻本に二分探索の例として載っている問題。問題はリンクを参照。xが切り取る紐の長さ、f(x)が作れる同じ長さ(長さ=x)の数だとすると、f(x)は広義単調減少関数なので二分探索法が使える。
つまり、f(x)が条件を満たしつつ、xが最大となる点を二分探索する(条件を満たさなければxを小さくし、満たすならまだ長くできると考えxを大きくする)。

#include <vector>
#include <algorithm>
#include <functional>
#include <iostream>
#include <iomanip>
#include <math.h>

class cable_master :
    public abstract_binary_search<double>
{
private:
    double tol;
    std::vector<double> len_list;
    int K;
    int N;
    bool is_condition(double x);
    virtual bool is_in_right(double x); //求めているのがxより大きいか
    //virtual bool is_in_left(double x); //求めているのがxより小さいか
    virtual double ret_val(double L, double R); //見つからなかった場合に返す値
    virtual double get_middle(double L, double R);
    virtual bool is_converge(double L, double R); //収束判定
public:
    cable_master(const std::vector<double>& len_list, int K, double tol);
    virtual double search();
    virtual bool is_match(double x);
    ~cable_master();
};

bool cable_master::is_condition(double x) {
    int count = 0;
    for (int i = 0; i < N; i++) 
        count += int(len_list[i] / x);

    if (count >= K)
        return true;

    return false;
}

bool cable_master::is_in_right(double x) {
    return is_condition(x);
}

double cable_master::ret_val(double L, double R) {
    return R;
}

double cable_master::get_middle(double L, double R) {
    return (L + R) / 2.0; 
}

bool cable_master::is_converge(double L, double R) {
    return R - L < tol;
}

cable_master::cable_master(const std::vector<double>& len_list, int K, double tol) :
    len_list(len_list),
    K(K),
    tol(tol)
{
    N = len_list.size();
}

double cable_master::search(){
    return search_rec(0.0, *std::max_element(len_list.begin(), len_list.end()));

}

bool cable_master::is_match(double x) {
    return false;
}

cable_master::~cable_master()
{
}



int main()
{

    std::vector<double> cable_list;
    int N, K;
    std::cin >> N >> K;

    for (int i = 0; i < N; i++) {
        double len;
        std::cin >> len;
        cable_list.push_back(len);
    }
    cable_master cm(cable_list, K, 0.0001);
    std::cout << std::fixed << std::setprecision(2) << floor(cm.search() * 100) / 100.0 << std::endl;
}

2次元累積和

用途

  • \displaystyle \sum_{i = p_x}^{q_x} \sum_{j = p_y}^{q_y}f(x_i, y_j)の計算

アルゴリズムのポイント

  1. あらかじめ、全てのn,mについて\displaystyle  S_{nm} = \sum_{i = 1}^{n} \sum_{j = 1}^{m} f(x_i, y_j)  を求めておく
    • これは、\displaystyle S_{nm} =f(x_n, y_m) + S_{n-1, m} + S_{n, m-1} - S_{n-1, m - 1}で効率的に計算できる
  2. 目的の値は\displaystyle  \sum_{i = p_x}^{q_x} \sum_{j = p_y}^{q_y}f(x_i, y_j) =  S_{q_x, q_y} - S_{q_x-1, q_y} - S_{q_x, q_y-1} + S_{q_x - 1, q_y-1}  で求まる

抽象クラスの実装

このアルゴリズムはストラテジーパターンで実装する。具体的には、f(x_i, y_j)を計算する関数オブジェクト受け取り、それに応じた部分和を求める。

#include <vector>
#include <functional>

template <typename T>
class func_sum
{
private:
    std::vector<std::vector<T>> sum_table;
protected:
    virtual T acc_sum(int i, int j);
public:
    func_sum(std::function<T(int, int)> func, int nx, int ny) ;
    ~func_sum();
    //目的の区間挿話を求める関数
    T calc_section_sum(int px, int qx, int py, int qy);
};

template <typename T>
T func_sum<T>::acc_sum(int i, int j) {
    return sum_table[i][j];
};

template <typename T>
func_sum<T>::~func_sum() {}

template <typename T>
func_sum<T>::func_sum(std::function<T(int, int)> func, int nx, int ny) {
    sum_table = std::vector<std::vector<T>>(nx + 1, std::vector<T>(ny + 1, T()));
    //j=0の場合の総和の計算
    for (int i = 1; i <= nx; i++) {
        sum_table[i][0] += func(i - 1, 0);
    }
    //i=0の場合の総和の計算
    for (int j = 1; j <= ny; j++) {
        sum_table[0][j] += func(0, j - 1);
    }
    //総和の計算
    for (int i = 1; i <= nx; i++) {
        for (int j = 1; j <= ny; j++) {
            sum_table[i][j] += func(i - 1, j - 1) +
                sum_table[i - 1][j] +
                sum_table[i][j - 1] -
                sum_table[i - 1][j - 1];

        }
    }
}

template <typename T>
T func_sum<T>::calc_section_sum(int px, int qx, int py, int qy) {
    return acc_sum(qx, qy) - 
           acc_sum(qx, py - 1) - 
           acc_sum(px - 1, qy) + 
           acc_sum(px - 1, py - 1);
}

具体的な問題を解く実装

\sum_{i = p_x}^{q_x} \sum_{j = p_y}^{q_y} \frac{x_i}{y_i}の計算

f(x_i, y_j)= \frac{x_i}{y_i}  として、上記アルゴリズムに関数オブジェクトを渡す。

int main()
{
    int N, M;
    int px, py, qx, qy;
    std::cin >> N >> M ;
    std::cin >> px >> qx >> py >> qy;
    std::vector<double> x_list(N, 0);
    std::vector<double> y_list(M, 0);
    for (int i = 0; i < N; i++) {
        double x;
        std::cin >> x;
        x_list[i] = x;
    }
    for (int i = 0; i < M; i++) {
        double y;
        std::cin >> y;
        y_list[i] = y;
    }
    auto func = [](int i, int j, 
                   const std::vector<double>& x_list,
                   const std::vector<double>& y_list) { return x_list[i] / y_list[j]; };

    std::function<double(int, int)>  func_curried = std::bind(func, 
                                                              std::placeholders::_1, 
                                                              std::placeholders::_2, 
                                                              x_list, 
                                                              y_list);

    func_sum<double> ds(func_curried, x_list.size(), y_list.size() );
    std::cout<<ds.calc_section_sum(px, qx, py, qy)<<std::endl;
    return 0;
}
ABC106 D:AtCoder Express 2

問題はリンクを参照。この問題はf(x_i, y_j)  x_i = L_i, y_i = R_iのときそのペア(L_i, R_i)の存在する個数を返す関数とすれば2次元累積和で解ける。

#include <iostream>
int main()
{
    int N, M, Q;
    std::cin >> N >> M >> Q;
    std::vector<std::vector<int>> table(N, std::vector<int>(N, 0));
    for (int i = 0; i < M; i++) {
        int L, R;
        std::cin >> L >> R;
        table[L - 1][R - 1]++;
    }

    auto func = [](int i, int j,
        const std::vector<std::vector<int>>& table) { return table[i][j]; };

    std::function<int(int, int)>  func_curried = std::bind(func,
        std::placeholders::_1,
        std::placeholders::_2,
        table);

    func_sum<int> ds(func_curried, N, N);

    for (int i = 0; i < Q; i++) {
        int p, q;
        std::cin >> p >> q;
        std::cout << ds.calc_section_sum(p, q, p, q) << std::endl;
    }
    return 0;
}

しゃくとり法

用途

  • 解きたい問題が以下のいずれかの性質を満たす時、条件C(l, r)を満たすl,rの数え上げ、もしくは評価関数f(l,r)の最大化・最小化
    • l_p \leq l\leq r \leq r_qとして、C(l, r)が条件を満たすならばC(l_p, r_q)も条件を満たす(以降、これを上位区間条件と呼ぶ)
      • 例えば、整数の数列x_1, x_2, \cdots, x_nの中からその和がS以上となる連続する部分数列を数え上げる問題
    • l \leq l_p \leq r_q \leq rとして、C(l, r)が条件を満たすならばC(l_p, r_q)も条件を満たす(以降、これを部分区間条件と呼ぶ)
      • 例えば、整数の数列x_1, x_2, \cdots, x_nの中からその和がS以下となる連続する部分数列を数え上げる問題

アルゴリズム

上位区間条件を満たす問題と部分区間条件を満たす問題とで若干異なる(うまいことやれば、共通にできるかもしれない)

  1. 条件を満たさない間rをインクリメント(右端まで来てインクリメントできない場合はそのまま)
  2. 条件を満たしたところで、インクリメントを止め、その時点のl, rを記録(この時点でlp = l, r \leq r_pを満たす(l_p, r_p)は全て条件を満たす)
  3. 条件が満たさなくなるまでlをインクリメント
  4. この1~3の操作をrがインクリメントできなくなり、かつ、条件C(l, r)を満たさなくなるまで繰り返す
  1. 条件を満たす間rをインクリメント(右端まで来てインクリメントできない場合はそのまま)
  2. 条件を満たさなくなったところで、インクリメントを止め、その時点のl, rを記録(この時点でlp = l, r_p \leq rを満たす(l_p, r_p)は全て条件を満たす)
  3. 条件を満たすまでlをインクリメント
  4. この1~3の操作をrがインクリメントできなくなり、かつ、lrと同じになるまで繰り返す

抽象クラスの実装

こちらもテンプレートメソッドパターンで実装する。上述の通り上位区間条件と部分区間条件を分けて考えるので、それぞれ抽象クラスを作る。

#include <functional>
 
//ある区間が条件を満たすならば、
//それを包含する上位区間も条件を満たす問題を解くクラス
template <typename T>
class abstract_inchworm_superset
{
private:
    virtual void move_left_pointer() = 0;//lを右に動かす処理
    virtual void move_right_pointer() = 0;//rを右に動かす処理
    virtual bool is_cond() = 0; //条件成立か否かを返す
    virtual bool is_right_end() = 0;//rが右端に到達したかを返す
    virtual std::pair<T, T> no_next() = 0;//次の要素がない場合に終了値を返す
    virtual std::pair<T, T> get_LR() = 0;//lとrを返す
public:
    virtual std::pair<T, T> next(); //条件を満たすデータを返し、ポインタを次に進める
};
 
template <typename T>
std::pair<T, T> abstract_inchworm_superset<T>::next(){
    while (is_cond() == false && is_right_end() == false) {
            move_right_pointer();
    }
    if (is_cond() == false)
        return no_next();
    std::pair<T, T> ret = get_LR();
    move_left_pointer();
    return ret;
}
 
//ある区間が条件を満たすならば、
//それが包含する部分区間も条件を満たす問題を解くクラス
template <typename T>
class abstract_inchworm_subset
{
private:
    virtual void move_left_pointer() = 0;//lを右にインクリメントする処理
    virtual void move_right_pointer() = 0;//rを右にインクリメントする処理
    virtual bool is_cond() = 0; //条件成立か否かを返す
    virtual bool is_right_end() = 0;//rが右端に到達したかを返す
    virtual std::pair<T, T> no_next() = 0;//次の要素がない場合に終了値を返す
    virtual std::pair<T, T> get_LR() = 0;//lとrを返す
    virtual bool is_left_eq_right() = 0;//lとrが一致しているかを返す
public:
    virtual std::pair<T, T> next(); //条件を満たすデータを返し、ポインタを次に進める
};
 
template <typename T>
std::pair<T, T> abstract_inchworm_subset<T>::next() {
    std::pair<T, T> ret;
    while (is_right_end() == false && is_cond() == true) {
        move_right_pointer();
    }
    ret = get_LR();
    if (is_left_eq_right() == true) {
        return no_next();
    }
    
    move_left_pointer();
    return ret;
}

具体的な問題を解く実装

POJ 3061:Subsequence

蟻本に例として載っている和がS以上になる部分数列を数え上げる問題。上位区間条件型のしゃくとり法で解ける。

#include <iostream>
#include <algorithm>
#include <vector>

class subsequence :
    public abstract_inchworm_superset<int>
{
private: 
    int Left, Right;
    long long S;
    int N;
    long long LRsum;
    std::vector<long long> sequence;
    virtual void move_left_pointer();
    virtual void move_right_pointer();
    virtual bool is_cond();
    virtual bool is_right_end();
    virtual std::pair<int, int> no_next();
    virtual std::pair<int, int> get_LR();
public:

    void set_S(long long S);
    subsequence(std::vector<long long>& sequence, int N, int start);
};

subsequence::subsequence(std::vector<long long>& sequence, int N, int start):
    sequence(sequence),
    N(N),
    Left(start),
    Right(start),
    S(0),
    LRsum(0)
{
}

void subsequence::set_S(long long S) {
    this->S = S;
}

void subsequence::move_left_pointer() {
    LRsum -= sequence[Left];
    Left++;
}

void subsequence::move_right_pointer() {
    if (Right < N) {
        LRsum += sequence[Right];
        Right++;
    }
}

bool subsequence::is_cond() {
    if (LRsum >= S)
        return true;
    return false;
}

bool subsequence::is_right_end() {
    if (Right >= N)
        return true;
    return false;
}

std::pair<int, int> subsequence::no_next() {
    return std::make_pair(-1, -1);
}

std::pair<int, int> subsequence::get_LR() {
    return std::make_pair(Left, Right);
}

int main()
{

    int N;
    std::cin >> N;
    for (int k = 0; k < N; k++) {
        std::vector<long long> sequence;
        int n;
        long long S;
        std::cin >> n >> S;
        for (int i = 0; i < n; i++) {
            long long a;
            std::cin >> a;
            sequence.push_back(a);
        }

        subsequence seq(sequence, sequence.size(), 0);
        seq.set_S(S);
        int min_len = INT_MAX;
        while (1) {
            std::pair<int, int> LR = seq.next();
            if (LR.first == -1) break;
            min_len = std::min(min_len, LR.second - LR.first);
        }
        if (min_len == INT_MAX)
            std::cout << 0 << std::endl;
        else
            std::cout << min_len << std::endl;
    }
    return 0;
}
ARC098 D:Xor Sum 2

部分区間条件型のしゃくとり法で解けると気付くのが難しい問題。気づいてしまえば、普通にしゃくとり法を適用するだけ(部分区間を数え上げる問題なので、何も考えずしゃくとり法を適用してACしちゃったって人も結構いるかも)。なぜ、しゃくとり法で解けるかは解説参照。
ここで注意するのは、関数is_cond()のtrueを返す条件が現在指しているrの更に一つ先のrにおいて、条件が成立している場合としている。つまり、right pointer(r)を動かしても条件成立するなら、実際にright pointerを動かすという処理になっている。なお、この実装では区間[l, r)という半区間で扱っているため、一つ先のrrが指示しているインデックスの値を足すこと(XORを取ること)に当たる。

#include <vector>
#include <stdio.h>
class xor_sum :
    public abstract_inchworm_subset<int>
{
private:
    int Left, Right;
    
    int N;
    long long LR_XOR;
    long long LRsum;
    std::vector<long long> sequence;
 
    virtual void move_left_pointer();
    virtual void move_right_pointer();
    virtual bool is_cond();
    virtual bool is_right_end();
    virtual std::pair<int, int> no_next();
    virtual std::pair<int, int> get_LR();
    virtual bool is_left_eq_right();
public:
 
    xor_sum(std::vector<long long>& sequence, int N, int start);
};
 
xor_sum::xor_sum(std::vector<long long>& sequence, int N, int start) :
    sequence(sequence),
    N(N),
    Left(start),
    Right(start),
    LR_XOR(0),
    LRsum(0)
{
}
 
 
void xor_sum::move_left_pointer() {
    LRsum -= sequence[Left];
    LR_XOR ^= sequence[Left];
    Left++;
 
}
 
void xor_sum::move_right_pointer() {
    if (Right < N) {
        LRsum += sequence[Right];
        LR_XOR ^= sequence[Right];
        Right++;
    }
}
 
bool xor_sum::is_cond() {
    if (Left >= N) 
        return false;
    //Rightが1つ先に行っても条件を満たすか
    if ((LRsum + sequence[Right]) == (LR_XOR ^ sequence[Right]))
        return true;
    return false;
}
 
bool xor_sum::is_right_end() {
    if (Right >= N)
        return true;
    return false;
}
 
std::pair<int, int> xor_sum::no_next() {
    return std::make_pair(-1, -1);
}
 
std::pair<int, int> xor_sum::get_LR() {
    return std::make_pair(Left, Right);
}
 
bool xor_sum::is_left_eq_right() {
    return (Left == Right);
}
 
 
int main()
{
    int N;
    std::vector<long long> A;
    scanf("%d", &N);
 
    for (int i = 0; i < N; i++) {
        long long a;
        scanf("%lld", &a);
        A.push_back(a);
    }
    xor_sum xs(A, N, 0);
    long long count = 0;
    while (1) {
        std::pair<int, int> LR = xs.next();
        if (LR.first == -1) break;
        count += LR.second - LR.first ;
    }
    printf("%lld\n", count);
    return 0;
}

おわりに

今回は勉強のため、本質的なアルゴリズムを記述した抽象クラスと実際の問題を解く具象クラスを分けましたが、実際のコンテスト時はそんなことせずに、普通に書いた方が速いと思います。

競プロとかに使うアルゴリズム実装メモ(幅優先・深さ優先探索、union-find、最小全域木)

はじめに

以前の記事にて最短経路問題を解くアルゴリズムの実装を書きましたが、今回はその続きとしてグラフアルゴリズムの中でも幅優先探索深さ優先探索、union-find、最小全域木問題を解くアルゴリズム2種について実装を書いていきます。
例によって、自分の脳内整理メモの側面が強い上、検索すればそれぞれのアルゴリズムについて素晴らしい解説記事が出てくるので、まじめな解説を見たい場合はそちらを当たった方が良いと思います。

幅優先探索

用途

  • コストが全て同じグラフ(重みなしグラフ)の最短経路探索。
  • 連結成分の列挙。
  • 重みなしグラフの最小全域木

アルゴリズムのポイント

  • 隣接する頂点から順々に探索していく。
  • 距離を求める場合は、遷移元の距離に+1をする。

実装のポイント

  • キューを使う。
  • まずスタートの頂点をキューに入れ、キューの先頭から頂点を取り出し、それに隣接する頂点を探索。
  • 未探索の全て頂点はキューに追加して、またキューの先頭から頂点を取り出して・・・ということを目的の頂点が見つかるか、キューが空になるまで繰り返していく。
  • 既に探索済みの頂点かというフラグを持たせておき、既に見た頂点だった場合はキューに追加しない(下記実装では距離が初期値-1でなければ探索済みと判定)。

実装

以前と同様、メイン関数は省略。

#include <unordered_map>
#include <vector>
#include <queue>

class breadth_first_search
{
private:
    int node_num;
    int abs_time;
    std::unordered_map<int, std::vector<int>> adj_list;
    std::vector<int> dists;//始点からの距離

public:
    breadth_first_search(const int node_num,
        const std::unordered_map
        <int, std::vector<int>>& adj_list);
    void exec_search(int start);
    int get_dist(int v);
};

breadth_first_search::breadth_first_search(const int node_num,
    const std::unordered_map
    <int, std::vector<int>>& adj_list) :
    node_num(node_num),
    adj_list(adj_list),
    dists(node_num + 1, -1)
{
}

void breadth_first_search::exec_search(int start) {
    std::queue<int> v_queue;
    v_queue.push(start);
    dists[start] = 0;
    int d = 0;
    while (1) {
        if (v_queue.empty() == true)
            break;
        int v = v_queue.front();
        v_queue.pop();
       
        for (int &adj : adj_list[v]) {
            if (dists[adj] == -1) {
                dists[adj] = dists[v] + 1;
                v_queue.push(adj);
            }
        }
    }
}

int breadth_first_search::get_dist(int v) {
    return dists[v];
}

AOJのALDS1_11_Cが通ることを確認。

深さ優先探索再帰呼び出し版)

用途

  • コストが全て同じグラフ(重みなしグラフ)において、与えられたコスト以内でたどり着ける頂点の列挙。
  • 現実における迷路探索(現実だと幅優先探索は分身かワープでもできないと無理)。
  • 深さを時系列的な情報と関連づけて探索するケース。

(基本、幅優先探索と同じことができるが、それぞれ解きたい問題に応じて使い分ける)

アルゴリズムのポイント

  • 行けるところまでどんどん深くまで探索していく。
  • これ以上探索できないというところまで探索したら根まで戻り、更に別なルートの探索を進める。

実装のポイント

  • 幅優先探索と同じノリでスタックを使うパターンと再帰呼び出しを使うパターンがある(今回の実装は再帰版)。
  • 再帰呼び出しで行う場合、分岐して最深部まで探索した後、分岐点まで戻って来たときに既に探索済みか否かという情報をクリアしたいケースの実装が楽になる(上記の深さと時系列情報を関連付ける場合、分岐点に戻る=時系列的にも戻る=探索済み情報をクリアする、とするべき問題もある。これとか)。
  • 今回の実装では、AOJの問題に対応するべく、探索で移動したらカウントアップ(後戻り時もカウントアップ)するabs_time変数を用意し、exec_search_sub関数に入った時と出る時にカウントアップさせている。また、find_abs_timesには各頂点を見つけたときのabs_timeを保持しておく。
  • find_abs_timeが初期値(-1)でなければ、確定頂点として探索しない。
  • このタイムスタンプのつけ方は、問題に応じて適宜カウントアップするタイミング、回数、場合によっては初期値にクリアする等の変更が必要。

実装

#include <unordered_map>
#include <vector>
class depth_first_search
{
private:
    int node_num;
    int abs_time;
    std::unordered_map<int, std::vector<int>> adj_list;
    std::vector<int> depths;//ルートからの深さ(距離ではない)
    std::vector<int> find_abs_times;//頂点に到達した時刻
    std::vector<int> finish_abs_times;//サーチが完了した時刻
    
    void exec_search_sub(int start, int depth);

public:
    depth_first_search(const int node_num, 
                       const std::unordered_map
                        <int, std::vector<int>>& adj_list);
    
    void exec_search(int start);
    int get_depth(int v);
    int get_find_abs_time(int v);
    int get_finish_abs_time(int v);
};

depth_first_search::depth_first_search(const int node_num, 
                                       const std::unordered_map
                                        <int, std::vector<int>>& adj_list):
    node_num(node_num),
    abs_time(0),
    adj_list(adj_list),
    depths(node_num + 1, -1),
    find_abs_times(node_num + 1, -1),
    finish_abs_times(node_num + 1, -1)
{
}

void depth_first_search::exec_search_sub(int start, int depth) {
    
    if (find_abs_times[start] != -1) return;
    depths[start] = depth;
    abs_time++;//探索進みのカウントアップ
    find_abs_times[start] = abs_time;
    for (auto &v : adj_list[start]) {
        exec_search_sub(v, depth + 1);         
    }
    abs_time++;//探索後戻りのカウントアップ
    finish_abs_times[start] = abs_time;
}

void depth_first_search::exec_search(int start) {
    for (int i = start; i < node_num + 1; i++) {
        exec_search_sub(i, 0);
    }
}

int depth_first_search::get_depth(int v) {
    return depths[v];
}
int depth_first_search::get_find_abs_time(int v) {
    return find_abs_times[v];
}

int depth_first_search::get_finish_abs_time(int v) {
    return finish_abs_times[v];
}

AOJのALDS1_11_Bが通ることを確認。

union-find

用途

  • 集合を扱うデータ構造。
  • 値xとyが同じ集合に属するか否かの判定や2つの集合の結合が高速で実行可能。
  • 後に説明する最小全域木問題を解くアルゴリズムであるクラスカル法でも使われる。

アルゴリズムのポイント

  • 各集合、木構造でデータを保持する。
  • 値xとyが同じ集合か否かを判断するには(以下、find操作)、木の根を見て同じだったら同じ集合、異なれば別の集合に属すると判断すればよい。
  • そのため、なるべく木の高さが低くなるようにデータを保持した方が、高速に判断できるようになる。
  • このことから、2つの集合を結合する操作(以下、unite操作)では木の高さが低い方の根の親を高い方の根になるように連結させる。こうすれば、連結後の木の高さは元の集合のうち高い方の高さと同じになるので、高さは増えない(両方の木の高さが同じ場合は、木の高さが1だけ増える)。

実装のポイント

  • 木は親ノードの値を保持する1次元配列と所属する木の高さを保持する1次元配列で表現する。
  • こうすれば、unite操作時、低い方の木の根の親を書きかえるように親を保持する配列を操作だけで良くなる。
  • find操作時に全てのノードを直接親ノードとつなぐことで木の縮約を行う。
  • 木の縮約時に木の高さの更新を行わないので実際の木の高さと配列に保持している高さが合わなくなるが、あまり気にしなくてよい(重要なのはノードの親子関係であり、それさえ合っていればunion-findは問題なく動作する。木の高さの配列に入っている値がおかしいとunite操作時に、接続する関係が逆になり木が高くなる恐れがあるが、それよりもfind操作で縮約する方がメリットが大きい)。

実装

#include <vector>

class union_find
{
private:
    std::vector<int> parent;
    std::vector<int> rank;
    int node_num;
public:
    union_find(int node_num);
    void unite(int x, int y);
    int find(int x);
};

union_find::union_find(int node_num):
    parent(node_num + 1, 0),
    rank(node_num + 1, 0)
{
    for (int i = 0; i <= node_num; i++) {
        parent[i] = i;
    }
}

void union_find::unite(int x, int y) {
    int x_root = find(x);
    int y_root = find(y);
    //既に同じ根なら、同じ所属なので何もしない
    if (x_root == y_root) return;
    //ランクが大きい方に小さい方をくっつける
    //結合後のランクは大きい方になる。
    if (rank[x_root] > rank[y_root]) {
        parent[y_root] = x_root;
    }
    else {
        parent[x_root] = y_root;
        if(rank[x_root] == rank[y_root]){
            //ランクが等しい場合は、必ずランクが1上がる。
            rank[y_root]++;
        }
    }
}

int union_find::find(int x) {
    if (x != parent[x]) {
        // findの結果を各parentに入れることで経路を圧縮。
        //ただし、変更コストが大きいのでランクは変えてない。
        parent[x] = find(parent[x]);
    }
    return parent[x];
}

後に説明するクラスカル法の中で一緒に確認。

プリム法

用途

  • 最小全域木問題を解く。
  • 最小全域木問題自体は与えられた無向グラフから、ループがなくなるように枝を取り除き、辺のコストの合計が最小になる全域木を求める問題。
  • 例えば、全ての街間で行き来できるように最小コストで道を作りたい的な問題に使える(完全グラフを構成して最小全域木を求める等*1)。

アルゴリズムのポイント

  • ダイクストラ法と似ており、まず探索の初期位置となる頂点を決め、そこから出ている辺を探索中リストに追加する。
  • 探索中リストの中で最小のコストの辺を取り出す。この辺は最小全域木の構成要素として確定する。
  • 上記で取り出した辺につながっている頂点を次の探索対象の頂点とし、その頂点を確定頂点とする。また同様にそこから出ている辺を探索中リストに追加する(ただし、確定頂点へつながる辺は登録しない)。
  • これらの操作を全ての頂点が確定するまで繰り返す。

実装のポイント

  • ダイクストラ法と同様に優先度付きキューを使う。
  • 各頂点に対して現在どのコストで優先度付きキューに登録しているかを保持する外部変数keyを用意する。
  • 優先度付きキューから最小コストの辺(と接続先の頂点)を取り出し、この辺を最小全域木の構成要素として確定させる。
  • 同時に接続先の頂点も確定頂点として記録し、次はその頂点につながっている辺を見ていく。
  • 辺のコストを見て、接続頂点のkeyの値より小さいコストの辺だったらkeyの値を更新し、優先度付きキューに追加する。もし、大きいコストの辺だったら、その辺は最小全域木の構成要素にはなり得ないので、優先度付きキューには追加しない。これ結構大事。これをやらないと、どんどん優先度つきキューに登録されていくので、サイズがばかでかくなり、確定ノードかの判定回数が大幅に増えるので遅くなる可能性がある。

実装

#include <vector>
#include <unordered_map>
#include <algorithm>
#include <vector>
#include <unordered_map>
#include <queue>
#include <functional>
#include <climits>


class edge {
public:
    edge(int to, long long cost);
    int get_to() const;
    long long get_cost() const;
private:
    int to;
    long long cost;
};

edge::edge(int to, long long cost) : to(to), cost(cost) {}

int edge::get_to() const { return to; }

long long edge::get_cost() const { return cost; }

class prim
{
private:
    //unordered_mapを使って隣接リストを構築
    //first:接続元、second:辺のコストと接続先
    std::unordered_map<int, std::vector<edge>> adj_list;
    std::vector<int> pred;
    std::vector<long long> key;
    int node_num;
    std::vector<bool> is_done;//確定ノードか否かを保持
    long long MST_cost;
public:
    prim(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list);
    long long get_MST_cost() const;
    //std::unordered_map<int, std::vector<edge>> calc_MST();
    std::unordered_map<int, std::vector<edge>> get_MST_list();
    void calc_MST();
};

prim::prim(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list):
    adj_list(adj_list),
    node_num(node_num),
    is_done(node_num + 1, false),
    pred(node_num + 1, INT_MAX),
    key(node_num + 1, LLONG_MAX),
    MST_cost(0)
{

}

void prim::calc_MST() {
    MST_cost = 0;
    key = std::vector<long long>(node_num + 1, LLONG_MAX);
    pred = std::vector<int>(node_num + 1, INT_MAX);
    //優先度付きキューpairで頂点番号と最短距離を保持
    //firstの要素で比較されるので、firstが距離、secondを遷移先の頂点とする
    std::priority_queue<std::pair<long long, int>,
        std::vector<std::pair<long long, int>>,
        std::greater<std::pair<long long, int>>> p_queue;
    p_queue.push(std::make_pair(0, adj_list.begin()->first));
    key[adj_list.begin()->first] = 0;
    pred[adj_list.begin()->first] = -1;

    while (1) {
        if (p_queue.empty() == true) break; //キューが空なら終了
        std::pair<long long, int> cost_node = p_queue.top();
        p_queue.pop();
        if (is_done[cost_node.second] == false) {
            MST_cost += cost_node.first;
        }
        
        int from_node = cost_node.second;
        is_done[from_node] = true; //キューから取り出した頂点を確定頂点とする
        for (auto &it : adj_list[from_node]) {
            int adj_node = it.get_to();
            if (is_done[adj_node] == false) {
                long long adj_cost = it.get_cost();
                if (adj_cost < key[adj_node]) {
                    p_queue.push(std::make_pair(adj_cost,adj_node));
                    key[adj_node] = adj_cost;
                    pred[adj_node] = from_node;
                }
            }
        }
    }
}


long long prim::get_MST_cost() const {
    return  MST_cost;
}

std::unordered_map<int, std::vector<edge>> prim::get_MST_list() {
    std::unordered_map<int, std::vector<edge>> MST_list;
    for (int i = 0; i < node_num + 1; i++) {
        if (pred[i] >= 0 && pred[i] < node_num + 1) {
            MST_list[i].push_back(edge(pred[i], key[i]));
            MST_list[pred[i]].push_back(edge(i, key[i]));
        }
    }
    return MST_list;
}

AOJのGRL_2_Aが通ることを確認。

クラスカル

用途

アルゴリズムのポイント

  • 最初は全ての頂点間が未連結の状態(辺がない状態)からスタート。
  • 辺をコストでソート後、コストが小さい辺から順に全域木の構成要素に追加していく。ただし、ループができる場合はその辺は追加しない。
  • 辺を追加してループができるか否かは、その辺で結んでいる両頂点が既に連結関係にあるか否かで判断可能(既に連結関係にあれば、その辺を追加すると必ずループができる)。

実装のポイント

  • 連結されている頂点同士が同じ集合になるようunion-findで管理し、頂点間の連結関係を効率よく判定および操作できるようにする。
  • 実装では隣接リストをコンストラクタで渡しているが、辺を管理する変数を用意する方が楽なので、隣接リストから辺のリストに変換させている(最初から辺のリストを構築してもよい)。
  • 辺はstd::pair<int, std::pair<int, int>>型で(コスト, (頂点, 頂点))という変数にし、リストにしたとき第1要素のコストでソートされるようにする。

実装

上述のunion-findを使用する。

#include <algorithm>
#include <vector>
#include <unordered_map>
#include <queue>
#include <functional>
#include <climits>

class edge {
public:
    edge(int to, long long cost);
    int get_to() const;
    long long get_cost() const;
private:
    int to;
    long long cost;
};

edge::edge(int to, long long cost) : to(to), cost(cost) {}

int edge::get_to() const { return to; }

long long edge::get_cost() const { return cost; }

class kruskal
{
private:
    //unordered_mapを使って隣接リストを構築
    //first:接続元、second:辺のコストと接続先
    std::unordered_map<int, std::vector<edge>> adj_list;
    //最小全域木の隣接リスト
    std::unordered_map<int, std::vector<edge>> MST_list;
    int node_num;
    long long MST_cost;
public:
    kruskal(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list);
    long long get_MST_cost() const;
    std::unordered_map<int, std::vector<edge>> get_MST_list();
    void calc_MST();
};


kruskal::kruskal(const int node_num, const std::unordered_map<int, std::vector<edge>>& adj_list):
    adj_list(adj_list),
    node_num(node_num),
    MST_cost(0)
{
}

void kruskal::calc_MST() {
    //辺のリスト。firstがコスト、secondがつながっている両頂点
    std::vector<std::pair<long long, std::pair<int, int>>> edge_list;

    //隣接リストから辺のリストへ変換
    for (auto &v : adj_list) {
        for (auto &e : v.second) {
            edge_list.push_back(std::make_pair(e.get_cost(), 
                        std::make_pair(v.first, e.get_to())));
        }
    }

    std::sort(edge_list.begin(), edge_list.end());

    union_find uf(node_num + 1);
    MST_cost = 0;
    MST_list.clear();
    for (auto &ed : edge_list) {
        int v1 = ed.second.first;
        int v2 = ed.second.second;
        long long cost = ed.first;
        //同じ根を持つなら、既に両頂点は同じ木の中
        //=辺を追加するとループができるのでなにもしない。
        if (uf.find(v1) != uf.find(v2)) {
            //根が異なる場合、両頂点間に辺を追加
            uf.unite(v1, v2);
            MST_cost += cost;
            MST_list[v1].push_back(edge(v2, cost));
            MST_list[v2].push_back(edge(v1, cost));
        }
    }
}

long long kruskal::get_MST_cost() const {
    return MST_cost;
}

std::unordered_map<int, std::vector<edge>> kruskal::get_MST_list() {
    return MST_list;
}

プリム法同様、AOJのGRL_2_Aで確認。

おわりに

前回、今回とグラフアルゴリズムについて見てきましたが、これら以外にもネットワークフロー問題やマッチング問題など様々な問題とそれらを解くアルゴリズムがあります。また、前回紹介した最短経路問題や今回扱った最小全域木問題、他にもマッチング問題など多くの問題は離散凸解析の枠組みで抽象化され、その枠組みの中で最適化アルゴリズムが提案されています。この離散凸解析についてはまだ全然理解していませんが、いつか勉強して記事を書きたいところです。
あと、このシリーズはもう少し続きます。次回は尺取り法や累積和法などを書いていこうと思います。

*1:競プロだと、完全グラフを考えると大抵の場合TLEになるので問題の性質から元のグラフをどのように構築するかというところが肝になる。

競プロ関係の雑メモ 2018/9/9

はじめに

以前と同様、atcoderの400~500点問題穴埋め時のメモに加え、ABC109参加のメモです。

ARC088-D Wide Flip

区間[l - 1, r]を反転させた後、区間[l, r]を反転させると、l - 1番目の1文字のみ反転させることができる。同様に区間[l, r + 1]を反転させた後、区間[l, r]を反転させると、r + 1番目の1文字のみ反転させることができる。ということは、選べる区間の最小の長さKが与えられたとして、左端からN-K個、右端からもN-K個、合計2(N-K)個の要素は自由に0,1を選べる、ということになる。
じゃあ、左端もしくは右端から何番目の要素まで自由に変えられるようにする必要があるかと考えると、真ん中の要素がいくつ連続で同じ文字かを見ればよいという結論に至る。つまり、中央N/2番目の要素から左右にC個が全て同じ要素だとすると、そのC個は変化させる必要がなく、左端から(N-C)/2番目までの要素もしくは右端から(N-C)/2番目までの要素のどれかが自由に0, 1を選べるようKを定めればよい。よって、N - K = (N - C)/2という方程式が出てくるので、これをKについて解いてK = N - (N - C)/2となる。
Cについては、普通に真ん中から見て同じ要素をカウントすればいいのだけど、文字列の長さが偶数か奇数かで話が変わってくるので、少し注意が必要(実際、1度for文回す回数間違えてWAになった)。

ARC064-D An Ordinary Game

まず個数に着目。1個ずつ文字を削除していくゲームなので、負けるときの文字数の法則性を見つけ出せれば、後は逆算して元の文字数から勝者がわかると考えた。じゃあ、負けるときの文字列とは?ってところについては、2パターン考えられる。ひとつはabのような文字数が偶数での負け、もうひとつはabaのような文字数が奇数での負け(この辺り、大体の予想は立ったが、じゃあ本当?っていうところが解説見るまで分からなかった)。では、どういうときに偶数負けパターンになるか?奇数負けパターンになるか?
ポイントはゲームの性質上、先頭の文字と後尾の文字は必ず最後まで残るということ。このことから、先頭の文字と後尾の文字が同じであれば奇数負けパターン、異なるのであれば偶数負けパターンになる。なぜなら、例えばa**b***a*b(*はaとb以外の文字)のように先頭と後尾が異なるなら、必ず*は取り除くことができ、最終的な負けパターンとしてはabかababになる(そうなるように、両プレイヤーがコントロール可能)。どちらになるかは、一意に定まらないが偶数の文字数になることはかわらない。これと同じような考えで、a*b**a*b**a等のように、先頭と後尾が同じ場合は、abaかa*aかababaかのいずれかが最終形になるので必ず奇数の文字数になる。ここまで来れば、上述の通り逆算して元の文字列数からどちらが勝つかはわかる。
どうもゲーム系は苦手意識があったが(未だに「互いが最適に行動を取る」ってことのイメージがつかめてない)、最終的に詰み(負け)となる形を考えて、そこから逆算していくっていう方法が有効そう。
この問題では一発でACを取れたが、根拠の薄い予想がたまたま当たったっていうパターンなので、あまりよろしくない*1

ARC067-D Walk and Teleport

まず、いくつか移動パターンを考える。

  1. 西から東へ順々に移動するパターン(テレポートしても町は飛び越えない)。
  2. テレポートを使って、まだ移動していない町を飛び越え、その飛び越えた町へ行くため東から西へ逆走をするパターン。
  3. 全てテレポートで移動するパターン。

この中で2の移動パターンは使う必要がない。なぜなら、最終的には全ての町を訪れる必要があるため、逆走しても順方向に移動しても距離は変わらず、町を飛び越すメリットがないから。また、1のパターンはB<A(X[i + 1] - X[i])であればテレポートを使うし、そうでなければ徒歩で移動する。3のパターンは1の条件で、全ての町間の距離がテレポート使用条件に合うなら自動的に適用される。ということで、min(A(X[i + 1] - X[i]), B)を足しこんでいけばよい。
500点問題の中では明らかに簡単なような・・・

ABC 109参加メモ

2018/9/8開催のABC109に参加。4完達成もD問題の実装に時間がかかりすぎたのがよろしくない。

A問題

A*Bが奇数か偶数か。

B問題

setを使って既出かの判定と、前回の文字列の最後の文字と今見てる文字列の頭文字が同じかという判定を行う。

C問題

Xも数列に組み込んでソートしてx[i + 1] - x[i]の最大公約数を求めたが、解説読むとソートなんかしなくても|X - x[i]|の最大公約数でよかった。

D問題

着目しているマスのコイン数が奇数なら隣のマスにコインを1枚ずらし、着目するマスもそのずらした先に移動。逆に、着目しているマスのコイン数が偶数なら何もせず、着目するマスのみ隣にずらす。この操作を左端からジグザグでもグルグルでも何でもいいから順々に全てのマスに対して行えば、目的が達成される。
その理由としては、まず奇数のマスから奇数のマスへコインを一つずらすと両方偶数になる。その上で、ずらし元が奇数、ずらし先が偶数の場合、操作後はずらし先が奇数になるが、直後に奇数になったずらし先に対して操作を行えば、そのマスも偶数になる。この操作を続けていくと、どこかかでずらし元、ずらし先ともに奇数という個所が出てきて、両方偶数になる。もしそのような個所が出てこない場合、最後のマスまでコインが伝搬され最後のマスのみが奇数になる。従って、最適な最終形としては、コインの総枚数が偶数の場合は全てのマスが偶数になり、コインの総枚数が奇数の場合、1つのマスのみが奇数でそれ以外全手のマスが偶数になる。
ただ、いざ実装ってなると、配列のジグザグ操作またはグルグル操作の実装が面倒だなーと億劫になり、結局、既に操作済みのマスか?という情報を記録する配列を別途用意。ずらせるマスを探し、ずらした先の座標を次のずらし元とする作戦を取る。ACだったが実装に時間がかかりすぎた。

おわりに

最近、競プロ関連チラ裏ブログと化しているので、次回はもう少しまともな内容を書きます。

*1:もう少し正確に言うと、予想した負けパターンは、abのような2文字か、abaのような3文字、もしくは最初からabababのようなFirstプレイヤーが取れず負けるというパターンの3種類を考え、最後のパターンだけ(無駄な)例外処理をする実装でACになった。

競プロ関係の雑メモ 2018/9/1

はじめに

競プロの問題解いてて、解けなかった問題は何故解けなかったのか、何故その発想が出なかったのかを記していく自分用メモです。atcoderの400、500点くらいの問題なら、解けないなりにもボチボチ近いところまでは行けてるっていうパターンが多いので(それ以上の問題は今のところ大体歯が立たない)、あと一歩何が足りなかったのかを探って書いていきます。

ARC074-D 3N Numbers

とりあえず、境界を決めて左側の数列から大きい順にN個、右側の小さい順にN個取ればよく、後はスコアが最大になる境界位置を全探索すればよい。ただし、これだとO(N^2logN)だからTLEになる。じゃあ、どうする?優先度付きキューか?っていうところまでは発想が至った。
ということで、左側数列用優先度付きキューと右側数列用優先度付きキューを用意して、境界位置を左から順にずらし、両優先度付きキューを更新していくってことを行おうとしたところ、右側数列用優先度付きキューから指定要素を取り除く操作が出てきて、その操作には1回O(N)かかるぞ、どうすればいいんだーとなる。
境界位置探索について1回で済まそうと、左側と右側の更新を同時(同じループ内)で行おうと考えたのが間違い。
まず境界位置を左からずらして、左側数列優先度付きキューを更新していく。それが終わった後、今度は境界位置を右からずらして、右側数列の優先度付きキューを更新していく。こうすれば優先度付きキューの操作は、先頭要素取り出しと要素追加だけになるので1回の操作はO(logN)で可能。最後にどの境界位置が最適かを見ればよい。

ARC083-D Restoring Road Network

与えられている最短距離表を各頂点間を結ぶ隣接行列だとみなしてワーシャル・フロイドで解いた後、矛盾がないか探る(矛盾があれば、ワーシャルフロイドで求めた最短距離と最短距離表の値が異なる部分が出てくる)。矛盾がなければ、頂点i,j間を直接結ぶ辺が必要かを考えればよい。つまり、i,j間を直接結ばずともどこかを経由して最短距離となる経路があるかを探す。
っていうところまでは出てきたけど、ここで躓く。i, j間の最短経路で直接結ぶ以外のパスがあるかを求める方法がわからなかった(i, jを直接結ぶパスのみを削除したグラフG'を用意して、これに対してダイクストラ法でi, j間の最短距離を求め、元のグラフで求めた最短距離と変わらないかってことをやろうとしたけどTLE)。
これを求めるには、頂点kを経由して、最短距離が等しくなる経路が存在するかを調べればよい。つまり、i-k間の最短距離+k-j間の最短距離が最短距離表にあるAijと等しくなるkが存在するかみればよい。

ARC076-D Built?

x座標、y座標それぞれでソートして、隣り合う街同士にのみ道を作ることを考えればOK。x座標、y座標ごった煮にして、コストの小さい道から順に追加していけばいいってところまでは考えた。ところが、やってみると各頂点をが全て連結されているのか効率よく判断する方法がわからず、言葉だけは知っていたUnion-Findを使うのかな?って思ったところで、なんとなく敷居が高い気がして断念。
というか、ここまで来れば、まんまクラスカル法。もし、クラスカル法知らなくても、ソートして隣り合う街同士のみ考えればいいってことを鑑みると、考慮する辺数が2Vだけなので、O(VlogV)で解ける最小全域木問題だって事気付けばプリム法で解いても良かった(これ解いた時、プリム法は知ってた)。いずれにしても、問題そのものが最小全域木の典型問題にもかかわらず、ソートとか隣り合う街同士のみ考慮とか色々考えていくうちに、最小全域木問題を解くアルゴリズム使うという発想がすっぽり頭から抜けていたのが反省点。

おわりに

もう少しレベルを上げたいので、時々こういう反省ポエムも書いていきます。