0. はじめに
前回の記事では、Grevilleの平滑化について解説した。この記事は続編である。Grevilleの平滑化で用いる係数$L_t$と$a_t$がどのようにして得られるかということについて解説する。
1. 数列空間
平滑化を行う系列データを数学的に理想化したものとして、数列空間$\mathbb R^{\mathbb Z}$を以下で定義しておく。
\mathbb R^{\mathbb Z} = \{ \{x_t\}_{t\in\mathbb Z} \;\vert\; x_t \in \mathbb R\}
ここで、添字の$t=i$の箇所のみが１で他が０であるような$\mathbb R^{\mathbb Z}$の要素を$e_i$とする。
また、シフト作用素$E$を、$Ee_i=e_{i-1}$で定義する。
さらに、$\mathbb R^{\mathbb Z}$の２つの元$x=\{x_t\}$と$y=\{y_t\}$に対して、内積とL2ノルムを以下で定義する。
\begin{align}
\langle x,y \rangle &= \sum_{t\in \mathbb Z} x_ty_t\\
\lVert x \lVert^2 &=\langle x,x\rangle
\end{align}
これらは収束するとは限らないが、収束するときに限って意味があるものとして定義する。
2. 中央部分の理論
スパンが$2m+1$の平滑化作用素$G$を以下で定義する。
G = \sum_{t=-m}^mL_tE^t
実際、$Y=\sum_iy_ie_i$に対して$GY=\sum_i\sum_tL_{t}y_{i+t}e_i$となり、Grevilleの平滑化の中央部分の公式を表す。以下、パラメータ$L_t$の特徴付けのための議論を行う。
2.1. 多項式の再生
ここで、「多項式の再生」を定義しておく。「作用素$G$が多項式$P(t)$を再生する」とは、$p=\sum_tP(t)e_t$としたときに、$Gp=p$が成り立つことをいう。また、作用素$G$が高々$r$次の多項式をすべて再生して、かつ、$r+1$次の多項式で再生しないものが存在するとき、$G$の再生次数は$r$であるという。
今、$G$の再生次数は$r$であることを要求しよう。この条件は、以下の方程式で書ける。
P_{(m,r)} L=e_{(r)} \tag{1}
ただし、$L$は$2m+1$次元列ベクトル,$P_{(m,r)}$は$(2m+1)\times (r+1)$次の行列,$e_{(r)}$は$r+1$次元列ベクトルであり、各成分は以下のとおり。
\begin{align}
L &= [L_{-m},\cdots,L_m]^{\mathrm{T}}\\
\; \\
P_{(m,r)} &= \left(
\begin{array}{ccccccc}
1&1&\cdots&1&\cdots&1&1\\
(-m)&(-m+1)&\cdots&0&\cdots&m-1&m\\
(-m)^2&(-m+1)^2&\cdots&0&\cdots&(m-1)^2&m^2\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\
(-m)^r&(-m+1)^r&\cdots&0&\cdots&(m-1)^r&m^r\\
\end{array}\right)\\
\;\\
e &= [1,0,\cdots,0]^{\mathrm{T}}
\end{align}
2.2. 平滑化指数
前進差分作用素を$\Delta=E-1$で定義する。$y_t$を任意の数列、$z_t$を独立同分布の$z_t\sim N(0,\sigma)$の列として、$Y=\sum_t(y_t+z_t)e_t$で与えられるデータを考える。この$Y$の平滑化後のデータ$GY$の$k$階差分$\Delta^kGY$を平均的に最小化することを考えよう。任意の$s\in\mathbb Z$に対して、$\mathbb E[\langle e_s,\Delta^kGY\rangle^2]$を計算すると、
\mathbb E[\langle e_s,\Delta^kGY\rangle^2] = \sigma^2\Vert D_{(m,k)}L\Vert^2+c_s
となる。ここで、$L=[L_{-m},\cdots,L_m]^{\mathrm{T}}$であり、$c_s$は$L$には無関係の定数部分である。また、$D_{(m,k)}$は$(2m+1+k)\times(2m+1)$の行列で、その成分は以下のとおり。
\begin{align}
D_{(m,k)} &= \left(
\begin{array}{ccccccc}
C_0^{(k)}&C_1^{(k)}&\cdots&C_k^{(k)}&0&\cdots&0\\
0&C_0^{(k)}&\cdots&C_{k-1}^{(k)}&C_{k}^{(k)}&\cdots&0\\
\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots\\
0&\cdots&C_0^{(k)}&C_1^{(k)}&\cdots&C_{k-1}^{(k)}&C_k^{(k)}\\
\end{array}\right)^T\\
\; \\
&\text{ただし、} \;C_j^{(k)} =(-1)^{k-j}\frac{k!}{j!(k-j)!}
\end{align}
ここで、$\Vert D_{(m,k)}L\Vert^2$を平滑化公式$G=\sum_tL_tE^t$の$k$次平滑化指数という。
2.3. 最適化の解
上記(1)式の制約のもと、$\Vert D_{(m,k)}L\Vert^2$を最小化する。このためには、未定乗数$\lambda^{\mathrm{T}}=[\lambda_1\cdots,\lambda_{r+1}]$を導入して、損失関数$\Psi(L,\lambda)$を、以下で定義すればよい。
\begin{align}
\Psi(L,\lambda) = \frac{1}{2}\Vert D_{(m,k)}L\Vert^2 + \lambda^{\mathrm{T}}(P_{(m,r)} L-e_{(r)})
\end{align}
この解は、以下のように求められる。
L = (D_{(m,k)}^TD_{(m,k)})^{-1} P_{(m,r)}^T (P_{(m,r)} (D_{(m,k)}^TD_{(m,k)})^{-1} P_{(m,r)}^T)^{-1}e_{(r)}
ここで、$\mathcal P L = [L_{m},\cdots,L_{-m}]^{\mathrm{T}}$とすると、$\Psi(L,\lambda)=\Psi(\mathcal P L,\lambda)$となるから、$L_t=L_{-t}$が成り立つことがわかる。
3. 両端部分の理論
所与の有限長データ$y_0,\cdots,y_{N-1}$について、その左端の外挿値$y_{-1},y_{-2},\cdots$と右端の外挿値$y_{N},y_{N+1},\cdots$を考える。ここでは左右の外挿値をパラメータとみなそう。そこで、$k=1,2,\cdots$に対して、
\begin{align}
\theta^{-}_k &= y_{-k}\\
\theta^{+}_k &= y_{N-1+k}
\end{align}
とおいて、これらをまとめて$\theta$と書くことにし、外挿後の系列データはこの$\theta$によってパラメータづけられていると考えて、以下のように表示する。
Y_{\theta} = \sum_i y_ie_i
$\theta$を決めるための損失関数の候補として考えられるのは、内積が２つのベクトルの類似度を表すことから、原系列$Y_\theta$と平滑化後の系列$GY$の内積のマイナス倍
\Phi_0(\theta)=-\langle Y_\theta,GY_\theta\rangle
が考えれられる。しかし、$\Phi_0(\theta)$は下に有界とは限らないから、制約条件$\Vert{Y_\theta}\Vert=c$を追加する。ラグランジュ未定乗数$\mu$も導入して、
\begin{align}
\Phi_\mu(\theta) &= -\langle Y_\theta,GY_\theta\rangle + \mu(\Vert{Y_\theta}\Vert^2-c^2)\\
&=\langle Y_\theta,(\mu-G)Y_\theta \rangle + \mathrm{CONST}
\end{align}
ここで、$\mu$と$c$はお互いに関係しているので、制約条件$\Vert{Y_\theta^2}\Vert=c$のことは忘れて、$\mu$をモデルパラメータとして解釈しなおしてもよい。このとき、$\mu$はL2正則化のパラメータと見做される。
Grevilleは$\mu=1$を採用した。以下、$\mu=1$とする。
3.1. Cramer-Wold分解
ここで$\Phi_1(\theta)$を評価するため、$1-G$を変形する。$1-G$は実係数ローラン多項式$f(z)$を用いて、以下のようにかける。
1-G = f(E)
ここで、ラプラス作用素$\delta^2$を$\delta^2 = E+E^{-1}-2=-(E-1)(E^{-1}-1)$で定める。すると、ある多項式$g$があって$f(E)=g(\delta^2)$とできる。このことは次のようにわかる。いま、$L_t=L_{-t}$なので、$f(E)$は$E^t+E^{-t}$という項の線型結合となる。これらが$\delta^2$の多項式で書けることは$t$に関する数学的帰納法で示せる。したがって、$f(E)$は$\delta^2$の多項式で書ける。
さて、$\delta^{2s}$は高々$2s-1$次の多項式$P(t)$がつくるデータ$p=\sum_tP(t)e_t$に対して、$\delta^{2s}p=0$となる。
一方で、$G$の再生次数は$r$であるから、高々$r$次の多項式$P(t)$がつくるデータ$p=\sum_tP(t)e_t$に対して、$(1-G)p=0$となる。
このことから、$2s>r$が成り立つ最小の$s$をとると、適当な係数$c_k$を用いて、$g(E)=\sum_{k\geq s}c_k \delta^{2k}$という形に書ける。なお、$r$が奇数の時は$r=2s-1$,rが偶数のときは$r=2s-2$である。
以上より、$1-G$は次のように書けることがわかる。
1-G = \delta^{2s}\hat f(E)
いま、$\hat f$は$\hat f(z)=\hat f(z^{-1})$が成り立つ実係数ローラン多項式であり、$z^{-(m-2s)}$から$z^{m-2s}$までの項を含む。ここで、$\hat f$は単位円上にゼロ点を持たないと仮定する（このことは数値的な確認が必要になる。）。このとき、$z^{m-2s}\hat f(z)$の因数分解を以下のように表示することができる。
z^{m-2s}\hat f(z) = c \prod_{i=1}^{m-2s} (z-\alpha_i)(1-\alpha_iz)
ただし、$\vert\alpha_i\vert>1$とし、$c$は適当な実数である。ここで多項式$R$を、
R(z) = c^{1/2}(z-1)^s(z+1)^s\prod_{i=1}^{m-2s} (z-\alpha_i)
で定めると、$1-G=(-1)^sR(E)R(E^{-1})$となる。これを、Cramer-Wold分解という。
3.2. 解の表示
上述のCramer-Wold分解を用いると、
\Phi_{1}(\theta) = \Vert {R(E)Y_\theta }\Vert^2
となる。
ここで、$Q^{(-)}$を$e_{-1},e_{-2},\cdots$が張る空間への射影,$Q^{(+)}$を$e^{N},e^{N+1},\cdots$が張る空間への射影とする。つまり、$e_t$への作用を以下で定義する。
Q^{(-)}e_t = \left\{
\begin{array}{ll}
e_t & (t \lt 0) \\
0 & (t \geq 0)
\end{array}
\right.
Q^{(+)}e_t = \left\{
\begin{array}{ll}
e_t & (t \geq N) \\
0 & (t \lt N)
\end{array}
\right.
このとき、
\Phi_{1}(\theta) = \Vert {Q^{(-)}R(E)Y_\theta }\Vert^2 + \Vert {(1-Q^{(-)})R(E)Y_\theta }\Vert^2
となり、第２項は$\theta^{(-)}$には依存しない。また、
\Phi_{1}(\theta) = \Vert {Q^{(+)}R(E^{-1})Y_\theta }^2\Vert + \Vert {(1-Q^{(+)})R(E^{-1})Y_\theta }\Vert^2
とすれば、第２項は$\theta^{(+)}$には依存しない。よって、
\Phi_{1}^{(\mp)}(\theta) = \Vert {Q^{(\mp)}R(E^{\pm 1})Y_\theta }\Vert^2
とおけば（復号同順）、これらをそれぞれ最小化すれば$\Phi_{1}(\theta) $の最小化が得られる。実際、以下に述べるとおり、
Q^{(\mp)}R(\pm E)Y_\theta = 0
には、解があることが確かめられるので、このようにして得られた$\theta$こそが求める解である。
それでは、解を具体的に表示しよう。まず、係数$\rho_k$を導入して、$R(z)=\sum_{k=0}^{m}\rho_kz^k$とする。すると、$\theta^{(-)}$については
\begin{align}
0 &= \rho_0 \theta^{(-)}_{1} + \rho_1y_0 + \rho_2 y_1 + \cdots + \rho_my_{m-1} \\
0 &= \rho_0 \theta^{(-)}_{2} + \rho_1\theta^{(-)}_{1} + \rho_2 y_0 + \cdots + \rho_my_{m-2} \\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\ \vdots
\end{align}
であり、$\theta^{(+)}$については
\begin{align}
0 &= \rho_0 \theta^{(+)}_{1} + \rho_1y_{N-1} + \rho_2 y_{N-2} + \cdots + \rho_my_{N-m} \\
0 &= \rho_0 \theta^{(+)}_{2} + \rho_1\theta^{(+)}_{1} + \rho_2 y_{N-1} + \cdots + \rho_my_{N-m+1} \\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\ \vdots
\end{align}
となる。これらの方程式はそれぞれ上から順に解いていけばよい。これが、Grevilleの平滑化の両端部分の公式である。Grevilleの数表の$a_t$との対応は、$a_t=-\rho_t/\rho_0$である。
最後に、Cramer-Wold分解の非一意性について簡単に言及しておく。上記の議論では、$R(z)$の１以外の根がすべて単位円の外側にあるように選んだ。一方で、１以外の根がすべて単位円の内側にあるように選ぶことも可能である。このことは、$\Phi_1(\theta)$の極値がひとつではないことを意味する。
ここで、$\theta^{(\pm)}_N$を決めるための漸化式の特性方程式を考えればすぐにわかることだが、単位円の外側に選んだ場合、$\theta^{(\pm)}_N$の$N\to\infty$における振る舞いが、高々$\mathcal{O}(N^{s-1})$程度となる。逆に、単位円の内側に選んだ場合は、指数関数的な発散が生じうる。理論的には、この理由で単位円の外側が選ばれている。
実際に、単位の内側を選んでみると、$a_1<0$となったり、係数の絶対値が大きくなったりするなど、外挿の係数としては採用しにくい解が出てくることが確かめられる。
4. 係数を求めるコード(Python-Numpy)
import numpy as np
import scipy.special
def calc_Greville_center(m,k,r):
#calc D
C = np.zeros(2*m+1+k)
for j in range(k+1):
C[j] = ( (-1)**(k-j) ) * scipy.special.binom(k,j)
D = np.array([np.roll(C,i) for i in range(2*m+1)]).T
#calc P
P = np.zeros([1+r,2*m+1])
for i in range(0,r +1):
for j in range(-m,m +1):
P[i][j+m] = (j**i)
#calc e
e = np.zeros([r+1,1])
e[0][0]=1
#calc L
INV = np.linalg.inv
L = INV(D.T@D) @ (P.T) @ INV(P@INV(D.T@D)@(P.T)) @ e
L = L.T[0]
return L
def calc_Greville(m,k,r):
if r%2==0 :
s = int((r+1)/2 ) # r=2s-1
else:
s = int((r+2)/2 ) # r=2s-2
L = calc_Greville_center(m,k,r)
# G(z) = Lm*z^-m + ... + L0 + ... + Lm*z^m
# coef of H(z) = z^m*(1-G(z))
H = -L
H[m] +=1
# calc roots outside the unit circle in C
roots = np.roots(H)
idx = np.argsort(-np.abs(roots))[0:m-s]
roots = roots[idx]
# calc poly R such that 1-G(z) = R(z)R(z^-1)
R = ((np.poly1d([1,-1]))**s) * np.poly1d( np.poly(roots) )
a = -np.flip(R.coef/R.coef[-1])[1:]
return L,a
for m in range(2,7):
L,a = calc_Greville(m,k=3,r=3)
print(m,"L",np.round(L[m:],6))
print(m,"a",np.round(a,6))
※このコードでは、Grevilleの論文(1979)の数表の小数第６位を再現できない。Grevilleは、端数処理の段階で何らか最適化をした可能性がある。
5. 参考文献
[1979] T.N.E. Greville. Moving-weighted-average smoothing extended to the extremities of tha data, MRC Technical Summary Report # 2025. Mathematics Research Center, University of Wisconsin.
[1980] 藤田輝昭 「Grevilleの補正公式について-両端の処理に関する最近の話題-」 社団法人アクチュアリー会会報 第33号 第2分冊
筆者はあまり目を通していないが、Grevilleの一連の論文を挙げておく。
[1947] T.N.E. Greville. Acturial note : Adjusted average graduation formulas of maximum smoothness, R.A.I.A. 36, 249-264.
[1948] T.N.E. Greville. Acturial note : Tables of coefficients in adjusted average graduation formulas of maximum smoothness, R.A.I.A. 37, 11-30.
[1974] T.N.E. Greville. Parts 5 Study Notes-Graduation (53. 1. 73). Education and Examination Committee of the Society of Actuaries.