MLPG：Maximum Likelihood Parameter Generationとは

#MLPG
Maximum Likelihood Parameter Generationについて書きます。

DeepLearningなど機械学習では

```math
\begin{align}
 \boldsymbol{y} &= f_\theta (\boldsymbol{x}) \\
 \theta &= argmin_\theta[ Loss(\boldsymbol{y}, \hat{\boldsymbol{y}}) ]
\end{align}
```
のように、入力$\boldsymbol{x}$から関数$f$を通して$\boldsymbol{y}$を出力するモデルを考えたりします。
このとき教師データ$\hat{\boldsymbol{y}}$と出力$\boldsymbol{y}$間のLossを最小化する関数パラメータ$\theta$を求めるのが学習ブロックです。

ただ、この特徴量$\boldsymbol{x}$は大抵が時刻$t$における特徴量です。
できれば前後の時刻との相関というか時間変化も捉えて学習した方が良さそうですよね？？？
（とは言え、MLPGでは「生成した特徴量の、時間方向の連続性が破綻してて音質が悪い！」のを防ぐことが目的で、音声の時間変化成分に含まれる個人性を捉えれるかどうかは別だと思います。せいぜい近傍10〜20msしか捉えられないので。）

なので例えば

```math
\begin{align}
\boldsymbol{x} &= [\boldsymbol{x_1}^T, \boldsymbol{x_2}^T, \cdots, \boldsymbol{x_T}^T]^T \\
&(\boldsymbol{x_t} = [x_{1t}, x_{2t}, \cdots, x_{Dt}]^T)
\end{align}
```

という静的特徴量$\boldsymbol{x}$に動的特徴量を含んだ特徴量$\boldsymbol{X}$

```math
\begin{align}
\boldsymbol{X} &= [\boldsymbol{X_1}^T, \boldsymbol{X_2}^T, \cdots, \boldsymbol{X_T}^T]^T \\
&(\boldsymbol{X_t} = [\boldsymbol{x_t}^T, \Delta\boldsymbol{x_t}^T])
\end{align}
```

を扱います。
それで

```math 
\begin{align}
\boldsymbol{Y} &= f_\theta (\boldsymbol{X}) \\
\theta &= argmin_\theta[ Loss(\boldsymbol{y}, \hat{\boldsymbol{y}}) ] \\
       &= argmin_\theta[ Loss(\boldsymbol{RY}, \hat{\boldsymbol{y}}) ]
\end{align}
```

というように、静的・動的特徴量$\boldsymbol{X}$を入力したら出力$\boldsymbol{Y}$を求めるモデルを考えます。
そのモデルのパラメータ$\theta$ は $\boldsymbol{RY} (= y)$ と $\hat{\boldsymbol{y}}$のロスを最小化するように学習します。

ここで重要なのは静的・動的特徴量と静的特徴量間の変換行列$\boldsymbol{R}$をどうやって産むかということです。

次節以降、もう少し具体的に必要な行列とかに関する算出方法のメモとコードを記します。
#
#
#####################################################################
#
#
#ちょいと数学
ここでやるのは$\boldsymbol{X} = \boldsymbol{Wx}$ なる変換行列$\boldsymbol{W}$ と
$\boldsymbol{x} = \boldsymbol{RX}$ なる変換行列$\boldsymbol{R}$を求めることです。


#X = Wx について
※今回は$\Delta\boldsymbol{x}$までしか扱いません。記事や図、コードを書くのが楽なので。さらに$\Delta\Delta\boldsymbol{x}$も考慮したって理屈は同じです。

では本題

![スクリーンショット 2021-03-03 23.25.01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/156754/f3ca84b1-6e15-770b-8f7a-81437872b1ab.png)

>記事では$\boldsymbol{x}$とか$\boldsymbol{y}$とか入り混じってますが、
大文字と小文字間の変換なので意味は変わりません。ごめんなさい。

図にしてみました！

この行列$\boldsymbol{W}$が言いたいところは

```math
\Delta\boldsymbol{y_t} = -0.5\boldsymbol{y_{t-1}} + 0.5\boldsymbol{y_{t+1}}
```

だよっていうだけです。
#
###x = RX について

```math
\boldsymbol{R} = (\boldsymbol{W}^T\Sigma^{-1}\boldsymbol{W})^{-1}\boldsymbol{W}^T\boldsymbol{\Sigma}^{-1}
```

$\boldsymbol{\Sigma}$は$\boldsymbol{X}$の共分散行列です。

#
###################################################
#
#「W」算出の実装しまーす

```python
import scipy.sparse
import numpy as np

def computeW(T, D):
  for t in range(T):
    w0 = scipy.sparse.lil_matrix((D, D * T))
    w1 = scipy.sparse.lil_matrix((D, D * T)) 
    w0[0:, t * D:(t + 1) * D] = scipy.sparse.diags(np.ones(D), 0)
    
    if t >= 1:
      tmp1 = np.zeros(D)
      tmp1.fill(-0.5)
      w1[0:, (t - 1) * D:t * D] = scipy.sparse.diags(tmp1, 0)

    if t < T-1:
      tmp1 = np.zeros(D)
      tmp1.fill(0.5)
      w1[0:, (t + 1) * D:(t + 2) * D] = scipy.sparse.diags(tmp1, 0)
            
        
     W_t = scipy.sparse.vstack([w0, w1])
     
     if t == 0:
       W = W_t
     else:
       W = scipy.sparse.vstack([W, W_t])

  W = scipy.sparse.csr_matrix(W) 
  return W                          
```

では解説しまーす。
方針は、2D行DT列の行列を１まとまりとして考えます。
最初のD行が静的特徴量、次のD行が動的特徴量に関連します。

```python
w0 = scipy.sparse.lil_matrix((D, D*T))
```

これは、D行DT列のゼロ行列w0を生成します。時刻tにおける静的特徴量に関連する行ブロックです。
次のw1は動的特徴量の算出を担当します。

例えば音声だとDはメルケプストラムなどで３９次元とかありますし、Tは発話単位のフレーム数で１０００次元くらいあるので素直にnumpy.arrayを使うと１５２１０００個の要素を持ちます。大変な事ですね。

なのでscipyのスパース行列を利用します。
 - どのインデックスに０以外の値を持つか
 - そこにどの値が入ってるか
という情報だけを格納します。なので今回みたいな０が大半を占める行列を扱う時にはもってこいなわけです。

```python
w0[0:, t*D:(t+1)*D] = scipy.sparse.diags(np.ones(D), 0)
```

さっき乗せた図の各マスはD行D列の対角行列で、対角要素はマス目の数字です。
まず、「１」が入ってるところに注目してください。
時刻tにおけるブロック（2D行TD列）の上段で、t*D:(t+1)*D列に「１」が割り当てられてます。
つまり、そこのマス（D行D列のサブ行列）が単位行列になります。
上記のコードはそれを意味してます。

```python
    if t >= 1:
      tmp1 = np.zeros(D)
      tmp1.fill(-0.5)
      w1[0:, (t - 1) * D:t * D] = scipy.sparse.diags(tmp1, 0)
```
今度は「-0.5」に注目します。
まぁやりたいことは、「１」のときと同じなので、さっきの図を見ながら呼んでくれればわかるかと思います。

```python
    if t < T-1:
      tmp1 = np.zeros(D)
      tmp1.fill(0.5)
      w1[0:, (t + 1) * D:(t + 2) * D] = scipy.sparse.diags(tmp1, 0)
```

今度は「0.5」です。詳細は今までと同じ感じです。

次、

```python
 W_t = scipy.sparse.vstack([w0, w1])
     
     if t == 0:
       W = W_t
     else:
       W = scipy.sparse.vstack([W, W_t])

  W = scipy.sparse.csr_matrix(W) 
```

w0とv1を垂直結合します。最初に「2D行DT列を１まとまりとして」と言いました。
これで１まとまりのブロックの完成です。

あとはそれをT回繰り返してがっちゃんこすれば
2DT行DT列の行列$\boldsymbol{W}$の完成です。やったー！！！

#
#
###########################################################
#
#
#「R」算出の実装しまーす

```python
def covar_matrix(Y, T, D):
  Y_mat = Y.reshape(T, int(len(Y)/T))
  S_vec = np.var(Y_mat, axis=0)
    
  M = len(S_vec)
  diagonal_t = 1.0/S_vec
  S_inv = np.zeros((T, M, M))
  for t in range(T):
    S_inv[t] = np.diag(diagonal_t)

  S_inv = scipy.sparse.block_diag(S_inv, format='csr')
  return S_inv
```

まず$\Sigma^{-1}$を求めます。共分散行列の逆行列です。ただし簡易化するために今回は対角行列化したただの分散行列にします。
対角行列の逆行列は、対角成分の逆数をとれば良いだけなのでめっちゃ簡単です。

次、

```python
def computeR(W, S_inv):
  R1 = W.T.dot(S_inv).dot(W)
  R1 = scipy.sparse.linalg.inv(R1)
  R2 = W.T.dot(S_inv)
  R = R1.dot(R2)

  return R
```

```math
\boldsymbol{R} = (\boldsymbol{W}^T\Sigma^{-1}\boldsymbol{W})^{-1}\boldsymbol{W}^T\boldsymbol{\Sigma}^{-1}
```
を計算してます。


#
#
###############################################################
#
#
#全体コード
```python
import scipy.sparse
import numpy as np

def computeW(T, D):
  for t in range(T):
    w0 = scipy.sparse.lil_matrix((D, D * T))
    w1 = scipy.sparse.lil_matrix((D, D * T)) 
    w0[0:, t * D:(t + 1) * D] = scipy.sparse.diags(np.ones(D), 0)
    
    if t >= 1:
      tmp1 = np.zeros(D)
      tmp1.fill(-0.5)
      w1[0:, (t - 1) * D:t * D] = scipy.sparse.diags(tmp1, 0)

    if t < T-1:
      tmp1 = np.zeros(D)
      tmp1.fill(0.5)
      w1[0:, (t + 1) * D:(t + 2) * D] = scipy.sparse.diags(tmp1, 0)
            
        
     W_t = scipy.sparse.vstack([w0, w1])
     
     if t == 0:
       W = W_t
     else:
       W = scipy.sparse.vstack([W, W_t])

  W = scipy.sparse.csr_matrix(W) 
  return W   


def covar_matrix(Y, T, D):
  Y_mat = Y.reshape(T, int(len(Y)/T))
  S_vec = np.var(Y_mat, axis=0)
    
  M = len(S_vec)
  diagonal_t = 1.0/S_vec
  S_inv = np.zeros((T, M, M))
  for t in range(T):
    S_inv[t] = np.diag(diagonal_t)

  S_inv = scipy.sparse.block_diag(S_inv, format='csr')
  return S_inv             


def computeR(W, S_inv):
  R1 = W.T.dot(S_inv).dot(W)
  R1 = scipy.sparse.linalg.inv(R1)
  R2 = W.T.dot(S_inv)
  R = R1.dot(R2)

  return R


def super2static(Y, W, T, D):
  S_inv = covar_matrix(Y, T, D)
  R = limitationR(W, S_inv)
  return R.dot(Y)     


def static2super(y, T, D):
  y_ = y.reshape(T*D, 1)
  W = computeW(T, D)
  Y = W.dot(y_)
  return Y, W
```

#
#######################################################################
#
#DeepLearningに組み込む
よく見るのはMSE基準による最適化です。

```math
  L = \frac{1}{T}(\boldsymbol{y} - \hat{\boldsymbol{y}})^T(\boldsymbol{y} - \hat{\boldsymbol{y}})
```

単に教師データ$\boldsymbol{y}$と生成データ$\hat{\boldsymbol{y}}$の誤差を
フレーム方向に（データ数方向に）平均してます。

ただこのとき、$\boldsymbol{y}$はある関数$f$を介して入力$\boldsymbol{x}$から生成されたデータとします。

$\boldsymbol{x}$はミニバッチ学習などを採用していれば、無造作に選ばれたN個のデータからなるので
n番目のデータとn+1番目のデータに繋がりはありません。

#
#######################################################################
#
#MGE基準の最適化
こいつはMinimum-Generation-Errorのことで、

```math
\begin{align}
  L &= \frac{1}{T}(\boldsymbol{RY} - \hat{\boldsymbol{y}})^T(\boldsymbol{RY} - \hat{\boldsymbol{y}}) \\
  \boldsymbol{Y} &= f(\boldsymbol{x})
\end{align}
```

のように、静的・動的特徴量を出力するモデルに対して、本来欲しい静的特徴量間の誤差を用いて最適化します。
※もし音声変換ならモデルの入力は$\boldsymbol{X}$など

なので、この状況下では
- ミニバッチ単位では無く、１塊のデータ単位での学習（例えば１発話など）
　※その都度$\boldsymbol{Y}$から制約行列$\boldsymbol{R}$を推定
になります。
これが「Generation」の由来でしょうかね。

今回はscipyのスパース行列でMLPGを実装しましたが、DeepLearningをPyTorchで実装するならMLPGもそれに対応させる必要があります。
それは疲れたのでまたいつかやります。

#
#################################################################
#
#まとめ
ということで今回はMLPGとMGE基準についてメモを書きました。
細かいミスはごめんなさい。
勘違いや誤りがあれば教えていただけると幸いです。。。