甲斐性なしのブログ

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

再構成誤差における平均ベクトルとの差はどのように評価されるか

はじめに

異常検知などのタスクでは、あらかじめ取得した学習データにより形成される主部分空間にテストデータを射影し、元のテストデータと射影後のデータの距離を異常度とすることがしばしば行われる。これを再構成誤差という(詳細は下記書籍を参照)。

もう少し具体的に書くと、学習データを {\bf X} = \left\{{\bf x}_1, \cdots {\bf x}_N | {\bf x}_i \in {\bf R}^d, i = 1, \cdots, N \right\}、テストデータを{\bf x}_{*} \in {\bf R}^dとすると、学習データの平均ベクトルと共分散行列は

 \displaystyle {\bf m} = \frac{1}{N} \sum_{i = 1}^N {\bf x}_i \tag{1}  \displaystyle {\bf C} = \frac{1}{N} \sum_{i = 1}^N ({\bf x}_i - {\bf m}) ({\bf x}_i - {\bf m})^T \tag{2}

として推定される。そして、この{\bf C}固有値分解{\bf C} = {\bf U} {\bf \Lambda} {\bf U}^Tにより、d個の固有ベクトル{\bf U} = \left[{\bf u}_1, \cdots , {\bf u}_d \right] を得るが、このうち固有値が大きい方からp個だけ対応する固有ベクトルを取ってくる。このp個の固有ベクトルが射影先の主部分空間の基底である。これを列ベクトルに並べた行列を{\bf U}_p = \left[{\bf u}_1, \cdots , {\bf u}_p \right] とすると射影行列は{\bf U}_{p} {\bf U}_{p}^{T}なので、テストデータ{\bf x}_{*}の再構成誤差は、

 \displaystyle a({\bf x}_{*} ) = \parallel ({\bf x}_{*} - {\bf m}) - {\bf U}_{p} {\bf U}_{p}^{T} ({\bf x}_{*} - {\bf m}) \parallel_{2}^2 \tag{3}

として計算される。さて、式(3)では平均ベクトル {\bf m}で引いているが、このテストデータと平均ベクトルとの差({\bf x}_{*}-{\bf m})は再構成誤差にどのように効いてくるだろうか?直感的には{\bf x}_{*}{\bf m}から遠ければ遠いほど再構成誤差は大きくなりそうだが*1実はそうではない、というのが今回のお話。結論を言ってしまうと、

  • テストデータ{\bf x}_{*}が主部分空間の元である場合、平均ベクトルとは無関係に一定値になる
  • テストデータ{\bf x}_{*}が主部分空間の元でない場合、平均ベクトルとの距離が大きくても再構成誤差は大きくならないことがある
  • 主部分空間を平行移動させてできるアフィン部分空間内のテストデータは、全て同じ再構成誤差になる

である。以下、これを詳しく見ていく。

テストデータ{\bf x}_{*}が主部分空間の元である場合

このケースはたまたま{\bf x}_{*}が主部分空間の基底の線形結合で表せる、つまり {\bf x}_{*} = {\bf U}_p {\bf a}となる{\bf a} \in {\bf R}^{p}が存在するケースである。 式(3)を少し変形すると、

 \displaystyle
 a({\bf x}_{*} ) =\parallel ({\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf x}_{*} ) - ({\bf m} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} ) \parallel_{2}^2 \tag{4}

となる。式(4)の中で({\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf x}_{*} )に着目し{\bf U}_{p}^{T} {\bf U}_{p} = {\bf I}であることに注意すると、

 \displaystyle
 {\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf x}_{*}   = {\bf x}_{*} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf U}_{p} {\bf a} = {\bf x}_{*} - {\bf U}_{p} {\bf a} =   {\bf x}_{*}  -  {\bf x}_{*}  = {\bf 0}

なので、再構成誤差は {\bf x}_{*}によらない一定値\parallel - ({\bf m} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} ) \parallel_{2}^2になる。従って、 {\bf x}_{*}が主部分空間の元の場合 {\bf m}からどれだけ離れていようが近くにあろうが全く関係ない。

このことを以下のpythonコードで確かめてみる。

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

if __name__ == '__main__':
    # 学習データの生成とプロット
    np.random.seed(1234)
    N = 1000
    mu = [3, 10]
    Sig = [[8, 4],[4, 7]]
    X = np.random.multivariate_normal(mu, Sig, N).T
    plt.scatter(X[0, :], X[1, :], color = 'r', label="Train Data")

    # 学習データのサンプル平均推定とそのプロット
    m = np.average(X, axis = 1)
    plt.scatter(m[0], m[1], color = 'c', 
                s = 70, marker = "*", label = "Train Sample Mean")

    # 主成分分析で主部分空間推定とそのプロット
    Z = X - m[:, np.newaxis]
    C = Z @ Z.T/N
    lam, U = np.linalg.eig(C)
    Up = U[:, 0][:, np.newaxis]
    x = [-100 * U[0, 0], 0, 100 * U[0, 0]]
    y = [-100 * U[1, 0], 0, 100 * U[1, 0]]
    plt.plot(x, y, color = 'b', label="Principal Axes")
    
    # 学習データの平均から離れたテストデータの再構成誤差
    test_x1 = Up.squeeze() * -20
    plt.scatter(test_x1[0], test_x1[1], color = 'g', 
                s = 70, marker = "^", label = "Test Data 1")
    a1 = np.linalg.norm((test_x1 - m).T - Up @ Up.T @ (test_x1 - m)) ** 2
    print("Test Data 1 Reconstruct Error", a1)

    # 学習データの平均に近いテストデータの再構成誤差
    test_x2 = Up @ Up.T @ m 
    plt.scatter(test_x2[0], test_x2[1], color = 'g',
                s = 70, marker = "v", label = "Test Data 2")
    a2 = np.linalg.norm((test_x2 - m).T - Up @ Up.T @ (test_x2 - m)) ** 2
    print("Test Data 2 Reconstruct Error", a2)

    plt.xlim(-25, 25)
    plt.ylim(-25, 25)
    plt.legend()
    plt.show()

このコードにより、以下の図のような2次元のデータを生成し、主部分空間上にある2つのテストデータの再構成誤差を算出している。 f:id:YamagenSakam:20200415224630p:plain

そして、Test Data 1とTest Data 2の再構成誤差の計算結果は以下の通り。

Test Data 1 Reconstruct Error 28.670763684482974
Test Data 2 Reconstruct Error 28.670763684482964

Test Data 1はTest Data 2に比べて明らかに学習データの平均値から遠いにもかかわらず再構成誤差は同じ値になっている。

テストデータ{\bf x}_{*}が主部分空間の元でない場合

一般の場合を考える。式(3)を展開すると以下を得る。

 \displaystyle 
\begin{aligned}
a({\bf x}_{*} ) &= ({\bf x}_{*} - {\bf m})^T ({\bf I}- {\bf U}_{p} {\bf U}_{p}^{T})^T ({\bf I}- {\bf U}_{p} {\bf U}_{p}^{T}) ({\bf x}_{*} - {\bf m}) \\
&= ({\bf x}_{*} - {\bf m})^T ({\bf I}- {\bf U}_{p} {\bf U}_{p}^{T}) ({\bf x}_{*} - {\bf m}) \\
&=  ({\bf x}_{*} - {\bf m})^T({\bf x}_{*} - {\bf m}) - ({\bf x}_{*} - {\bf m})^T{\bf U}_{p} {\bf U}_{p}^{T}({\bf x}_{*} - {\bf m}) \\
&=  {\bf x}_{*}^T ({\bf I} -{\bf U}_{p} {\bf U}_{p}^T ){\bf x}_{*} - 2 {\bf m}^T  ({\bf I} -{\bf U}_{p} {\bf U}_{p}^T ){\bf x}_{*} +{\bf m}^T ({\bf I} -{\bf U}_{p} {\bf U}_{p}^T ){\bf m}
\end{aligned}
\tag{5}

第1項はテストデータの主部分空間からのはみだし具合だけで決まるので平均ベクトル{\bf m}は無関係、第3項は定数、従って第2項に着目する。この第2項が大きければ再構成誤差としては小さく、逆に第2項が小さければ再構成誤差は大きくなる。({\bf I} -{\bf U}_{p} {\bf U}_{p}^T )は主部分空間の直交補空間への射影行列である。すなわち、第2項はテストデータ{\bf x}_{*}を直交補空間へ射影したベクトルと平均ベクトルの内積を計算している。従って、どんなにテストデータが平均ベクトルから離れていても直交補空間へ射影されてから内積を計算するので射影された先が同じなら第2項は同じ値になる。また内積なので、射影後のベクトルと平均ベクトルの方向が重要で、方向が同じであるほど第2項は大きな値となる(再構成誤差は小さくなる)。

このことを確かめるために以下のようなpythonコードで実験。なお、主部分空間推定までは上のコードと共通なので省略。

    # 主部分空間推定までは上のコードと同じなので省略

    # 学習データの平均から離れているが学習データの平均と同じ方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x1 = (Up.squeeze() * -20) + 5 * U[:, 1]
    plt.scatter(test_x1[0], test_x1[1], color = 'g', 
                s = 70, marker = "^", label = "Test Data 1")
    a1 = np.linalg.norm((test_x1 - m).T - Up @ Up.T @ (test_x1 - m)) ** 2
    print("Test Data 1 Reconstruct Error", a1)

    # 学習データの平均から離れているが学習データの平均と逆方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x2 = (Up. squeeze() * -20) - 5 * U[:, 1]
    plt.scatter(test_x2[0], test_x2[1], color = 'g', 
                s = 70, marker = "v", label = "Test Data 2")
    a2 = np.linalg.norm((test_x2 - m).T - Up @ Up.T @ (test_x2 - m)) ** 2
    print("Test Data 2 Reconstruct Error", a2)

    # 学習データの平均と近く学習データの平均と同じ方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x3 = Up @ Up.T @ m + 5 * U[:, 1]
    plt.scatter(test_x3[0], test_x3[1], color = 'g', 
                s = 70, marker = ">", label = "Test Data 3")
    a3 = np.linalg.norm((test_x3 - m).T - Up @ Up.T @ (test_x3 - m)) ** 2
    print("Test Data 3 Reconstruct Error", a3)

    # 学習データの平均と近く学習データの平均と逆方向に
    # 主部分空間から飛び出ているテストデータの再構成誤差
    test_x4 = Up @ Up.T @ m - 5 * U[:, 1]
    plt.scatter(test_x4[0], test_x4[1], color = 'g', 
                s = 70, marker = "<", label = "Test Data 4")
    a4 = np.linalg.norm((test_x4 - m).T - Up @ Up.T @ (test_x4 - m)) ** 2
    print("Test Data 4 Reconstruct Error", a4)

    # 以降のプロットする部分も上のコードと同じなので省略

この実験では以下のように、主部分部分空間からそれと直交する方向に平行移動させた4つのテストデータで実験。 f:id:YamagenSakam:20200415230950p:plain

そして、この再構成誤差の計算結果は以下の通り。

Test Data 1 Reconstruct Error 0.12567643599399253
Test Data 2 Reconstruct Error 107.21585093297189
Test Data 3 Reconstruct Error 0.12567643599399056
Test Data 4 Reconstruct Error 107.21585093297192

このように、Test Data 1とTest Data 3は平均ベクトルと同じ方向に同じ量だけ主部分空間から平行移動させたため、再構成誤差としてかなり小さい同じ値になっている(Test Data 3はほとんど平均ベクトルと一致しているので再構成誤差は小さくなるのは当然だが、平均ベクトルから離れているTest Data 1もそれと全く同じ値になっている)。逆にTest Data 2とTest Data 4については、平均ベクトルと逆方向に主部分空間から平行移動させているため、再構成誤差は大きな値になっている(Test Data 4はTest Data 1よりも平均ベクトルとの距離という意味では小さいのに、再構成誤差はTest Data 4の方が大きい)。

ちなみに平均値{\bf m}が主部分空間の元である場合、直交補空間の元との内積{\bf 0}になるので、平均値とは無関係に元のデータの主部分空間からのはみ出具合のみで評価される。 これは式(4)からも明らか。

再構成誤差の等高線

ここまでの結果を見ると主部分空間からどれだけ、どの方向に平行移動させたかが再構成誤差に効いてくるように思える。これを確かめるため再構成誤差の等高線を引いてみる。そのコードと実行結果は以下の通り。

    # 主部分空間推定までは上のコードと同じなので省略

    # 再構成誤差の等高線の計算と描画
    x1 = np.arange(-35, 35, 1) 
    x2 = np.arange(-35, 35, 1)    
    X1, X2 = np.meshgrid(x1, x2)
    A = np.zeros(X1.shape)
    for i in range(X1.shape[0]):
        for j in range(X1.shape[1]):
            test_x = np.array([X1[i, j], X2[i, j]])
            # 2乗があると等高線の間隔が広がりすぎてわかりにくいので省く
            A[i, j] = np.linalg.norm((test_x - m).T - Up @ Up.T @ (test_x - m))
    plt.contour(X1, X2, A, levels=50, cmap="viridis_r")
    plt.colorbar()

    # 以降のプロットする部分も上のコードと同じなので省略

f:id:YamagenSakam:20200415233632p:plain

これを見ると確かに主部分空間と平行な線上に同じ再構成誤差の値が乗っているように見える。実際、主部分空間を平行移動させた平面(アフィン部分空間)上のテストデータは同じ再構成誤差になる。以下これを示す。

任意の{\bf x}_{*}は、主部分空間上の元{\bf U}_p {\bf a}とその直交補空間の元{\bf b}を使うと、

\displaystyle {\bf x}_{*} = {\bf U}_p {\bf a} + {\bf b} \tag{6}

と表せる。{\bf b}を加えるということは、主部分空間上の元を{\bf b}だけ平行移動させて主部分空間から飛び出させることに等しい。ここで式(6)を式(3)に代入し{\bf U}_{p}^T {\bf b} = {\bf 0}であることに注意して展開すると、

 \displaystyle
\begin{aligned}
a({\bf x}_{*} ) &= \parallel ( {\bf U}_p {\bf a} + {\bf b} - {\bf m}) - {\bf U}_{p} {\bf U}_{p}^{T} ( {\bf U}_p {\bf a} + {\bf b} - {\bf m}) \parallel_{2}^2 \\
&=  \parallel ( {\bf U}_p {\bf a} + {\bf b} - {\bf m}) -  ({\bf U}_{p} {\bf U}_{p}^{T} {\bf U}_p {\bf a} + {\bf U}_{p} {\bf U}_{p}^{T}{\bf b} - {\bf U}_{p} {\bf U}_{p}^{T} {\bf m}) \parallel_{2}^2 \\
&=  \parallel  {\bf U}_p {\bf a} + {\bf b} - {\bf m} -  {\bf U}_p {\bf a} +  {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} \parallel_{2}^2 \\
&= \parallel {\bf b} - {\bf m}+  {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} \parallel_{2}^2 \\
&= {\bf b}^T {\bf b}  - 2 {\bf b}^T{\bf m} + {\bf m}^T {\bf m} + 2 {\bf b}^T  {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} - 2{\bf m}^T {\bf U}_{p} {\bf U}_{p}^{T} {\bf m} +  {\bf m}^T {\bf U}_{p}{\bf U}_{p}^{T} {\bf m}\\
&= \parallel {\bf b} - {\bf m} \parallel_{2}^2- {\bf m}^T {\bf U}_{p} {\bf U}_{p}^{T} {\bf m}
\end{aligned}
\tag{7}

となり、再構成誤差は平行移動成分{\bf b}と平均ベクトル{\bf m}の距離のみに依存することがわかる。すなわち、平行移動成分{\bf b}が一定なら、主部分空間上の元をどれを取ってきても再構成誤差は同じである。つまりそれは、主部分空間を{\bf b}だけ平行移動させたアフィン部分空間上の元の再構成誤差は、全て等しい値になることを意味する。

まとめ

ここで示したように平均ベクトルとテストデータの距離が大きくても再構成誤差は小さくなりうる。再構成誤差を使って異常検知などを行う場合、このことを頭の片隅にでも置いておかないと落とし穴にはまるかもしれない。

*1:というか、異常検知の枠組みだとそうなって欲しいケースが多いと思う