特異スペクトル分解


今回の記事も、うすぼんやりとメモ程度です。

http://www.jspf.or.jp/Journal/PDF_JSPF/jspf2020_05/jspf2020_05-262.pdf

このペーパーを見て実装してみました。どんなご利益があるのかというと、データを複数の成分に分解してくれます。今回は、ノイズの乗った時系列のノイズ取りに使ってみました。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/377651/ae74c845-c52a-5716-4631-145571b1a345.png)

青線が元系列です。端っこの方で、結構なノイズが載っていますが点線の真の値に、結構な精度で一致する推定値が得られています。すごーい。ちなみに、これは分解した結果の１～５番目の成分を足し合わせたものを載せています。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/377651/bbeb6046-59dd-6a9a-94e3-09eef640dcb5.png)

これが成分ごとのプロットです。SVDを利用して時系列を分解して、成分ごとに分離するという手法です。なんか、PCAとかそんなんの親戚だと思われます。たぶん、スペクトルがどうとか、深遠な理論的根拠があるものだとは思いますが、私は時系列を何らかの個別成分＋ノイズに分解する方法に興味があっただけなので、面白半分で実装してみました。たぶん、ちゃんと動いていると思う。

SVDは固有値分解っぽい感じものですが、任意の行列に対して適用可能。任意の行列を直行する規定ベクトルと係数に分解してくれるというものだそうです。何らかの行列Aがあったとき、u@Σ@v.Tのような形に分解してくれます（np.linalg.svd関数の戻り値もこの順）。私たち数学音痴は思いついた数学者の人とnumpy（及び、下回りの数値計算ライブラリ、LINPACKかなんかだっけ？）を実装した人々に感謝しつつ、うすぼんやりの理解で利用するだけさせてもらいましょう。

ちなみに、np.linalg.svdの実装では、真ん中のΣに相当する奴が、行列じゃなくてベクトルで出てきます。また、uとかvも縦ベクトルとか、横ベクトルとか関係なしに出てくるので、戻すのに手間取りました。この辺の整合性はちゃんととって欲しいなぁ（まぁ、Σに関しては、元の行列の列数分、対角成分が並んでいるだけなので、ベクトルだけでも良いと思うけど）。

まーまー、論文の結果に近い結果が出ました。もちろん、参考にした論文は、こんな適当な感じではなくて、色々と評価を行って（ほかのフィルタにくらべて、振動の周期にバイアスがでてないよとか、なんとか）いますので、興味がある方は読んでみたらよいかもしれません。ガチ理系の論文ですが、でんでんむしもにょろにょろも出てこないので、数学音痴でも読めます。

以下、コードを晒します。

```python
#
# for time series, np.array(x(1), x(2), ..., x(T-1))
# return hankel matrix X(w,n).
# t is current time and it must be larger than W
#
# Hankel matrix:
# s(t-1) = (x(t-w),...,x(t-1)).T
# H(t) = [s(t-n),...,s(t-1)]
#
# 時系列からトラジェクトリー行列を作る。
def trajectory(x,t,w,n):    
    W = w + n - 1
    assert t <= len(x) and t >= W
    s = lambda t: x[t-w:t]
    return np.array([s(t-k) for k in range(n,-1,-1)])

# make timeseries from trajectory matrix
# トラジェクトリー行列を時系列に戻す。
def untrajectory(X):
    return np.concatenate([X[:,0],X[len(X)-1,1:]])

# return u,S,v
# 時系列→トラジェクトリー行列→SVDという操作をまとめる。
def trajectory_decomp(x,w):
    T = len(x)
    H = trajectory(x,T,w,T-w)
    return np.linalg.svd(H, full_matrices=True)

# これがめんどくさかったわ、マジで
# SVD分解すると、行列はU*S*V.Tという形に分解される。
# だからS[i]*u[i]@v[i].Tとやれば、分解された行列が取れると思うでしょう？
# それが、そうでもないんだな
def trajectory_restore(u, S, v, i):
    return S[i]*np.array([u.T[i]]).T @ np.array([v[i]])

# 元のペーパーに証明が載っている。
# ある行列Xに対して、|X-H|^2が最小となるようなHankel行列Hを計算する。
# compute hankel matrix H, where argmin(H,|X-H|^2)
def hankelize(X):
    N, M = X.shape
    H = np.zeros((N,M))
    As = [[] for i in range(N+M)]
    # count number of s = i + j
    for i in range(N):
        for j in range(M):
            As[i+j] += [(i,j)]
    # compute
    for i in range(N):
        for j in range(M):
            L = len(As[i+j])
            for k, l in As[i+j]:
                H[i][j] += X[k][l] / L
    # return matrix
    return H
            
# http://www.jspf.or.jp/Journal/PDF_JSPF/jspf2020_05/jspf2020_05-262.pdf
# 論文に載っているテスト用のデータ
f=1
a=20
b=0.8
c=0.08
d=0.1
e=1.8
ts=np.arange(0,30,30/150)
noise=np.random.normal(0.0, 3.0, 150)
y_hat=a*np.sin(2*np.pi*f*np.power(ts,b))*np.exp(-c*ts)+d*np.power(ts,e)
y=y_hat+noise

# 分解！！！！
u,S,v=trajectory_decomp(y,28)

# 成分ごとに行列に戻す
x0=trajectory_restore(u, S, v, 0)
x1=trajectory_restore(u, S, v, 1)
x2=trajectory_restore(u, S, v, 2)
x3=trajectory_restore(u, S, v, 3)
x4=trajectory_restore(u, S, v, 4)
x5=trajectory_restore(u, S, v, 5)

# ハンケル行列化して、時系列に戻す
X_hat=untrajectory(hankelize(x0+x1+x2+x3+x4))

```

この論文の手法のミソは、SVDでも特異スペクトル分解でもなくて、たぶん、ばらした行列を時系列に戻すのに「ハンケル行列化」している部分だと思われます。実際、特異スペクトル分解について記載された文献は山とあり、それほど珍しい手法でもなかったはずです（ICAとかそうじゃなかったっけ？）。計算としては至って単純なのですが、ハンケル行列Hを想定して、誤差をとりあえず計算し、それを最小にするにはどのような計算をすればよいのか？という論理構成で至ってシンプルに導出されているところが素敵でした。

まー、こんなのだれの参考にもならんだろうね。