ハイパーデクノボウブログ

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

近接勾配法とproximal operator

はじめに

前回前々回とl1ノルム正則化項をつけた離散フーリエ変換により、スパースな解が得られること、そして基底ベクトルを並べた行列{\bf W}がユニタリ行列(直交行列)のため解析的に解けることを見てきました。それでは、{\bf W}が一般の行列の場合どうするか?今度こそ解析的に解けないので、反復法により最適解を求めます。しかし、目的関数にはl1ノルムという微分できない項があるため、最急降下法などのアルゴリズムは使えません。
このような微分できない項を含む最適化問題を効率的に解く手法の1つに近接勾配法(Proximal Gradient Method)があります。今回はこの近接勾配法と近接勾配法において重要な役割を果たすproximal operatorについて書いていきます。また、pythonで近接勾配法の実装を行いましたので、そのソースコードも載せておきます。

近接勾配法、proximal operatorとは

まずは、近接勾配法で解ける最適化問題の定義です。

\displaystyle\min_{{\bf x}} f({\bf x} )+ g({\bf x}) \tag{1}

ただし、f({\bf x})微分可能な凸関数、g({\bf x})微分可能とは限らない凸関数*1です。
l1ノルム正則化項付きの問題を式(1)に当てはめると、g({\bf x}) = \lambda \|{\bf x}\|_1となります。

次にproximal operatorの説明です。とりあえず定義を見てみましょう。関数g({\bf x})のproximal operatorは以下のように定義されます。

\displaystyle {\rm prox}_{\gamma g}({\bf v}) = \arg \min_{{\bf x}} g({\bf x}) + \frac{1}{2 \gamma} \|{\bf x} - {\bf v}  \|^2_2 \tag{2}

ここで\gammaは後述の更新式のステップ幅です。
このproximal operatorは何者か?というのは、g({\bf x})として以下のような指示関数を考えると見えてきます。

\begin{eqnarray}  
g({\bf x}) = I_c(\bf{x}) &=& \left\{  
\begin{array}{ll}
0 & ({\bf x} \in C) \\
\infty & ({\bf x} \notin C) \\ 
\end{array}
 \right.
\end{eqnarray}  \tag{3}

ここでCは閉凸集合です。この関数は{\bf x}が集合Cの要素なら0、Cの要素でないなら\inftyに吹っ飛ばすという関数です。そして、この指示関数のproximal operatorは以下のようになります。

\displaystyle \begin{eqnarray} {\rm prox}_{\gamma g}({\bf v})
&=& \arg \min_{{\bf x}}  I_c(\bf{x})  + \frac{1}{2 \gamma} \|{\bf x} - {\bf v}  \|^2_2 \\
&=&  \arg \min_{{\bf x} \in C} \|{\bf x} - {\bf v}  \|^2_2 \\
&=&   \Pi_C({\bf v})\\
\end{eqnarray}  \tag{4}

{\bf v}が集合Cの要素の場合、そのまま{\bf x} = {\bf v}とすれば\|{\bf x} - {\bf v}  \|^2_2を最小化できます。一方で、{\bf v}が集合Cの要素ではない場合、集合Cの要素で{\bf v}とのユークリッド距離を最小とするものが選ばれます。すなわち、{\bf v}の集合Cへの射影です(このような射影作用素\Pi_C({\bf v})として定義しています)。
よって、proximal operatorは射影作用素の関数gへの一般化と考えることができ、{\bf v}gによって定められる領域からはみ出していたら、その領域に入るように射影する作用素と見なせます。
proximal operatorの具体例を挙げると、g({\bf x})g({\bf x}) = \lambda\| \bf{x}\|_1というl1ノルムであれば、前回述べたソフトしきい値作用素がproximal operatorになります。つまり、

 \begin{eqnarray}  
{\rm prox}_{\gamma g}({\bf v}) = S_{\gamma \lambda} ({\bf v}) &=& \left\{
\begin{array}{ll}
 v_i - \gamma \lambda & (v_i \geq \gamma \lambda) \\
0 & (- \gamma \lambda  < v_i < \gamma \lambda) \\
v_i + \gamma \lambda & (v_i \leq  -\gamma  \lambda) \\
\end{array}
 \right.
\end{eqnarray}  \tag{5}

です。そして、近接勾配法の更新式はこのproximal operatorを使って以下のようになります。

\displaystyle{\bf x}_{t+1} = {\rm prox}_{\gamma g}({\bf x}_t - \gamma \nabla f({\bf x}_t)) \tag{6}

この更新式の収束性や正しさなどの説明は他の文献を参照いただくとして、直観的には次のような解釈ができると思います。
まず、{\bf x}_t - \gamma \nabla f({\bf x}_t)最急降下法の更新式です。また、上述のようにproximal operatorはgによって定められる領域に射影する作用素です。従って、最急降下法の1ステップで得られる値がgの領域からはみ出してたら、その領域へ射影するという操作を繰り返して最適解を探すアルゴリズムと見ることができます。
このように、関数f({\bf x})の勾配と関数g({\bf x})のproximal operatorがわかれば、近接勾配法により最適解が得られます。
なお、\gammaの決め方も収束のためには重要ですが、こちらについても別の文献を参照してください。

Moreau decompositionと双対ノルム

さて、ここまで見たように近接勾配法を使うためにはproximal operatorを求める必要があります。しかし、関数g({\bf x})によっては、式(2)から直接的にproximal operatorを求めるのが難しいことがあります。例えば、g({\bf x})=\lambda \| {\bf x} \|_2とl2ノルムとした場合(2乗が付いていないことに注意)、式(2)を解こうとして{\bf x}微分しても、1/\| {\bf x} \|_2がかかる項が出てきてしまうので厄介です。
そこで便利なのが以下に示すMoreau decompositionが成り立つというproximal operatorの性質です。

\displaystyle {\bf v} = {\rm prox}_{g}({\bf v}) +{\rm prox}_{g^*}({\bf v})  \tag{7}

ここで、{\rm prox}_{g^*}({\bf v})g({\bf x})の共役関数g^*({\bf x})のproximal operatorです。この共役関数とは下記の式で表される関数のことを言います。

\displaystyle g^*({\bf z}) = \sup_{{\bf x}} {\bf z}^T {\bf x} - g({\bf x})  \tag{8}

よって、この共役関数のproximal operatorさえ分かれば、元の関数のproximal operatorも求められます。
では、g({\bf x}) = \lambda\| \bf{x}\|というノルム関数の場合どうなるかを見ていきます。まずノルム関数の共役関数は以下の指示関数となります。

 \begin{eqnarray}  
g^*({\bf z}) = I_{\it B} \left({\bf z} \right) &=& \left\{  
\begin{array}{ll}
0 & (\|z \|_{*} \leq \lambda) \\
\infty & (\|z \|_{*} > \lambda) \\ 
\end{array}
 \right.
\end{eqnarray}\tag{9}

ここで、\| \bf{x}\|_{*}は双対ノルムと呼ばれるもので、以下の式で定義されます。

\displaystyle\| {\bf z} \|_{*} = \sup_{{\bf x}} {\bf z}^T {\bf x} \quad  ({\rm s.t.} \quad \| {\bf x} \| \leq 1)  \tag{10}

ということで、共役関数のproximal operatorは式(4)でみたように下記の射影作用素になります。

 {\rm prox}_{g^*}({\bf v}) = \Pi_{\it B} ({\bf v}) \tag{11}

ただし、{\it B}{\it B}=\left\{ {\bf x} \, | \,   \| \bf{x} \|_* \leq \lambda \right\}で表される集合です。
つまり、{\bf v}\| {\bf v} \|_* \leq \lambdaの範囲に射影する作用素がノルム関数の共役proximal operatorとなります。
おそらく式だけをザッと書いているので何故ノルムの共役関数が指示関数なのか?双対ノルムってそもそも何?どうやって求めるの?などの疑問が湧くかと思います。その辺りを少し補足しておきます。

双対ノルムについての補足

説明が前後しますが、後の説明のためにまずここで双対ノルムについて補足を書きます。そのための準備としてヘルダーの不等式という不等式を紹介します。
 1 \leq p,q \leq \inftyという2つの整数に、1/p + 1/q = 1という関係があれば、

{\bf z}^T{\bf x} \leq \| {\bf x}\|_p  \| {\bf z} \|_q \tag{12}

が成り立つというのがヘルダーの不等式です。実は、 \| \cdot \|_p の双対ノルムは \| \cdot \|_qであり、その逆も成り立ちます。そのため、l2ノルムの双対ノルムはそのままl2ノルムですし、l1ノルムの双対ノルムはl\inftyノルムです。
このことはヘルダーの不等式から示すことができます。まずヘルダーの不等式より、任意の\|{\bf x}\|_p \leq 1に対して\| {\bf z} \|_q \geq {\bf z}^T {\bf x}が成り立ちます。lpノルムの双対ノルム\|{\bf z} \|_*\|{\bf x}\|_p \leq 1という条件下での{\bf z}^T {\bf x}の最大値です。ここで{\bf x}をうまくとることにより、等式を成立させることができます。つまり、

\displaystyle \begin{eqnarray}
\| {\bf z} \|_{*} &=& \sup_{{\bf x}} {\bf z}^T {\bf x} \quad ({\rm s.t.} \quad \| {\bf x} \|_p \leq 1) \\
& =& \| {\bf z} \|_q \\
\end{eqnarray}
\tag{13}

が成立します。

何故、ノルム関数の共役関数が指示関数なのか?

まず、\| {\bf z} \|_* > \lambdaの時を考えます。式(8)の定義よりg^*({\bf z}) = \sup_{{\bf x}} {\bf z}^T {\bf x} - \lambda \| {\bf x} \| なわけですが、双対ノルムの定義式(10)と前述の前提より、 {\bf z}^T {\bf x} > \lambdaかつ\|{\bf x}\| \leq 1となる{\bf x}が存在することになります。これを満たす{\bf x}{\hat{\bf x}}とした場合、その{\hat{\bf x}}に定数tをかけたt {\hat{\bf x}}も式(8)の最適解の候補として入ってきます。これはもちろん、t \rightarrow \inftyとした場合も同様です。これらのことから、

 t{\bf z}^T {\hat {\bf x}} - \lambda \| t{\hat {\bf x}} \| =  t({\bf z}^T {\hat {\bf x}} - \lambda \| {\hat {\bf x}} \| ) \rightarrow \infty \tag{14}

となるため、\| {\bf z} \|_* > \lambdaのときはg^*({\bf z}) = \inftyとなります。
一方で、\| {\bf z} \|_* \leq \lambdaの時は、式(12)のヘルダーの不等式より{\bf z}^T{\bf x} \leq \| {\bf x}\| \| {\bf z} \|_* なので、

{\bf z}^T{\bf x}  - \lambda \| {\bf x}\| \leq \| {\bf x}\| \| {\bf z} \|_* -  \lambda  \| {\bf x} \| = \| {\bf x}\| ( \| {\bf z} \|_* - \lambda) \tag{15}

です。前提より、\| {\bf z} \|_* \leq  \lambdaなので、左辺(つまりはg^*({\bf z}) )が0以上になることはありません。ここで、{\bf x}= {\bf 0}という取り方もできるので、式(15)の左辺を0にすることができます。よって、\| {\bf z} \|_* \leq \lambdaのときはg^*({\bf z}) = 0となります。

結局のところproimal operatorはどうやって求めればよいか?

2つ方法があると思います。1つ目は式(2)から直接求める方法です。例えば、前回記したようの方法でl1ノルム正則化項のproximal operator式(5)を求めるというものです。
もう1つの方法としては、式(7)のMoreau decompositionの性質を使う方法です。特に正則化項がノルムなら、共役関数のproximal operatorは{\it B}=\left\{ {\bf x} \, | \,   \| \bf{x} \|_* \leq \lambda \right\}という超球への射影なので簡単に求められます。例えば、l2ノルムの双対ノルムはl2ノルムなので、

\displaystyle\begin{eqnarray}
{\rm prox}_{g^{*}} ({\bf v})= \Pi_{\it B}({\bf v}) = \left\{
\begin{array}{ll}
\lambda \frac{\bf{v}}{\| \bf {v} \|_2} \quad (\|  {\bf v} \|_2 > \lambda) \\
{\bf v}  \quad (\|  {\bf v} \|_2  \leq \lambda) \\
\end{array}
\right. 
\end{eqnarray}
\tag{15}

であり、求めたいproximal operatorは、

\displaystyle\begin{eqnarray}
{\rm prox}_{g} ({\bf v})= {\bf v} - \rm {prox}_{g^*}  ({\bf v})= \left\{
\begin{array}{ll}
{\bf v} -\lambda  \frac{{\bf v}}{\| {\bf v} \|_2} \quad (\|  {\bf v} \|_2 > \lambda) \\
0 \quad (\|  {\bf v} \|_2  \leq \lambda) \\
\end{array}
\right. 
\end{eqnarray}
\tag{16}

となります。

なお、式(2)にはステップ幅の\gammaが付いています。\gammaが付いていても基本的には上式の\lambda\gamma \lambdaに置き換えればよいだけですが、\gammaを含めたMoreau decompositionの一般形を書いておきましょう。

\displaystyle {\bf v} = {\rm prox}_{\gamma g}({\bf v}) +\gamma{\rm prox}_{g^*/\gamma} \left(\frac{{\bf v}}{\gamma} \right)  \tag{17}

という形式になります。したがって、gがノルム関数の場合は、

\displaystyle {\rm prox}_{\gamma g}( {\bf v}) = {\bf v} - \gamma{\rm prox}_{g^*/\gamma} \left(\frac{{\bf v}}{\gamma} \right)  =  {\bf v} - \gamma \Pi_{\it B} \left(\frac{{\bf v}}{\gamma} \right) \tag{18}

になります。ちなみに、多くの文献ではg({\bf x })=\| {\bf x}\|としているので、{\it B}=\left\{ {\bf x} \, | \,   \| \bf{x} \|_* \leq 1 \right\}として理論を展開している点に注意してください。

様々な正則化項のproximal operator

ここまでノルム関数が正則化項の時のproximal operator導出方法を見てきました。ここで、いくつかの正則化項に対するproximal operatorを記載しておきます。

l2ノルムの2乗正則化g({\bf v}) = \lambda \| {\bf v} \|^2_2

リッジ回帰等でよく使われる正則化項です。普通にproximal operatorの式を微分して0とおけば求められます。2乗しているので、以下のように式(16)とは異なる結果が得られます。

\displaystyle prox_{\gamma g} ({\bf v})= \frac{1}{1+2\gamma\lambda} {\bf v} \tag{19}

Elastic Net:g({\bf v}) = \lambda_1 \| {\bf v} \|_1 + \lambda_2 \| {\bf v} \|^2_2

l1ノルム正則化とl2ノルム2乗正則化の混合です。l1ノルムは2つ変数が同程度の寄与度だったとしてもどちらかしか選ばれない(一方が0になってしまう)ケースが多いですが、Elastic Netではこの点を改善できるようです。このproximal operatorはl1ノルム正則化の時と同様、変数ごとに分離して求めることができます。

\displaystyle prox_{\gamma g} ({\bf v})= \frac{1}{1+2\gamma\lambda_2} S_{\gamma \lambda_1} ({\bf v}) \tag{20}

グループl1ノルム(グループに重複がない場合):\displaystyle g({\bf v}) = \lambda \sum_{g \in G}  \| {\bf v}_g \|_p

変数がいくつかのグループに分けられるとき、それぞれのグループのlpノルムを取って、更にその和を取るノルムです( ベクトル{\bf v}のうち、グループgに所属する変数を並べたのが{\bf v}_g)。これを正則化項に用いるとグループごとにスパースになる、ならないが決まります。このグループl1ノルムを正則化項として使った回帰をgroup lassoと言います。
これはグループに重複がなく完全に分離ができるため、グループごとにlpノルムのproximal operatorを求めればよいです。p=2の時のグループgのproximal operatorは下記の通りです。(式(16)とほぼ同じです。)

\displaystyle\begin{eqnarray}
{\rm prox}_\gamma{g} ({\bf v}_g)= \left\{
\begin{array}{ll}
{\bf v}_g  -\gamma \lambda  \frac{{\bf v}_g}{\| {\bf  v}_g \|_2} \quad (\|  {\bf v}_g \|_2 > \gamma \lambda) \\
0 \quad (\|  {\bf v}_g \|_2  \leq \gamma \lambda) \\
\end{array}
\right. 
\end{eqnarray}
\tag{21}

ちなみに、グループに重複がある場合も提案されていますが、ここでは割愛します(まだ追いきれてない)。

トレースノルム :g({\bf Z}) = \lambda \| {\bf Z} \|_{{\rm trace}} = {\rm Tr}\left(\sqrt{{\bf Z}^T {\bf Z}}\right)

最後は行列関数に関する正則化項です。トレースノルムは\sqrt{{\bf Z}^T {\bf Z}}のトレースを正則化項とするものなので、{\bf Z}の特異値を\sigma({\bf Z})=\left[\sigma_1, \sigma_2, \cdots, \sigma_m \right]とすると、 \| {\bf Z} \|_{{\rm trace}} = \sum_{i=1}^{m} \sigma_iです。そのため、このトレースノルムを正則化項として使うと特異値が0になりやすく、結果低ランクな解が得られるようになります。
それではproximal operatorを見ていきます。{\bf Z}特異値分解{\bf Z} = {\bf U} {\rm diag}(\sigma ({\bf Z})) {\bf V}^Tとすると、proximal operatorは以下のようになります。

\displaystyle {\rm prox}_{\gamma g} ({\bf Z}) =  {\bf U} {\rm diag}(S_{\gamma \lambda}(\sigma ({\bf Z}))) {\bf V}^T \tag{22}

特異値分解して特異値に対してソフトしきい値作要素を施した後に再構築すればよいわけです。

近接勾配法を実装

ということで式(1)を解く近接勾配法を実装しました。

近接勾配法を行う関数

ということで、まずは近接勾配法そのものの実装です。

def proximal_gradient(grad_f, prox, gamma, objective, init_x, tol = 1e-9):
    x = init_x
    result = objective(x)
    while 1:
        x_new = prox(x - gamma * grad_f(x), gamma)
        result_new = objective(x_new)
        if (np.abs(result - result_new)/np.abs(result) < tol) == True :
            break;
        x = x_new
        result = result_new
    return x_new, result_new

引数について説明します。
grad_fは\nabla f({\bf x})にあたります。つまり、grad_fはxを引数として、その勾配を返す関数です。
proxはproximal operatorです。proxも関数で、入力変数vとステップ幅のgammaの2引数をとり、proximal operatorの計算結果を返します。
objectiveは目的関数です。objectiveはxを引数として、目的関数の計算結果を返す関数です。
init_xは初期値、tolは収束の許容誤差です。
やっていることは、式(6)の更新式をそのまま素直に実装しているだけです。

l1ノルム正則化項のproximal operator

次にproximal operatorの実装です。

def prox_norm1(v, gamma, lam):
    return soft_thresh(v, lam * gamma)

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

prox_norm1がl1ノルムのproximal operatorを計算する関数です。補助としてsoft_threshという関数を作っています。

lasso回帰問題を解く関数

以下は、\min_{{\bf x}}\| {\bf y} - {\bf W}{\bf x} \|_2^2 + \lambda \| {\bf x}\|_1というlasso回帰問題を解く関数です。なお、l1ノルム正則化項に対する近接勾配法はISTA(Iterative Shrinkage Thresholding Algorithm)と呼ばれるため、関数名をISTAとしています。

def ISTA(W, y, lam):
    objective = lambda x:np.sum(np.power(y - np.dot(W, x), 2))/2.0 + lam * np.sum(np.abs(x)) #(1)
    grad_f = lambda x: np.dot(W.T, np.dot(W, x) - y)  #(2)
    (u, l, v) = np.linalg.svd(W)   #(3)
    gamma = 1/max(l.real*l.real)   #(4)
    prox = lambda v, gamma:prox_norm1(v, gamma, lam)  #(5)
    x_init = randn(W.shape[1])  #(6)
    (x_hat, result) = proximal_gradient(grad_f, prox, gamma, objective, x_init, 1e-5)  #(7)
    return x_hat, result

(1)ではラムダ式を使ってobjectiveをxを引数としてとり、上記の目的関数を計算する関数として定義しています。
また、\nabla f({\bf x})={\bf W}^T ({\bf W} {\bf x} - {\bf y})であり、(2)はこの計算を行うgrad_fをxに引数とる関数として定義しています。
(3)、(4)は\gammaを定めています。最大特異値の逆数より小さい値を\gammaとして定めると収束するようですが、その辺りはまだ理解ができておりません。とりあえず、文献の通りに実装します。
(5)がproximal operatorです。prox_norm1関数はv,gamma,lamの3引数を受け取りますが、proximal_gradient関数に渡すproxはvとgammaの2引数のみ取る関数なので、ラムダ式を使って再定義しています。
(6)はランダムに初期値を決めて、(7)でproximal_gradient関数にて問題を解きます。

proximal operatorを別なものにする

上記ISTA関数のproximal operatorとobjectiveの部分だけ変えれば他の正則化にも対応可能です。例としてグループl1ノルム正則化(group lasso)の場合を示します。
まずは、グループl1ノルムのproximal operatorの実装です。

def prox_group_norm(v, groups, gamma, lam):
    u = np.zeros(v.shape[0])
    for i in range(1, max(groups) + 1):
        u[groups == i] = block_soft_thresh(v[groups == i] , gamma * lam)
    return u

def block_soft_thresh(b, lam):
    return max(0, 1 - lam/np.linalg.norm(b)) * b

groupsは1~G(Gはグループ数)を要素に持つリストで、各変数がどのグループに属しているかを表します。後はグループごとにblock_soft_thresh関数にて、l2ノルムのproximal operatorを計算しています。
そして、上記のISTA関数内の(5)でproxとして、prox_group_normを指定すればgroup lassoになります。

    prox = lambda v, gamma:prox_group_norm(v, groups, gamma, lam) #(5)'

もし、Elastic Netを使いたかったらここにElastic Netのproximal operatorを持ってくればよいし、他のノルムの場合も同様です(本当はobjectiveも変更する必要がありますが上記では省略しています)。
また、f({\bf x})を別なものにしたかったら\nabla fを別途定義して、grad_fとすればよいです。
行列関数の場合も同様で、grad_fとproximal operator (あと、objective)さえ定義してあげたら機能します。そのため、proximal operatorとしてトレースノルムを指定してあげると、低ランク解を得ることができます。

まとめ

今回は近接勾配法とproximal operatorについて見てきました。上記の実装については次の機会に応用事例で実験したいと思います(ちなみに、疑似データで実験したところそれっぽく動いていました)。
また、今回は取り上げませんでしたが、この界隈ではISTAの高速版FISTAや拡張ラグランジュ法、その改良であるAlternating Direction Method of Multipliers (ADMM) など様々なアルゴリズムが提案されており、非常に奥が深いです。

参考文献

最後にこの記事を書くにあたって参考にした文献です(基本的に前回とほぼ同じです)。

また日本語の書籍で、proximal operator界隈が詳しく書かれているのは以下2冊かなと思います(他にもあったら教えてください)。

*1:もうひとつ、「単純」という条件も付きます。

l1ノルム正則化フーリエ変換を複素数のまま解く

はじめに

前回の記事では、離散フーリエ変換を係数をLasso回帰として求めスパースな解を得る方法を書きましたが、この際に複素数を取り扱いたくないため実数のみの問題に変換して、ソルバに解かせるということを行いました。今回は、複素数のまま問題を解くことを考え、複素数の最適解を得る方針で考えます。この方針で考察することで、以下のことが見えてきたので記事に残そうと思った次第です。

  • 前回の記事では「この問題は解析的には解けず」などと書いているが、実は{\bf W}が離散フーリエの基底ベクトルを並べた行列であれば解析的に解けること
  • 結局は、離散フーリエ変換してしきい値未満の要素を0にしていることと(ほぼ)等価であること
  • 複素数としてのl1ノルムを考えると、別な解が得られること

ちなみに、自分は「この方法だとなくうまくいきそうだな」と思ったことはまず試してみて、後から理屈を考える(特にうまくいかなかったとき)というスタンスなので、理論的、数学的な厳密性はあまり考慮していません。したがって、本来理論的にはあり得ない式展開、考え方をしている場合もあるかもしれませんがご了承ください。

定式化

前回と同様、正則化項をつけた回帰問題として定式化します。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x}) \tag{1}

{\bf y}は元の時間軸の信号で、{\bf y} = \left[ y(0), y(1), \cdots, y(T - 1) \right]^T{\bf x}が求めたい離散フーリエ係数であり、{\bf x} =  \left[ x(0), x(1), \cdots, x(F - 1) \right]^Tです。ここで{\bf x}の各要素は複素数です。{\bf y}は実信号なので実数ベクトルですが、虚部がすべて0の複素数と見なしてもよいでしょう。また、行列{\bf W}の各列ベクトル{\bf w}_k ( k = 0, 1, \cdots, T-1)
{\bf w}_k = \frac{1}{T}\left[\exp(i\frac{2\pi k \times (0)}{T}) ,\exp(i\frac{2\pi k \times (1)}{T}), \cdots, \exp(i\frac{2\pi k \times (T - 1)}{T} ) \right]^T
とします。前回の定義とほぼ一緒ですが、各要素を要素数Tで割っています。こうすると、{\bf W}はユニタリ行列になります。つまり、{\bf W}^* {\bf W} = {\bar {\bf W}}^T {\bf W}  = {\bf I}です({\bf W}^*は随伴行列、 {\bar {\bf W}}複素共役を表します)。
後は正則化\lambda \Omega({\bf x})です。普通に考えると
\Omega({\bf x}) = \| {\bf x} \|_1 = | x(0)| + |x(1)|+ \cdots + |x(T-1)| \tag{2}

(ただし、|x|xの絶対値を表し |x|  = \sqrt{{\bar x} x})としたいところですが、実はこの定義だと前回と同じ結果にはなりません。このことは後ほど説明するとして、今は
\Omega({\bf x}) = | {\rm  Re}( x(0))| +  | {\rm  Im}( x(0))| + \cdots +  | {\rm  Re}( x(T - 1))| +  | {\rm  Im}( x(T - 1))|  \tag{3}

とします({\rm Re(x)}xの実部、{\rm Im(x)}xの虚部を表す)。式(2)と式(3)の大きな違いは、実部と虚部が関連し合っているか、完全に独立なものとして扱うか、です。そのため、式(3)を正則化項として用いることは、実部と虚部を分けて解を求める前回の考え方と一致します。

解析解の導出と実装

当初は「スパースモデリングのための凸最適化」の解説論文を読んで近接勾配法で解こうなんて考えていましたが、解説を読んでいくと実は、{\bf W}がユニタリ行列であるため解析的に解が得られることがわかりました。
ということでその導出です。まず、解説に従うと、式(1)は以下のように式を変形できます。
\begin{eqnarray} 
\displaystyle (1) &=& \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x})  \\
\displaystyle  &=& \min_{\bf x} \| {\bf W}^*{\bf y} -  {\bf W}^*{\bf W}{\bf x}\|^2_2 + \lambda \Omega({\bf x}) \\
\displaystyle  &=&  \min_{\bf x} \| {\bf b} -  {\bf x}\|^2_2 + \lambda \Omega({\bf x}) \ \ \  (\because {\bf b} =  {\bf W}^*{\bf y})\\
\displaystyle  &=&  \min_{\bf x}  \sum_{i = 0}^{T - 1} \left\{ |(b(i) -  x(i))|^2+ \lambda( | {\rm  Re}( x(i))| +  | {\rm  Im}( x(i))| )  \right\} \\
\displaystyle  &=&  \min_{\bf x}  \sum_{i = 0}^{T - 1} \left[  \left\{ {\rm Re} (b(i)) -  {\rm Re} ( x(i) ) \right\} ^2 + \lambda | {\rm  Re}( x(i))|   +  \left\{ {\rm Im} (b(i)) -  {\rm Im} ( x(i) ) \right\} ^2 + \lambda | {\rm  Im}( x(i))|   \right]   \\ 
\displaystyle  &=&\sum_{i = 0}^{T - 1}   \min_{x(i)}   \left[  \left\{ {\rm Re} (b(i)) -  {\rm Re} ( x(i) ) \right\} ^2 + \lambda | {\rm  Re}( x(i))|   +  \left\{ {\rm Im} (b(i)) -  {\rm Im} ( x(i) ) \right\} ^2 + \lambda | {\rm  Im}( x(i))|   \right]    
\end{eqnarray}
このように式(1)の問題は各要素独立に最適解を選ぶ問題に変形されます。しかも、実部と虚部も分離しているので、こちらについても独立して最適解を選ぶことができます。ようするに、以下の最適化問題を各要素の実部、虚部それぞれ解けばよいわけです。
\displaystyle\min_{x}    (b-  x)^2 + \lambda | x|  \tag{4}
この解は次のような式で得られます。
\begin{eqnarray} 
{\hat x} = S_{\lambda} (b) &=& \left\{
\begin{array}{ll}
 b - \frac{\lambda}{2} & (b \geq \frac{\lambda}{2}) \\
 0 & (- \frac{\lambda}{2}  < b < \frac{\lambda}{2}) \\
b + \frac{\lambda}{2} & (b \leq  - \frac{\lambda}{2}) \\
\end{array}
 \right.
\end{eqnarray} \tag{5}
解説論文と異なり式(1)および式(4)の第1項に1/2がかかっていないため、最適解も\lambda1/2かけていることに注意してください。ここで式(5)のS_{\lambda}(b)はソフトしきい値作用素と呼ばれる作用素です。
このソフトしきい値作用素を用いると、今回得たい式(1)の最適解は下記のように表せます。
{\hat{\bf x}} = S_{\lambda}({\rm Re}({\bf W}^* {\bf y}))+ i S_{\lambda}({\rm Im}({\bf W}^* {\bf y})) \tag{6}

さて、ここで通常の離散フーリエ変換の係数は{\bf W}^* {\bf y}で得られることを思い出すと、式(6)は通常の離散フーリエ変換の結果の実部、虚部それぞれに対してソフトしきい値作用素を作用させていることになります。
これは式(5)により離散フーリエ係数の絶対値が\lambda/2未満であれば0にすることに対応します。l1ノルム正則化項の離散フーリエ変換なんて言っていますが、紐解いていくと離散フーリエの結果に対ししきい値未満ならば0にするという割りとよくやる操作とほぼ同じになるわけですね(しきい値を超えている要素には\pm \lambda/2を加える点は異なりますが)。
ということで、以下は上式によるスパース離散フーリエ変換を行う関数の実装です。

def sparse_dft_with_thresh(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (soft_thresh((np.dot(W.T.conj(), y)).real, lam/2)
            +1j * soft_thresh((np.dot(W.T.conj(), y)).imag, lam/2))

def soft_thresh(b, lam):
    x_hat = np.zeros(b.shape[0])
    x_hat[b >= lam] = b[b >= lam] - lam
    x_hat[b <= -lam] = b[b <= -lam] + lam
    return x_hat

sparse_dft_with_threshがスパースな離散フーリエ変換を行う関数で、soft_threshがソフトしきい値作用素を施す関数です。soft_threshを解説と合わせるために、引数で渡す時にlam/2をしています。

複素数l1ノルム正則化

ここで1つ疑問が残ります。正則化として式(2)のノルムを使うとどうなるか?です。これまでの議論から、式(2)のノルムを使った場合、式(1)は次のように変形できます。
\displaystyle \sum_{i = 0}^{T - 1} \min_{x(i)}  \left\{ |(b(i) -  x(i))|^2+ \lambda| ( x(i)| \right\} \tag{7}
この式を見ると要素同士は独立して最適解を選ぶことができますが、実部と虚部は独立させることができません(|x(i)|  = \sqrt{{\bar x(i)} x(i)}なので)。
この最適解を得るために、まず以下のような式を考えます。
\displaystyle {\rm prox}_{\lambda}({\bf b})=  \min_{\bf x}  \frac{1}{2} \| {\bf b} -  {\bf x} \|_2^2+ \lambda \Omega({\bf x}) \tag{8}
この式(8)は近接作用素(proximal operator)と呼ばれるもので、近接勾配法などで重要な役割を持つ作用素です。今、式(7)について複素数を実部と虚部を並べた2次元ベクトルとして考えると(つまり、{\bf x} = \left[{\rm Re}(x) , {\rm Im}(x)\right]^T )、式(7)と式(8)の問題は(第1項の1/2を除いて)等価になります(\Omega({\bf x}) = \| x \|_2)。
ということで、式(8)を解くことを考えますが、実はこれはgroup lassoにおけるproximal operatorを求めることと等価です。group lassoのproximal operatorの導出やその基礎となる劣微分、凸解析、双対理論などは「ptimization with sparsity-inducing penalties」「Proximal Algorithms」にまとまっています。正直、まだ理解できていない部分が多いですが、これらの資料に記載されているgroup lassoにおけるproximal operatorを使うと式(7)の解は下記のようになります(たぶん)。
\begin{eqnarray} 
{\hat x} = S_{\lambda} (b) &=& \left\{
\begin{array}{ll}
 b -   \frac{\lambda b}{2|b|}& (|b| > \frac{\lambda}{2}) \\
0 & (|b| \leq  - \frac{\lambda}{2}) \\
\end{array}
 \right.
\end{eqnarray} \tag{9}

式(6)では、実部、虚部それぞれでしきい値未満だと0になりましたが、こちらは複素数の絶対値がしきい値未満だと0になります。つまり、実部と虚部は運命共同体で、しきい値未満となる要素は実部も虚部一緒に0になり、その逆もまたしかりです。ということで実装です。

def group_lasso_prox(b, lam):
    x_hat = np.zeros(b.shape[0], "complex")
    x_hat[np.abs(b) >= lam] = (b[np.abs(b) >= lam] -
                                lam *  b[np.abs(b) >= lam]/np.abs(b)[np.abs(b) >= lam])
    return x_hat

def sparse_dft_with_complex_l1_norm(y, T, lam):
    Z = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (Z.T * u) * (2 * np.pi /T)
    W= np.exp( 1j * V)/T
    return (group_lasso_prox(np.dot(W.T.conj(), y), lam/2))

一応、実験した結果も載せておきます。実験には前回同様sin波とcos波とホワイトノイズを重ね合わせた疑似データを用います。

def main():
    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T))

    x1 = sparse_dft_with_thresh(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x1.real, "r")
    plt.hold(True)
    plt.plot(x1.imag, "b")

    x2 = sparse_dft_with_complex_l1_norm(y, T, 0.5)
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-2, 2))
    plt.plot(x2.real, "r")
    plt.hold(True)
    plt.plot(x2.imag, "b")

    plt.show()


実部と虚部を独立の変数として扱った正則化の結果
f:id:YamagenSakam:20180204165656p:plain


複素数l1ノルムによる正則化の結果
f:id:YamagenSakam:20180204165749p:plain


こうやって見ると、前者の方がスパース性という意味だと高いです。実部、虚部を運命共同体にする理由がない限り、前者のように独立して扱った方がよいように思えます。

最後に

ここまで見て感じたかと思いますが、結局のところ単なるしきい値での枝刈りに帰着するので、実装まで落としこんでしまうと正直あまり面白いものではありませんし、応用としてもあまり意味のあるものでもありません(これやるくらいなら、普通にFFTしてしきい値未満を0にする方が断然早いです)。
ただ、ここまで導出する過程で、ソフトしきい値作用素やproximal operator、劣微分、凸解析などなど最適化理論の色々なトピックスに触れることができました(まだ全然理解できていないですが)。
冒頭に書いた思いついたことはまず試してみるというスタンスは、試すために知らないトピックスを調べるきっかけになるというメリットがあります。やってみてダメなら何がダメか更に深追いするモチベーションにもなります。後はこういう試行錯誤は何より楽しいですね。ということで、今後は近接勾配法や凸解析の理論などの知識を深めていければと思います。

l1ノルム正則化でスパースな離散フーリエ変換

はじめに

離散フーリエ変換は離散化された時間軸の信号を周波数軸に写像する基底変換の一種で、周波数解析などによく使われます。この離散フーリエ変換の係数は、基底ベクトルの重み付き線形和で表したときの誤差の2乗和を最小にする重みと見なせます。つまり、離散フーリエ係数は一種の線形回帰問題の解として得られます。ということは、l1ノルム正則化項つけたLasso回帰として解けばスパースな離散フーリエ変換ができるのでは?ということでやってみました。

離散フーリエ変換の線形回帰的定式化

元の信号のデータ点数をT個、基底ベクトルの個数をF個とし、上述の通り回帰問題として離散フーリエ係数を求める問題を定式化すると下記のようになります。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 \tag{1}

ここで、{\bf y}は元の時間軸の信号で、{\bf y} = \left[ y(0), y(1), \cdots, y(T - 1) \right]^Tです。
{\bf W} = \left[{\bf w}_0, {\bf w}_1, \cdots, {\bf w}_{F - 1} \right]F個の基底ベクトルを劣ベクトルとして並べたT \times Fの行列です。
一般的に離散フーリエ変換では{\bf W}{\bf w}_k = \left[\exp(i\frac{2\pi k \times (0)}{T}) ,\exp(i\frac{2\pi k \times (1)}{T}), \cdots, \exp(i\frac{2\pi k \times (T - 1)}{T} ) \right]^Tという複素数の要素を持つ行列で、各列ベクトル同士は直交します({{\bf w}_i^{*} {\bf w}_j} = 0 \quad i \neq j)。
そして、\bf{x}が求めたい離散フーリエ係数であり、{\bf x} =  \left[ x(0), x(1), \cdots, x(F - 1) \right]^Tです。
通常、F = Tです。そのため、(1)の最適解は{\hat{\bf x}} = {\bf W}^{-1} {\bf y} =  {\bf W}^{*} {\bf y}で表され、最小2乗誤差は0になります。

l1ノルムによるスパース化

ということで、回帰問題としてみなせるならl1ノルムの正則化項をつける(つまりは、Lasso回帰を行う)とスパースな解が得られるだろうという考えです。(1)の目的関数に正則化項を付け加えます。

\displaystyle \min_{\bf x} \| {\bf y} - {\bf W}{\bf x}\|^2_2 + \lambda \|  {\bf x} \|_1 \tag{2}

こうなると最早この問題は解析的には解けず、正則化項は\|{\bf x} \|_1 \leq \gammaという制約条件を付けるのと等価なので、最適解による最小2乗誤差も一般的には0にはなりません。ということで、cvxpyというソルバを用いて上記最適化問題を解きます。
なお、この定式化は入口こそ離散フーリエ変換を回帰問題と見なしてLasso回帰を行うという考えですが、結局のところこれは、この資料等に記載されいるノイズあり圧縮センシングの一種です。 

複素数最適化問題から実数のみの最適化問題

離散フーリエ変換は上述のように、基底ベクトルが複素数で得られる解も一般的に複素数となります。ところが、このままだと実装上面倒くさい(というより、ソルバが動いてくれるかわからない)ので実数のみの式に変換します。といっても、\exp({i \theta}) = \cos\theta + i\sin\thetaより実部と虚部に分けてそれぞれ個別の変数として考えるだけです。
具体的には、{\bf W} {\bf W} = \left[ {\bf w}_0^c, {\bf w}_1^c, \cdots, {\bf w}_{T - 1}^c,  {\bf w}_0^s, {\bf w}_1^s \cdots, {\bf w}_{T - 1}^s  \right]というT \times 2Tの行列として再定義します。ここで、
{\bf w}_k^c = \left[\cos(\frac{2\pi k \times (0)}{T}) ,\cos(\frac{2\pi k \times (1)}{T}), \cdots, \cos(\frac{2\pi k \times (T - 1)}{T} )\right]^T
{\bf w}_k^s = \left[\sin(\frac{2\pi k \times (0)}{T})  ,\sin(\frac{2\pi k \times (1)}{T}), \cdots, \sin(\frac{2\pi k \times (T - 1)}{T} )\right]^T
です。あわせて説明変数{\bf x}を要素数2Tのベクトルとして式(2)を解けばよいです。この変換で基底ベクトルの数も説明変数の要素数も倍になったように見えますが、元々も複素数でそれぞれが実部と虚部を持っていたため本質的には変わりません。
また、離散フーリエ変換の場合、x(k) + {\bar x(k + \frac{T}{2}})という関係があるため、実質0\sim\frac{T}{2}-1番目の要素までを求めれば、残りの\frac{T}{2}\sim T-1番目は自動で決定されます。これは、l1ノルム正則化項をつけても変わりません(実際に試しました)。そのため、説明変数の要素数2T\rightarrow Tへとシュリンクすることも可能ですが、今回は理解を深めるためにもそのまま解くようにしました。

cvxpyを用いた実装

上述の通り、今回は最適化問題を解くためにcvxpyというツールを利用しました。このcvxpyの使い方はこちらの記事を参考にしました。cvxpyの使い方は他の解説記事に譲るとして、早速式(2)を解く関数のソースコードを示します。

from cvxpy import *
import numpy as np
from numpy.random import *
import matplotlib.pyplot as plt

#y:原信号、T:信号のデータ点数、lam:ハイパーパラメータλ
def sparse_dft(y, T, lam):
    #行列Wを生成
    V = np.meshgrid(np.arange(T),np.zeros((T,1)))[0]
    u = np.arange(T)
    V = (V.T * u) * (2 * np.pi /T)
    W = np.c_[np.cos(V), np.sin(V)]

    #変数xを定義
    x = Variable(T * 2)

    #ハイパーパラメータを定義
    lamda = Parameter(sign="positive")
    lamda.value =lam

    #目的関数を定義
    objective = Minimize(sum_squares((y - W * x))/T + lamda*norm(x, 1))
    p = Problem(objective)

    #最適化計算
    result = p.solve()

    return result, x

特に特別なことはしておらず、最初に{\bf W}を生成して、後はcvxpyの解説に従って変数、ハイパーパラメータ、目的関数を定義し、p.solve()で問題を解きます。

実験

以下のようなコードで疑似データを発生させて実験してみました。

    T = 2**10
    y = (3 * np.cos(28* 2*np.pi * np.arange(T)/T) +2 * np.sin(28* 2*np.pi * np.arange(T)/T) +
            np.cos(400* 2*np.pi * np.arange(T)/T) + np.sin(125 * 2*np.pi * np.arange(T)/T) + 5 * randn(T) )

いくつかの周波数が異なるsin波とcos波を重み付きで重ね合わせて、更に正規分布に従うノイズを乗せた1024点のデータです。疑似データの波形は以下の通り。
f:id:YamagenSakam:20180129000732p:plain

通常の離散フーリエ変換

まずは上記の疑似データにFFTにより離散フーリエ変換を行いました。

    Y = np.fft.fft(y)/T #スケーリングを合わすため要素数で割る
    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(Y.real , "r")
    plt.hold(True)
    plt.plot(Y.imag , "b")

結果は下記の図です。
f:id:YamagenSakam:20180129000744p:plain
重ねたsin波、cos波に対応する周波数成分が大きな値になるのは当たり前ですが、ノイズが重畳しているせいで全周波数成分まんべんなく値が入ります。

Lassoによる離散フーリエ変換

次は今回記した手法での離散フーリエ変換\lambda = 0.5で実験)。上述のsparse_dft関数で計算します。

    (result, x) = sparse_dft(y, T, 0.5)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-1.5, 1.5))
    plt.plot(x.value[0:T,0], "r")
    plt.hold(True)
    plt.plot(x.value[T:T*2,0], "b")

f:id:YamagenSakam:20180129000758p:plain
通常の離散フーリエ変換とは異なりスパースな解が得られ、余分な周波数成分は0になっていることがわかります。

逆変換

得られたフーリエ係数から逆変換を行い元の信号を再構築してみました。こうすることで、余計なノイズ成分が除去された信号が得られるはずです。
なお、逆変換のため、得られた{\bf x}から\cosにかかる係数を実部、\sinにかかる係数を虚部として複素数に変換してから、ifftメソッドで逆変換を行っています。

    X=(x.value[0:T] + (x.value[T:T*2] * 1j))*T
    y_hat = np.fft.ifft(X.T)

    plt.figure()
    plt.xlim((0, T))
    plt.xticks(np.arange(0, T , 100))
    plt.ylim((-20, 20))
    plt.plot(y_hat.real.T, "r")

以下が原波形から今回の方法でノイズを除去したデータです。
f:id:YamagenSakam:20180129000809p:plain
また、以下が上述の疑似データ生成のコードからノイズ項を取り除いたデータの波形です(つまり、真に得たい波形)。
f:id:YamagenSakam:20180129002851p:plain
確かにノイズが除去された波形が得られています。実際のノイズ項を除いた信号と比べると若干スケーリングが小さくなっていますが、悪くない結果だと思います。

まとめ

離散フーリエ変換を回帰問題と見なし、目的関数にl1ノルムの正則化項をつけて誤差の最小化を行うことで、スパースな離散フーリエ係数が得られました。またそれを逆変換することで原信号からノイズを除去した信号を得ることもできました。また、最適化問題を解くためにcvxpyというツールを利用しました。
今回は何も考えずcvxpyに問題を突っ込んで解いただけでしたが、今後は、Lasso回帰の最適化等で使われる近接勾配法やcvxpyの詳しい使い方、凸解析など最適化の理論、実装両面で理解を深めていきたいと思います。

はてなダイアリーからはてなブログへ移行

このたび、本ブログをはてなダイアリーからはてなブログへ移行しました。旧はてなダイアリーの記事にアクセスすると、対応する本ブログの記事へのリダイレクトされるようなので、ブックマーク登録等の変更は不要です。
さて、この移行作業に当たって、元のダイアリーで書いていた数式がうまく表示されず修正に手間取った箇所がいくつかあったので、備忘録の意味でもここに記しておきます。(大体は元の書き方がよろしくなかっただけですが)

小中大括弧でエラー

  • 元々の書き方
[tex:\(x + y \)]
[tex:\{x + y \}]
[tex:\left{x + y \right}]
[tex:\[x + y \]]
[tex:\left[x + y \right]]

見え方
\(x + y \) :おかしい。
\{x + y \}:うまく表示される。
\left{x + y \right}:おかしい。
[x + y ]:うまく表示される。
\left[x + y \right]:おかしい。

小括弧は\エスケープだけではエラーになり、中大括弧は\left, \rightをつけるとエラーになります。そのため、小括弧には\left、\rightをつける、中大括弧は\left\{ \right}\}のように\left、\rightの後、更に\でエスケープする方法で対策しました。

  • 修正後の書き方:
[tex:\left(x + y \right)]
[tex:\{x + y \}]
[tex:\left\{x + y \right\}]
[tex:\[x + y \]]
[tex:\left\[x + y \right\]]

見え方
\left(x + y \right):\left, \right追加でOKになった。
\{x + y \}:変更なし。
\left\{x + y \right\}:\エスケープ追加でOKになった。
[x + y ]:変更なし。
\left[x + y \right]:\エスケープ追加でOKになった。

sum、min、max等の上置き,下置きが上付き、下付きになる

  • 元々の書き方
[tex:\sum_{i = 1}^{N}x_i]
[tex:\min_{x} f(x)]

見え方
\sum_{i = 1}^{N}x_i
\min_{x} f(x)

本当は\Sigma\minの真上、真下に添え字が来るようにしたいのに、右上、右下に来ています。これは\displaystyleにすることで解決しました。

  • 修正後の書き方
[tex:\displaystyle \sum_{i = 1}^{N}x_i]
[tex:\displaystyle \min_{x} f(x)]

見え方
\displaystyle \sum_{i = 1}^{N}x_i
\displaystyle \min_{x} f(x)

\bfにより意図しないところまでも太字化される

  • 元々の書き方
[tex:(\bf{x}_N)]

見え方
(\bf{x}_N):xの先にある添え字のNや)まで太字化される。

わかりにくいですが、xだけ太字にしたいのに、それ以降の文字まで太字になっています。これは元々の書き方がそもそも間違っており、中括弧の位置を\bf{x}→{\bf x}と変更すると修正できます。

  • 修正後の書き方
[tex:({\bf x}_N)]

見え方
({\bf x}_N):xのみが太字。


ということで、今回は以上。せっかくはてなブログに移行しましたが、次回更新は遠い未来になるかも知れません。

NIPS2017論文メモ

お久しぶりです。2014年9月以来、実に3年ぶりの更新です。しばらく機械学習関連の勉強は滞っていたのですが、思うところがあって再開。とりあえず、昨年12月に開催されたNIPS2017のプロシーディングス中から面白そうな記事をピックアップして斜め読みしました。備忘録&リハビリのためにも、記事を更新してみよう思った次第です。

Deep Subspace Clustering Networks

部分空間クラスタリング+オートエンコーダで非線形クラスタリングを行う手法を提案。

部分空間クラスタリングは、同じクラスなら同じ(低次元)部分空間埋め込まれている(つまり、同じクラスに所属している別のサンプルの線形結合で表せる)という仮定の下クラスタリングを行う枠組みで、SSCLRRが有名。サンプル集合を\bf{X}として、いずれの手法も\bf{X}=\bf{X}\bf{C}という制約条件のもと、線形結合の行列\bf{C}のスパース化や低ランク化を行う(前者はl1ノルムの最小化、後者はトレースノルムの最小化)。\bf{C}を求めたら|\bf{C}|+|\bf{C}^T|をaffinity matrixとしてスペクトラルクラスタリングを行うのが一般的(他にも色々やり方があるらしい)。

これらの方法では非線形クラスタリングは行えない。そこで、オートエンコーダを利用しようというのがこの論文の提案。オートエンコーダはニューラルネットを入力と出力を同じデータとして学習させることで、中間層の出力\bf{z}非線形な特徴抽出となるという手法。入力から特徴を抽出する中間層まで層をエンコーダ層、特徴を抽出する中間層から出力までの層をデコーダ層と呼ぶ。この論文では、このエンコーダ層とデコーダ層の間に\bf{Z}\bf{C}とするSelf-Expressive層を新たに設けることを提案(\bf{Z}はエンコーダ層で特徴抽出されたサンプルの集合)。これは、線形の全結合となる。この\bf{C}と他の層の重みを求めるため、下記の損失関数を最小化。
L(\Theta, {\bf C} ) = \frac{1}{2} \|{\bf X} - {\bf \hat{X}_{\Theta}} \|^2_{F}+\lambda_1 \|{\bf C} \|_p + \frac{\lambda_2}{2} \|{\bf Z_{\Theta_e}} - {\bf Z_{\Theta_e}} {\bf C} \|^2_{F}  \quad  {\rm s.t.}  \quad  {\rm diag}( {\bf C} )={\bf 0}
\Thetaニューラルネットの重み係数集合、{\bf Z_{\Theta_e}} はエンコーダ層の出力集合(バッチ)を表す。ノルム項はl1l2の両方を実験。

上述のように、Self-Expressive層は線形の全結合層とみなせるので、バックプロパゲーションにより一括で学習できる。この方法で\bf{C}を求めたら後は従来通りにスペクトラルクラスタリング。結局のところ、オートエンコーダによって特徴抽出したデータに対して部分空間クラスタリングを行うが、その部分空間クラスタリング用の\bf{C}ニューラルネット上に組み込んで、一緒に学習してしまおうという手法。

顔画像クラスタリングのベンチーマークデータセットExtended Yale BORL)や物体画像クラスタリングベンチマークデータセットCOIL20COIL100)で実験。結果、他の手法(カーネル部分空間クラスタリングなど非線形クラスタリング手法含む)よりも高精度なクラスタリングが可能。ノルム項ついては、何故かl2ノルムを使った方が良い結果に。

Mixture-Rank Matrix Approximation for Collaborative Filtering

協調フィルタリングでよく使われるLow-rank matrix approximationにおいて、低ランクと高ランクの近似行列を混合することにより精度を上げた手法。
実際、レイティング数が少ないアイテムやユーザに対する近似行列は低ランクにした方が推定結果は良く、逆にレイティング数に対する近似行列は高ランクにした方が推定結果は良い。

従来の手法はユーザ・アイテムの実際のレイティングマトリックス\bf{R}として、\bf{R}\approx \bf{\hat{R}} = \bf{U}\bf{V}^T\bf{U}, \bf{V}ともにk列の行列)となる\bf{U}, \bf{V}を求め、i番目のユーザのj番目のアイテムに対する嗜好スコアを\bf{U}_i\bf{V}^T_j(それぞれ、\bf{U}のi行目、\bf{V}のj行目のベクトル)として推定する。PMFという手法では、この問題を確率的にモデル化してMAP推定により\bf{U}, \bf{V}を求める手法を提供している。ここでkは固定の値とするため、ユーザおよびアイテムのレイティング数にかかわらず、ランクkの近似行列を推定する。

本論文では、このPMFをベースに複数のランクkを重み付きで混合させたモデルを提案。具体的には下記の確率モデルを構築している。
\displaystyle\Pr({\bf R}|{\bf U},{\bf V},{\bf \alpha},{\bf \beta},\sigma^2)=\prod_{i=1}^m \prod_{j=1}^n \left[\sum_{k=1}^K \alpha_i^k \beta_j^k N(R_{i,j}|{\bf U}_i^k{{\bf V}_j^k}^T,\sigma^2) \right]^{ \mathbb{1}_{i,j}}
ここで、{\textbf{U}_i^k}はランクkの行列\textbf{U}^kのi行目(ユーザiの特徴ベクトル)、{\textbf{V}_j^k}はランクkの行列\textbf{V}^kのj行目(アイテムjの特徴ベクトル)を表す。また、\alpha_i^kはユーザiにおけるランクkモデルの重み、{\beta_j^k}はアイテムjにおけるランクkモデルの重み。これにより、ランクkのモデルがユーザi、アイテムjにマッチするなら\alpha_i^k,\beta_i^kは大きくなるはずであり、マッチしないなら\alpha_i^k,\beta_i^kは小さくなるはず(期待的には0になる)。これらの\bf{U},\bf{V},\bf{\alpha},{\bf \beta}をMAP推定により求めようというのが提案手法。

MAP推定のため事前分布を与える。\textbf{U}_i^k\textbf{V}_j^kの事前分布はそれぞれ、平均\textbf{0}、共分散\sigma_U^2 \textbf{I}\sigma_V^2 \textbf{I}ガウス分布としている(行列\textbf{U},\textbf{V}の事前分布はこれらの総乗)。また、{\bf \alpha}^k, {\bf \beta}^kはそれぞれスパースなベクトルにしたいので、\alpha_i^k,\beta_i^kの事前分布はラプラス分布とする(それぞれ、平均\mu_{\alpha},\mu_{\beta}、分散b_{\alpha},b_{\beta}というパラメータを持つ)。これらを積算することによって導き出される事後分布のlogにマイナスをつけた損失関数を最小化する。実際には、そのままの損失関数を最小化するのは難しいので、イェンセンの不等式で導き出される下限を最小化する。

\textbf{U},\textbf{V},{\bf \alpha},{\bf \beta}\sigma,\sigma_U^2,\sigma_V^2,\mu_{\alpha},\mu_{\beta},b_{\alpha},b_{\beta}を固定して確率的勾配法で解く。一方、 \sigma,\sigma_U^2,\sigma_V^2,\mu_{\alpha},\mu_{\beta},b_{\alpha},b_{\beta}も学習したいので、今度は\textbf{U},\textbf{V},{\bf \alpha},{\bf \beta}を固定して最尤推定する(こちらは、解析的に求まる)。これらの操作を収束するまで繰り返す。

MovieLensのデータセットNetflix Prizeのデータセットで実験。結果、kを固定するPMFに比べて高い推定精度を得られるし、他の手法と比べても有意な差が得られる。また、混合するモデルの個数(k個数)はたくさんかつ細かくした方が良いようだ(計算時間はかかるけど)。


以下2つの論文については、まだほとんど理解できていない。とりあえず、備忘録として残す。

Kernel Feature Selection via Conditional Covariance Minimization

カーネル法を使って特徴とラベル間の非線形関係をとらえて特徴抽出する手法。条件付き共分散作要素のトレースを基準にすることを提案。(相互)共分散作要素は2つのヒルベルト空間間の共分散を表しており、そこから条件付き共分散作要素が導き出される。
論文では条件付き共分散作要素のトレースがラベルとの独立性を評価する指標となると主張し、その理由を説明している。また、一般に特徴選択の最適化をまともに行おうとするとNP困難な問題になってしまうので、ここでは離散最適化から連続最適化へ問題を緩和している。また、それでも計算量的にしんどいので、いくつか計算量を抑えるための目的関数、制約条件の変換も行っている。

共分散作要素をもう少し理解してから再チャレンジが必要。今のところ共分散行列のヒルベルト空間への拡張と理解。日本語だとこの辺りが良さげ。

Unsupervised Image-to-Image Translation Networks

2つのドメインの画像セットがあり、一方のドメインの画像をもう一方のドメインの画像に変換する手法を提案した論文。例えば、低解像画像から高解像画像へ変換するタスクや、朝の風景画像を夜の風景画像に変換するタスクなどがあげられる。Unspervisedなので、学習用に2つのドメイン間で対応する画像のセットが与えられるわけではなく、ただ単にドメインAの画像セットとドメインBの画像セットが与えられるのみ。
これを実現するために2つのドメインのデータは同じlatent spaceを共有するという仮定を置く。この仮定をもとに、variational autoencoder(VAE)とgenerative adversarial network(GAN)によりドメイン変換を行う。

まず、VAEでは2つのドメインのインプットが同じlatent spaceに変換されるように、エンコーダ層を学習。この同じlatent spaceへ変換というのを実現するために、エンコーダ層の最後の数層とデコーダ層の最初の数層はドメイン間で重みを共有するという制約をつける(なぜこれで、うまくいくのかがイマイチわからない)。
次にGANは、adversarial discriminatorが画像生成(ドメイン変換)のネットワークによって生成された画像か実際の画像かを高精度に判定できるように学習させる。画像生成ネットワーク側はそれに負けじと高精度なドメイン変換を行う。
これら2つを両立させる目的関数(VAEとGANと制約項の和)を最適化することにより、ドメイン変換ネットワークを学習する。

VAEやGANは名前しか聞いたことなかったが、面白い枠組みだし応用範囲も広そう。まだまだブームは続く?

所感

久々にNIPSの論文を読んでみましたが、やはりニューラルネット関連が多いという印象。それでもニューラルネットばかりだけではなく、未だにカーネル法に関する論文も結構あるので(「kernel」で論文タイトルの検索かけると20件ヒット)まだまだこの分野も廃れていない、今後盛り返す可能性も十分ありと感じました。
あと、VAEやGANなど生成モデル+ニューラルネットの枠組みは、今後いろいろな分野に波及しそうな気がします(GAN+トピックモデルとかVAE+HMMとか、もう誰かやってる?)。

【Python】生存報告がてら昔書いたFFTのコードを晒す

お久しぶりです。ちゃんと生きてます。色々と立て込んでたり何だりでモチベーションが低下していました。このまま、フェードアウトしていくのもよろしくないと思いつつもネタがない・・・なんて思っていたところ、PCのデータを整理してたら、以前作成したFFT高速フーリエ変換)のプログラムを発見。てことで、このソースコードを今回のネタにしようかなと思います。
なお、フーリエ変換やDFT(離散フーリエ変換)についての解説は割愛します。また、FFTについても多少は説明しますが、DFTからFFTアルゴリズムの導出過程など詳しい説明は省きます。また、データのサンプル数は2のべき乗個に限定します。

FFTソースコード

ということで早速ソースコードです。まずは、モジュールのインポートとFFTの関数です。

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

def fft(x):
    n = x.shape[0]
    X = np.zeros(n ,  dtype = 'complex')
    if n > 2:
        y = x[0:n:2] #(1)
        z = x[1:n:2]

        Y = fft(y) #(2)
        Z = fft(z)

        Wn = np.exp(-1j * (2 * np.pi * np.arange(0 , n/2))/float(n)) #(3)
        X[0:n/2] = Y + Wn * Z
        X[n/2:n] = Y - Wn * Z
        return X
    else:
        X[0] = x[0] + x[1] #(4)
        X[1] = x[0] - x[1]
        return X

モジュールはいつも通り、numpyをインポート。後は、後々の動作確認で使うモジュールです。
関数fftが今回の肝であるFFTを行う関数ですが、見ての通り、再帰を使って書いてます。通常、スタックオーバーフローとかを避けるために、ビットリバース処理なんかをして繰り返し構文で書くのですが、今回はそこまでやってません。こっちの方が簡単だし直感的なので再帰で書く方法を採用しています。
ソースコード内のコメントにあるように、ポイントは4つ。
まず、(1)で元の信号xから偶数番目だけ取り出した信号yと奇数番目だけ取り出した信号zを作ります。
(2)では、fft関数を再帰呼び出ししてyとzのFFTの結果Y,Zを得ています。このY,Zを利用すると元の信号xのFFTの結果Xは次の式で、表されます。
  \left\{  \begin{array}{l} 
X(i) =  Y(i) + \exp \left( -j \frac{2 \pi i}{n} \right) Z(i) \: (0 \leq i < \frac{n}{2}) \\
X(\frac{n}{2} + i) = Y(i) - \exp \left( -j \frac{2 \pi i}{n} \right) Z(i) \: (0 \leq i < \frac{n}{2}) 
\end{array} \right.

(3)はこの処理を行っています。n点のDFTをこの式のようにn/2点のDFTの和に書きなおすことができるってのが、FFTの勘所です。
そして、(1)、(2)の再帰呼び出し処理を繰り返していくと、最終的にn = 2に行きつきます(元のデータサンプル数が2のべき乗の場合)。ここが再帰の終着点で、2点のDFTは(4)のように凄く簡単な計算で求められます。
これだけです。numpy使えば簡単にFFTのプログラムが書けます。ただし、numpyには元々FFTを行うメソッドが用意されているので、普通はわざわざ書く必要はありません。処理時間もこっちの方が断然速いです。

動作確認

てことで、一応動作確認。

def main():
    #テスト用の信号生成
    n =2**10
    x = (3 * np.cos(250* 2*np.pi * np.arange(0,n)/n) +
         2 * np.sin(250* 2*np.pi * np.arange(0,n)/n) +
         np.cos(30* 2*np.pi * np.arange(0,n)/n) )

    tic = time.clock()
    X = fft(x) #FFTの計算
    toc = time.clock()
    print toc - tic #処理時間表示

    #以下、結果グラフ表示
    plt.figure()
    plt.plot(X.real , "r")
    plt.xlim((0, n))
    plt.xticks(np.arange(0, 1024 , 100))
    plt.hold(True)
    plt.plot(X.imag , "b")
    plt.show()

テスト用の信号は、データのサンプル数が1024個で、周波数250のcos波、sin波、周波数30のcos波をそれぞれ適当なゲインをかけて足し合わせたものを使います。で、結果はこちら。

赤が実数部、青が虚数部です。うん、それっぽい。numpyに用意されているfftメソッドも同様の結果を吐き出します。また処理時間は、1024点なら30[ms]程度で処理できます。しかし、既存のメソッドの方は0.2[ms]程度。速い・・・

まとめ

今回はFFTのプログラム(再帰版)について話しました。上述のように、FFTはループ処理で書くことの方が多いようで、僕が知っている参考書では、大抵そちらの方法のサンプルプログラムが記述されています。しかし、それらは全部C言語で書かれているものなので、他の言語のサンプルコードが載っているものだと、どうなのか気になるところです。
あと、今回はサンプル点数が2のべき乗個である場合という最も基本的なFFTのみ記述しましたが、サンプル点数が素数の場合のFFTアルゴリズムが提案されていたりと応用的なFFTアルゴリズムも数々あります。これらを調べてみるのも、面白いかもしれません。

関係ないですが、転職により愛知から大阪に引っ越ししてきました。新天地からの更新でした。

マルチ○○学習まとめ

機械学習の分野では、マルチ○○学習という名付けられた枠組み・手法が色々提案されています。僕は、接頭辞が共通だと、すぐごっちゃになって何が何だか分からなくなってしまうので、ちょっと整理したいと思っていました。ということで、今回は「マルチカーネル」、「マルチビュー」、「マルチタスク」、「マルチラベル」、「マルチインスタンス」をメモ書き程度にまとめました。

マルチカーネル学習

複数の特徴ベクトル表現、または、類似の尺度が考えられるタスクに対し、複数のカーネル関数を組み合わせたカーネル学習を考えましょうという手法の枠組みです。例えば、画像データの場合、特徴ベクトルとしてbag-of-keypointsを構築したり、色情報を利用したりと、様々な類似尺度を考えることができます。そこで、それぞれの類似尺度に対応するカーネル関数を構築して組み合わせれば、ハッピーになれるだろうという考え方です。
多くの場合、次のようにカーネル関数の線形和を考えます。
k({\bf x} ,{\bf y}) = \sum_{i=1}^{N} \beta_i k_i({\bf x} ,{\bf y})Nカーネル関数の数)
このカーネル関数を学習の目的関数に組み込んで、最適な\betaを求めるのが、メジャーな手法といえるでしょう。必要なカーネル関数の重み(\beta)だけ大きくなるので、どれが学習に寄与する類似尺度かもわかるというメリットもあります。
有名なのは、大規模データをマルチカーネル+SVM等で学習する手法を提案したLarge Scale Multiple Kernel Learningという論文でしょうか。また、Dimensionality Reduction for Data in Multiple Feature Representationsでは、マルチカーネルを用いた次元削減手法を提案しています。他にも「multiple kernel learning」とかで検索すれば様々な論文が引っかかり、Multiple Kernel Learning Algorithmsというサーベイ論文なども見つかるので、勉強はしやすいと思います。

マルチビュー学習

1つのサンプルが様々な特徴量(ビュー)の組み合わせで構成されるデータセットを学習する問題の枠組み。例えば、動画データは画像シーケンスと音声信号という2つのビュー組み合わせとみなせますし、画像が含まれるウェブサイトなんかはテキストデータ+画像+リンク情報というように考えることができます。このようなデータセットに対し、効果的な学習手法を考えるのがマルチビュー学習の枠組みです。
手法としては、Combining labeled and unlabeled data with co-trainingで提案された共訓練(co-training)と呼ばれる手法が代表的でしょうか。これは、元々半教師あり学習の一手法を提案した論文で、マルチビューとみなせるデータセットならば(ただし、ビュー間はラベル情報のもと条件付き独立である必要あり)、各ビュー別々に分類器を学習し、ラベルなしデータをそれに判別させた結果を用いて更に分類器を学習するという反復手法がうまく行くことが示されています。
一方で、各々のビューが潜在的な部分空間を共有しているという仮定の下、その部分空間を見つけて写像してあげようというアプローチもあります。その中でも、正準相関分析(CCA)を利用した手法がよく知られています。他にも、フィッシャー判別分析を利用した手法、Neighborhood Preserving Projectionsを利用した手法などがあります。
また、上述のマルチカーネル学習の手法を適用した論文も多数あります。というか、マルチビュー学習の問題に、類似尺度が複数考えられる場合に効果的なマルチカーネル学習の手法を適用するのは自然なアプローチと言えます。もう少し言うと、マルチカーネルを適用すること自体、一種のマルチビュー学習とみなせるんじゃないでしょうか*1。そのため、マルチカーネル学習を適用して解いた問題を他のマルチビュー学習手法で解けば、もっと良くなるってこともあるかも・・・
なお、この項目を書くにあたって、A Survey on Multi-view Learningというサーベイ論文を参考にしました。

マルチタスク学習

複数の関連するタスク同士で情報を共有することにより、予測精度を向上させようという考え方の枠組み。特徴ベクトルやラベルの定義域はタスク間で共通の場合が多いです。考え方は転移学習と類似しているので、その辺りの論文も参考になります。両者の違いは、マルチタスク学習の方は、タスク間で協調し合うことで全てのタスクの精度を向上させようという目的に対し、転移学習はある目標のタスクがあって、その目標タスクの精度向上が目的である点で違います。
マルチタスク学習の例は、書き手が異なる学習データを用いた手書き文字認識です。学習データに本人が書いたデータが少ない場合、他人の書いた文字データを利用することで精度を向上させようという試みは、マルチタスク学習とみなせます(タスク1はAさんの手書き文字認識器を学習、タスク2はBさんの手書き文字認識器を学習・・・)。それ以外にも、テキストデータを分類する問題で分野の違うテキストデータを用いたり、シーンの異なる画像データで物体検出を行ったりといった場面でマルチタスク学習が適用できます。
手法としては、異なるタスクのサンプルを学習する際、そのサンプルがどれだけ適合するかという重みを付けるアプローチや、正則化項などを用いてタスク間でパラメータが類似するように学習するアプローチなどがあります。例えば、前者には、タスク間の関係が共変量シフトという状況下であると仮定し確率密度比を重みづけした手法が、Covariate Shift Adaptation by Importance Weighted Cross Validation等の論文で提案されています*2。後者も行列ノルムの正則化を利用し、各タスクのパラメータがスパース、かつ、0になる要素がタスク間で共通(jointly sparse)になるように学習する手法Robust visual tracking via multi-task sparse learning等々様々な手法があります。また、上述のマルチビュー学習とマルチタスク学習の問題設定を組み合わせたA Graph-based Framework for Multi-Task Multi-View Learningなんてものもあります。
サーベイ論文としては、転移学習のサーベイですがA Survey on Transfer Learningが有名でしょう。検索すれば日本語の資料なんかもたくさん出てきます。

マルチラベル学習

その名の通り、1つのサンプルに複数のラベルが割り当てれる分類問題の枠組み。例えば、小説のカテゴリ分けを考えた場合、「SF+ミステリー」、「ホラー+恋愛+ファンタジー」等のように、1つのジャンルに定められないことが多々あります。このような問題をマルチラベル学習では考えており、学習用サンプルに複数ラベル(カテゴリ)が割り当てられ、テスト用データの分類結果も複数ラベルが出力されるような手法が様々提案されています。
手法としては、大まかに分けると既存の学習手法をマルチラベルに拡張するアプローチと、マルチラベル問題をシングルラベルの問題に変換するアプローチの2つがあります。前者はSVMの拡張やAdaBoostの拡張、k近傍法の拡張など様々なものが提案されています。
後者の方もたくさんのアルゴリズムが提案されていています。例を挙げると、そのラベルが付いているか否かの2値分類問題に変換する手法(ジャンルがSFか否か、ミステリーか否かといった分類器をそれぞれ構築)、存在するラベルの組み合わせを1つのラベルとしてみなす手法(と言う1つのラベル、という1つのラベル、<ホラー+恋愛+ファンタジー>という1つのラベル)、複数ラベルが付いているサンプルは特徴ベクトル同一だがラベルが異なる別のサンプルとして学習する手法、などなど
マルチラベル学習の手法はほとんど知らなかったので、朱鷺の杜Wikiのマルチラベルの項目やここで紹介されていたチュートリアル論文Mining Multi-label Dataが非常に参考になりました。

マルチインスタンス学習

他の枠組みより多少問題設定がややこしいです。まず、用語として「bag」と呼ばれるものがあります。これは、複数のサンプル*3がひと塊りになったものです。マルチインスタンス学習では、学習データとして複数のbagと、それぞれのbagに正例か負例かというラベル情報が与えられます。ここで言う正例とは、bag内に1つでも正例に属するインスタンスがあること(負例インスタンスがいくつあってもよい)、対して負例はbag内全てのインスタンスが負例に属することを言います。bag内個々のインスタンスのラベルは分かりません。このような、学習データを与えられた上で、正例か負例かわからないbagを分類するというのがマルチインスタンス学習問題の枠組みです。
このような問題設定は、画像の分類問題に利用できます。例として、与えられた1枚の画像がビーチの画像か否かを判定するタスクを考えます。教師あり学習なので、学習用データとしてビーチの画像であると分かっている画像と、ビーチの画像ではないと分かっている画像が複数枚事前に与えられているとします。このような場合に、画像1枚の画像から領域ごとにわけて特徴ベクトルを抽出すると、1枚の画像=bag、抽出した各特徴ベクトル=インスタンスとみなせます。学習用データには画像(bag)に対して正例か負例かというラベル情報が割り当てられているので、このような問題設定は、まさにマルチインスタンス学習と言えます。それ以外にも、ある分子が薬として適切か否かを判定するタスク、文書データの分類などにも使われるようです。
アプローチとしては、与えられたbagセットからインスタンス単位での分類器を構築する、新規のbagに対しては、各インスタンスを分類器で判別してその結果を統合してbagがどちらに属するか判断するというアプローチ、bagレベルで学習するアプローチ(bag同士の距離を定めたりカーネル法を利用したり)、bagから1つの特徴ベクトルを構築(bag-of-keypointsなど)して学習するアプローチなどがあります。
まだちゃんと読んでいませんが、Multiple instance classification: Review, taxonomy and comparative studyというレビュー論文を参考にこの項目を書きました。

*1:例えば、カーネルの超パラメータを変えるだけでも、ある意味でビューを変えているというのが僕の考えです

*2:ここにあげた論文自体は、マルチタスク学習や転移学習を想定しているわけではありませんが、それらの一種とみなせると思います

*3:マルチインスタンス学習の枠組みでは「インスタンス」と呼びます。