mambaの理論を理解する①：HiPPOフレームワークとLSSL

はじめに

この前のブログでは、mambaの論文を翻訳した。

https://izmyon.hatenablog.com/entry/2023/12/11/155551

本シリーズでは、mambaの理論的背景を理解するために、それらの先行研究を順々にまとめて解説していく。重要な先行研究としては、LMU, HiPPO, LSSL, S4, H3などがあり、主に以下の流れで解説する。

1. HiPPOフレームワークとLSSL
2. S4のアルゴリズム
3. H3とHyena Hierarchy
4. mamba

本記事は、第一段として、HiPPOフレームワークとLSSLについてまとめる。

元論文を以下に記す。

HiPPO

https://arxiv.org/abs/2008.07669

LSSL

https://arxiv.org/abs/2110.13985

この記事の続きはこちら

https://zenn.dev/izmyon/articles/c56a2fd6670546

読者の方へ

• 補足や訂正などがあれば、コメントにて優しく丁寧にご教示いただけると喜びます。
• 書くのけっこう大変だったのでバッチを送っていただけると嬉しいです。たくさんもらえると僕のやる気が上がって投稿頻度が上がるかもしれません。（逆にもらえないと下がるかもしれない。）

状態空間モデル

時刻$tt$における1次元の入力信号$u(t)u(t)$を、$NN$次元の潜在空間$x(t)x(t)$に投影し、出力信号$y(t)y(t)$を計算する以下の状態空間モデルで表される線形システムを考える。LSSL、S4、H3、mambaなどの一連の状態空間モデルは、下記の方程式を離散化し、長さ$LL$の系列を同じく長さ$LL$の系列に変換する$RL→RL\mathbb{R}^L \rightarrow \mathbb{R}^L$のsequence-to-sequenceモデルを用いる。

$$x˙(t)=Ax(t)+Bu(t)y(t)=Cx(t)+Du(t) \begin{align} \dot{x}(t) = A x(t) + B u(t) \quad \\ y(t) = C x(t) + D u(t) \quad \end{align}$$

この方程式で表される状態空間モデルは、制御工学でよく用いられるものである。
制御工学でどのように用いられるかについては以下のブログ（文献[8]）がわかりやすい。
https://blog.control-theory.com/entry/state-feedback

状態空間モデルは、Hidden Markov Models (HMM）のように、制御工学や機械学習などの幅広い分野でみられる潜在空間を利用するモデルであり、行列$A∈RN×NA \in \mathbb{R}^{N \times N}$、$B∈RN×1B \in \mathbb{R}^{N \times 1}$、 $C∈R1×NC \in \mathbb{R}^{1 \times N}$、$D∈R1×1D \in \mathbb{R}^{1 \times 1}$は、勾配降下によって求められる学習パラメータである。

これは連続時間のシステムであるが、コンピュータ上では離散時間信号を扱うため、式(1)、式(2)の離散化を行う。ここで、式(2)の離散化は簡単だが問題なのは式(1)の常微分方程式（ODE）を離散化することである。

ゼロ次ホールド（ZOH）と同様に、時刻$tt$における入力信号$u(t)u(t)$が次のタイムステップ$Δt\Delta t$まで一定である、つまり$u(t+Δt)=u(t)u(t+\Delta t) = u(t)$であると仮定する。(1)式の右辺を$tt$と$x(t)x(t)$の関数として、$x˙(t)=f(t,x(t))\dot{x}(t) = f(t, x(t))$と表す、離散化された時刻$tit_i$において、状態$x(t0),x(t1),…x(t_0), x(t_1), \ldots$は、ピカールの逐次近似法$x(ti+1)=x(ti)+∫titi+1f(s,x(s))dsx(t_{i+1}) = x(t_i) + \int^{t_{i+1}}_{t_i} f(s, x(s))ds$によって求められる。この右側積分を推定する方法には様々な方法があり、ここでは、Generalized bilinear transform (GBT)を用いる。この方法では、定数$α\alpha$を用いて、以下のように右側積分を近似する。

$$∫tt+Δf(s,x(s))ds≃αΔtf(t+Δ,x(t+Δt))+(1−α)Δtf(t,x(t))=αΔt(Ax(t+Δ)+Bu(t+Δ))+(1−α)Δt(Ax(t)+Bu(t))=αΔt(Ax(t+Δ)+Bu(t))+(1−α)Δt(Ax(t)+Bu(t))(∵u(t+Δ)=u(t))=αΔtAx(t+Δ)+(1−α)ΔtAx(t)+BΔtu(t) \begin{align*} \begin{split} &\int^{t+\Delta}_{t} f(s, x(s))ds \\ &\simeq \alpha \Delta t f(t+\Delta, x(t+\Delta t)) + (1-\alpha) \Delta t f(t, x(t)) \\ &= \alpha \Delta t (A x(t+\Delta) + B u(t+\Delta)) + (1-\alpha) \Delta t (A x(t) + B u(t)) \\ &= \alpha \Delta t (A x(t+\Delta) + B u(t)) + (1-\alpha) \Delta t (A x(t) + B u(t)) \quad (\because u(t+\Delta) = u(t)) \\ &= \alpha \Delta t A x(t+ \Delta) + (1-\alpha) \Delta t A x(t) + B \Delta t u(t) \end{split} \end{align*}$$

この結果を用いると、$x(t+Δt)=x(t)+∫tt+Δtf(s,x(s))dsx(t+\Delta t) = x(t) + \int^{t+\Delta t}_{t} f(s, x(s))ds$より、

$$x(t+Δt)=x(t)+αΔtAx(t+Δt)+(1−α)ΔtAx(t)+BΔtu(t)(I−αΔtA)x(t+Δt)=(I+(1−α)ΔtA)x(t)+BΔtu(t)∴x(t+Δt)=(I−αΔtA)−1(I+(1−α)ΔtA)x(t)+(I−αΔtA)−1BΔtu(t) \begin{align*} x(t+ \Delta t) &= x(t) + \alpha \Delta t A x(t+ \Delta t) + (1-\alpha) \Delta t A x(t) + B \Delta t u(t) \notag \\ (I- \alpha \Delta t A) x(t+\Delta t) &= (I + (1-\alpha) \Delta t A) x(t) + B \Delta t u(t) \notag \\ \therefore x(t+\Delta t) &= \left(I - \alpha \Delta t A \right)^{-1} \left(I + (1- \alpha) \Delta t A\right) x(t) + \left( I - \alpha \Delta t A \right)^{-1} B \Delta t u(t) \end{align*}$$

ここでは、$α=12\alpha = \frac{1}{2}$として、下式を用いる。

$$x(t+Δt)=(I−12ΔtA)−1(I+12ΔtA)x(t)+(I−12ΔtA)−1BΔtu(t) \begin{align} \begin{split} x(t+\Delta t) = \left(I - \frac{1}{2} \Delta t A \right)^{-1} \left(I + \frac{1}{2} \Delta t A\right) x(t) + \left( I - \frac{1}{2} \Delta t A \right)^{-1} B \Delta t u(t) \end{split} \end{align}$$

さらに、以下のように行列$A,BA, B$を行列$Aˉ,Bˉ\bar{A}, \bar{B}$に変換する。

$$Aˉ=(I−12ΔtA)−1(I+12ΔtA)Bˉ=(I−12ΔtA)−1BΔt \begin{align} \bar{A} = \left(I - \frac{1}{2} \Delta t A \right)^{-1} \left(I + \frac{1}{2} \Delta t A\right) \\ \bar{B} = \left( I - \frac{1}{2} \Delta t A \right)^{-1} B \Delta t \end{align}$$

すると、以下のように式(1), (2)の連続時間システムを離散化することが出来る。

$$xt=Aˉxt−1+Bˉutyt=Cxt+Dut \begin{align} x_t &= \bar{A} x_{t-1} + \bar{B} u_t \quad \\ y_t &= C x_t + D u_t \quad \end{align}$$

GBTは、$α=0\alpha = 0$とすると、オイラー法$x(t+Δt)=x(t)+Δtf(t,x(t))x(t+\Delta t)= x(t) + \Delta t f(t, x(t))$と一致し、$α=1\alpha = 1$とすると、後退オイラー法$x(t)=x(t+Δt)−Δtf(t+Δt,x(t+Δt))x(t) = x(t + \Delta t) - \Delta t f(t+\Delta t, x(t+\Delta t))$と一致する。実は、GBTはこれらの手法を定数$α∈[0,1]\alpha \in [0, 1]$を用いて一般化した手法である。

積分$∫tt+Δf(s,x(s))ds\int^{t+\Delta}_{t} f(s, x(s))ds$の近似方法に着目するならば、$α=0\alpha=0$のオイラー法の場合は、下図のように、高さ$f(t,x(t))f(t, x(t))$の長方形で近似し、$α=1\alpha=1$の後退オイラー法の場合は、高さ$f(t+Δt,x(t+Δt))f(t+\Delta t, x(t+\Delta t))$の長方形で近似する。一方で、$α=12\alpha = \frac{1}{2}$の場合には、左辺が$f(t,x(t))f(t, x(t))$、右辺が$f(t+Δt,x(t+Δt))f(t+\Delta t, x(t+\Delta t))$の台形で近似するということである。

HiPPO (High-Order Polynomial Projection Operator; 高次多項式投影演算子)

ここまで線形システム(1)、(2)を離散化し、新たなシステムである式(6),(7)を得た。次に、本節で説明するHiPPOと呼ばれるフレームワークを適用し、以下のように係数行列$AA$を初期化する。

$$Ank={(2n+1)1/2(2k+1)1/2if n>kn+1if n=k0if n<k A_{nk} = \left\{ \begin{array}{ll} (2n+1)^{1/2} (2k+1)^{1/2} & \text{if} \ n>k \\ n+1 & \text{if} \ n=k \\ 0 & \text{if} \ n<k \\ \end{array} \right.$$

なぜこのような初期化を行うのかをこれから説明化するが、長いのでLSSLの全体像を先につかみたい場合は、次の節「畳み込みによる高速化」に進んでもらって構わない。

HiPPOの定義

時間変化する$(−∞,t](-\infty, t]$でサポートされる測度族$μ(t)\mu^{(t)}$、$NN$次元の（直交）多項式が張る部分空間$G\mathcal{G}$、そして連続入力信号$u:R≥0→Ru: \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}$が与えられた時、$uu$を最適化された投影係数$x:R≥0→RNx: \mathbb{R}_{\geq 0} \rightarrow \mathbb{R}^N$に移す演算$hippo\text{hippo}$を定義する。この演算は、以下のように投影演算子$projt\text{proj}_t$と、係数抽出演算子$coeft\text{coef}_t$をすべての時間ステップ$tt$で求め、それらの合成$coeft∘projt\text{coef}_t \circ \text{proj}_t$である。（つまり、$(hippo(u))(t)=coeft(projt(u))(\text{hippo}(u))(t)=\text{coef}_t (\text{proj}_t (u))$である。）

1. $projt\text{proj}_t$は、時刻$tt$までの信号$uu$、つまり$u<t:={u(y)}y≤tu_{<t}:=\{u(y)\}_{y \leq t}$を、推論誤差$∥u≤t−g(t)∥L2(μ(t))\| u_{\leq t} - g^{(t)} \|_{L_2 (\mu^{(t)})}$を最小にする多項式$g(t)∈Gg^{(t)} \in \mathcal{G}$に射影する。

2. $coeft:G→RN\text{coef}_t: \mathcal{G} \rightarrow \mathbb{R}^N$は、多項式$g(t)g^{(t)}$を、測度$μ(t)\mu^{(t)}$に関して定義された直交多項式の基底関数の係数$x(t)∈RNx(t) \in \mathbb{R}^N$に射影する。

HiPPOを導くためのHiPPOフレームワーク

${Pn(t)}n∈N\{P_n^{(t)}\}_{n \in \mathbb{N}}$は、時間変化する測度$μ(t)\mu^{(t)}$に関する直交多項式列であり、$pn(t)p_n^{(t)}$を、$Pn(t)P_n^{(t)}$を正規した多項式、すなわち$pn(t)=Pn(t)/∥Pn(t)∥μ2p_n^{(t)} = P_n^{(t)}/ \| P_n^{(t)} \|_{\mu}^2$であるとする。正規直交基底の定義より、以下が成り立つ。

$$∫−∞∞pn(t)(y)pm(t)(y)ω(t)(y)dy=δm,n \begin{align} \int^{\infty}_{-\infty} p_n^{(t)}(y) p_m^{(t)}(y) \omega^{(t)}(y) dy = \delta_{m,n} \end{align}$$

ここで、$ω(t)(y)\omega^{(t)}(y)$は重みであり、$dμ(t)=ω(t)(y)dy,∫−∞∞dμ(t)=1d \mu^{(t)} = \omega^{(t)}(y) dy, \int^{\infty}_{-\infty} d \mu^{(t)} = 1$である。

ここでは直交多項式を基底関数に選択し、より簡易な設定で説明を試みている。しかし、HiPPOフレームワークはより一般なもので、必ずしも基底として直交基底関数を選択する必要はなく、さらには多項式である必要もないものとして定式化されている。詳しい証明や定式化については、以下の参考文献[5]のAppendix Cを参照のこと。

https://arxiv.org/abs/2206.12037

入力信号の履歴$u≤t:={u(y)}y≤tu_{\leq t}:=\{u(y)\}_{y \leq t}$が、ステップサイズ$Δt\Delta t$で連続関数$u(y)u(y)$からサンプリングされたと考える。N個の基底${pn(t)}n=0,1,…,N−1\{p_n^{(t)}\}_{n=0,1,\ldots, N-1}$を用いて、$u≤tu_{\leq t}$を通る全時刻にわたる連続関数$u(y)u(y)$を多項式$g(t)g^{(t)}$で近似する。すなわち、$u≤tu_{\leq t}$を${pn(t)}n=0,1,…,N−1\{p_n^{(t)}\}_{n=0,1,\ldots, N-1}$が張る空間へ射影し、多項式$g(t)g^{(t)}$が得られるとする。

$$u≤t≃g(t)=∑n=0N−1xn(t)pn(t) \begin{align} u_{\leq t} \simeq g^{(t)} = \sum_{n=0}^{N-1} x_n (t) p_n^{(t)} \end{align}$$

このように近似すると、時刻$tt$までの入力信号の履歴を、$NN$個の係数${xn(t)}n=0,1,…,N−1\{x_n^{(t)}\}_{n=0,1,\ldots, N-1}$で表現することが出来るようになる。係数$xn(t)x_n (t)$は、時刻$tt$までの入力信号の履歴$u≤tu_{\leq t}$を用いて、フーリエ係数を求める要領で以下のように求めることが出来る。

$$xn(t)=∫−∞∞g(t)pn(t)ω(t)(y)dy≃∫−∞∞u≤tpn(t)ω(t)(y)dy \begin{align} x_n (t) = \int^{\infty}_{-\infty} g^{(t)} p_n^{(t)} \omega^{(t)}(y) dy &\simeq \int^{\infty}_{-\infty} u_{\leq t} p_n^{(t)} \omega^{(t)}(y) dy \end{align}$$

なぜこのように求められるのかは、右辺から逆に求めると理解できる。

$$∫−∞∞g(t)pn(t)ω(t)(y)dy=∫−∞∞∑m=0N−1(xm(t)pm(t))pn(t)ω(t)(y)dy=∑m=0N−1xm(t)∫−∞∞pm(t)pn(t)ω(t)(Y)dy=∑m=0N−1xm(t)δm,n (∵式(8)を用いた。)=xn(t) \begin{align*} \int^{\infty}_{-\infty} g^{(t)} p_n^{(t)} \omega^{(t)}(y) dy &= \int^{\infty}_{-\infty} \sum_{m=0}^{N-1} (x_m (t) p_m^{(t)}) p_n^{(t)} \omega^{(t)}(y) dy \\ &= \sum_{m=0}^{N-1} x_m (t) \int^{\infty}_{-\infty} p_m^{(t)} p_n^{(t)} \omega^{(t)}(Y) dy \\ &= \sum_{m=0}^{N-1} x_m (t) \delta_{m,n} \quad (\because 式(8)を用いた。) \\ &= x_n (t) \end{align*}$$

実は、式(10)に従って、以下の入力信号$u≤tu_{\leq t}$の履歴から$NN$個の係数を計算し、出力することが$hippo\text{hippo}$演算に他ならない。

$$hippo(u≤t):=(x0,x1,…,xN−1)where u≤t≃g(t)=∑n=0N−1xn(t)pn(t) \begin{align*} \text{hippo}(u_{\leq t}) := (x_0, x_1, \ldots, x_{N-1}) \\ \text{where} \ u_{\leq t} \simeq g^{(t)} = \sum_{n=0}^{N-1} x_n (t) p_n^{(t)} \end{align*}$$

しかしながら、これらの係数を各時刻で計算するのは容易ではない。そこでHiPPOフレームワークでは、式(1)を離散化して式(6)の更新式を導いたように、$x(t)x(t)$が満たすODEを求めた後、それを離散化して更新式を求め、各時刻で$x(t)x(t)$を更新していくことを考える。

ということでHiPPOでは式(10)を用いて$xn(t)x_n(t)$を求めることはしないわけだが、式(10)の右辺を計算する方法を考えてみる。

まず、$u≤tu_{\leq t}$の$kk$番目の要素が$u[k],k=0,…,L−1u[k], k = 0, \ldots, L-1$で表されるとし、$u≤tu_{\leq t}$を以下のように定義し直す。

$$u≤t(y)={u[k]ify=kΔt,k=0,…,L−10otherwise \begin{align*} u_{\leq t}(y) = \left\{ \begin{array}{ll} u[k] & \text{if} \quad y = k\Delta t, \quad k = 0, \ldots, L-1\\ 0 & otherwise \\ \end{array} \right. \end{align*}$$

これを用いて、

$$∫−∞∞u≤tpn(t,y)ω(t,y)dy=∫−∞∞(∑k=0L−1u[k]δ(y−kΔt))pn(t,y)ω(t,y)dy=∑k=0L−1u[k]∫−∞∞δ(y−kΔt)pn(t,y)ω(t,y)dy=∑k=0L−1u[k]pn(t,kΔt)ω(t,kΔt) \begin{align*} \int^{\infty}_{-\infty} u_{\leq t} p_n (t, y) \omega (t,y) dy &= \int^{\infty}_{-\infty} ( \sum^{L-1}_{k=0} u[k]\delta(y-k \Delta t) ) p_n (t, y) \omega (t,y) dy \\ &= \sum^{L-1}_{k=0} u[k] \int^{\infty}_{-\infty} \delta(y-k \Delta t) p_n (t, y) \omega (t,y) dy \\ &= \sum^{L-1}_{k=0} u[k] p_n (t, k \Delta t) \omega (t, k \Delta t) \end{align*}$$

従って、以下の式で$xn(t)x_n(t)$を計算できるのではないかと個人的には考えている。

$$xn(t)≃∑k=0L−1u[k]pn(t,kΔt)ω(t,kΔt) \begin{align*} x_n (t) \simeq \sum^{L-1}_{k=0} u[k] p_n (t, k \Delta t) \omega (t, k \Delta t) \end{align*}$$

ちなみにここで、$nn$を$ω\omega$と入れ替えて、基底を$pω(t)(y)=e−jωyp_{\omega}^{(t)} (y) = e^{-j \omega y}$、測度を$ω(t,y)=1\omega (t, y) = 1$とし、求める係数を$x(t)(ω)x^{(t)} (\omega)$としてみると、

$$x(t)(ω)=∑k=0L−1u[k]e−jωkΔt \begin{align*} x^{(t)} (\omega) = \sum^{L-1}_{k=0} u[k] e^{-j \omega k \Delta t} \end{align*}$$

となり、これは離散フーリエ変換（DFT）と考えることが出来る。

係数が満たすODE

以上より、下式が成り立つことが分かった。

$$xn(t)=∫−∞∞u(y)pn(t,y)ω(t,y)dy \begin{align} x_n (t) = \int^{\infty}_{-\infty} u(y) p_n (t, y) \omega(t, y) dy \end{align}$$

両辺を微分して、係数$xn(t)x^{(t)}_n$が満たすODEを導く。ここで、微分と積分の可換性を認めた。

$$x˙n(t)=∂∂t∫−∞∞u(y)pn(t,y)ω(t,y)dy=∫−∞∞∂∂t(u(y)pn(t,y)ω(t,y))dy=∫−∞∞u(y)(∂∂tpn(t,y))ω(t,y)dy +∫−∞∞u(y)pn(t,y)(∂∂tω(t,y))dy \begin{align*} \dot{x}_n (t) &= \frac{\partial}{\partial t} \int^{\infty}_{-\infty} u(y) p_n (t, y) \omega(t, y) dy \\ &= \int^{\infty}_{-\infty} \frac{\partial}{\partial t} (u(y) p_n (t, y) \omega(t, y) ) dy \\ &= \int^{\infty}_{-\infty} u(y) (\frac{\partial}{\partial t} p_n (t, y)) \omega (t, y) dy + \int^{\infty}_{-\infty} u(y) p_n (t, y) (\frac{\partial}{\partial t} \omega (t, y) ) dy \end{align*}$$

以上より、以下の式(11)が成り立つ。

$$x˙n(t)=∫−∞∞u(y)(∂∂tpn(t,y))ω(t,y)dy +∫−∞∞u(y)pn(t,y)(∂∂tω(t,y))dy \begin{align} \dot{x}_n (t) = \int^{\infty}_{-\infty} u(y) (\frac{\partial}{\partial t} p_n (t, y)) \omega (t, y) dy + \int^{\infty}_{-\infty} u(y) p_n (t, y) (\frac{\partial}{\partial t} \omega (t, y) ) dy \end{align}$$

HiPPOフレームワークでは、具体的に用いる基底${Pn(t)}n∈N\{P_n^{(t)}\}_{n \in \mathbb{N}}$及び、測度$μ(t)\mu^{(t)}$を選択し、式(12)に代入することで、$x˙(t)\dot{x}(t)$とその時点での係数$x(t)x(t)$、入力信号$u(t)u(t)$との関係を導き、式(1)に当てはめて係数行列$AA$と$BB$を設定する。

HiPPO-LegS

HiPPOの論文では、様々な基底や測度を用いた場合について検証されており、LSSLではその中で最も良い結果を示したScaled Legendre Measure (LegS)という設定を用いて得られる行列$AA$に初期化する。この設定を用いたHiPPO-LegSは、基底としてルシャンドル多項式を用い、すべての履歴に対して一様に重みをつける$ωt=1t1[0,t]\omega_t = \frac{1}{t} \mathbb{1}_{[0, t]}$を用いるものである。

ここではまず、ルシャンドル多項式の性質や、基底として用いた場合の$ωt\omega_t$に関する直交性などを確認し、正規直交基底を求め、式(12)に代入してODEを求める。

正規直交規定および測度

ルシャンドル多項式は、以下の$nn$次多項式である。

$$Pn(x):=12nn!dndxn[(x2−1)n] \begin{align} P_n (x) := \frac{1}{2^n n!} \frac{d^n}{{dx}^n} [ (x^2 - 1 )^n ] \end{align}$$

測度$ωleg=1[−1,1]\omega^{\text{leg}} = \bold{1}_{[-1, 1]}$に関して、以下の直行性が成り立つ。

$$2n+12∫−11Pn(x)Pm(x)dx=δnm \begin{align} \frac{2n+1}{2} \int^1_{-1} P_n(x) P_m(x) dx = \delta_{nm} \end{align}$$

さらに、以下の性質を持つ。

$$Pn(1)=1Pn(−1)=(−1)n(2n+1)Pn=Pn+1′−Pn−1′Pn+1′=(n+1)Pn+xPn′ \begin{align} P_n (1) &= 1 \\ P_n (-1) &= (-1)^n \\ (2n+1) P_n &= P'_{n+1} - P'_{n-1} \\ P'_{n+1} &= (n+1) P_n + x P'_n \end{align}$$

以上のルシャンドル多項式は区間$x∈[−1,1]x \in [-1, 1]$で成り立つ性質であったが、区間[0, t]についても成り立つようにスケーリングするため、$y=t2(x+1)y = \frac{t}{2}(x+1)$、つまり$x=2ty−1x = \frac{2}{t} y - 1$と変数変換する。すると、$dx=2tdydx = \frac{2}{t} dy$であり、式(14)の左辺を変形すると、以下のようになる。

$$2n+12∫−11Pn(x)Pm(x)dx=2n+12∫0tPn(2ty−1)Pm(2ty−1)2tdy=(2n+1)∫0tPn(2ty−1)Pm(2ty−1)ωleg(2ty−1)1tdy \begin{align*} & \frac{2n+1}{2} \int^1_{-1} P_n(x) P_m(x) dx \\ =& \frac{2n+1}{2} \int^t_{0} P_n( \frac{2}{t} y - 1 ) P_m( \frac{2}{t} y - 1 ) \frac{2}{t}dy \\ =& (2n+1) \int^t_{0} P_n( \frac{2}{t} y - 1 ) P_m( \frac{2}{t} y - 1 ) \omega^{\text{leg}} (\frac{2}{t} y - 1 ) \frac{1}{t} dy \end{align*}$$

従って、以下が成り立つ。

$$(2n+1)∫0tPn(2ty−1)Pm(2ty−1)ωleg(2ty−1)1tdy=δnm \begin{align} (2n+1) \int^t_{0} P_n( \frac{2}{t} y - 1 ) P_m( \frac{2}{t} y - 1 ) \omega^{\text{leg}} (\frac{2}{t} y - 1 ) \frac{1}{t} dy = \delta_{nm} \end{align}$$

この式より、新たな測度$ωt=1t1[0,t]\omega_t = \frac{1}{t} \mathbb{1}_{[0, t]}$を用い、さらに以下の正規直交規定を用いることにする。

$$pn(t,y)=(2n+1)1/2Pn(2ty−1) \begin{align} p_n (t, y) = (2n+1)^{1/2} P_n ( \frac{2}{t} y - 1 ) \end{align}$$

すると、以下が成り立つ。

$$∫0tpn(t,y)pm(t,y)ωtdy=δnm \begin{align} \int^t_{0} p_n(t, y) p_m(t, y) \omega_t dy = \delta_{nm} \end{align}$$

ODEの導出

HiPPO-LegSが用いるODEを導出する前に、導出で用いるルシャンドル基底の微分が満たす性質について述べる。まず、式(17)から、以下が導かれる。

$$Pn+1′=(2n+1)Pn+(2n−3)Pn−2+…, \begin{align} P'_{n+1} = (2n+1) P_n + (2n-3)P_{n-2} + \ldots, \end{align}$$

次数を一つずらすと、

$$Pn′=(2n−1)Pn−1+(2n−5)Pn−3+…, \begin{align} P'_{n} = (2n-1) P_{n-1} + (2n-5)P_{n-3} + \ldots, \end{align}$$

これらの性質を用いて$(x+1)Pn′(x)(x+1)P'_n (x)$を計算すると、以下が成り立つ。

$$(x+1)Pn′=xPn′+Pn′=Pn+1′+Pn′−(n+1)Pn(∵式(18)より、xPn′=Pn′−(n+1)Pn)=nPn+(2n−1)Pn−1+(2n−3)Pn−2+…(∵式(22), (23)を代入) \begin{align} (x+1) P'_{n} =& x P'_{n} + P'_{n} = P'_{n+1} + P'_n - (n+1) P_n \notag \\ &(\because \text{式(18)より、} xP'_n = P'_n - (n+1) P_n) \notag \\ =& nP_n + (2n-1) P_{n-1} + (2n-3)P_{n-2} + \ldots \\ &(\because \text{式(22), (23)を代入}) \notag \\ \end{align}$$

まず、$ωt\omega_t$と、$pn(t,y)p_n(t, y)$の微分をそれぞれ求める。

$$∂∂tω(t)=−t−21[0,1]+t−1δt=t−1(−ω(t)+δt)∂∂tpn(t,y)=−(2n+1)122yt−2Pn′(2ty−1)=−(2n+1)12t−1((2ty−1)+1)Pn′(2ty−1)=−(2n+1)12t−1(x+1)Pn′(x) (∵x=2ty−1)=−(2n+1)12t−1[nPn(x)+(2n−1)Pn−1(x)+(2n−3)Pn−2(x)+…] (∵式(24)を用いた。)=−(2n+1)12t−1[n(2n+1)−12pn(t,y)+(2n−1)12pn−1(t,y)+(2n−3)12pn−2(t,y)+…](∵式(20)より、Pn(x)=(2n+1)−12pn(t,y)) \begin{align} \frac{\partial}{\partial t} \omega(t) =& -t^{-2} \mathcal{1}_{[0,1]} + t^{-1} \delta_t \notag \\ =& t^{-1} ( - \omega (t) + \delta_t) \\ \frac{\partial}{\partial t} p_n (t, y) &= - (2n+1)^{\frac{1}{2}} 2y t^{-2} P_n' (\frac{2}{t} y - 1 ) \notag \\ =& - (2n+1)^{\frac{1}{2}} t^{-1} ((\frac{2}{t} y - 1) + 1 ) P_n' (\frac{2}{t} y - 1 ) \notag \\ =& - (2n+1)^{\frac{1}{2}} t^{-1} (x + 1 ) P_n' (x) \ (\because x = \frac{2}{t} y - 1) \notag \\ =& - (2n+1)^{\frac{1}{2}} t^{-1} [nP_n (x) + (2n-1) P_{n-1}(x) + (2n-3) P_{n-2}(x) + \ldots ] \ (\because 式(24)を用いた。) \notag \\ =& - (2n+1)^{\frac{1}{2}} t^{-1} [n (2n+1)^{-\frac{1}{2}} p_n (t, y) + (2n-1)^{\frac{1}{2}} p_{n-1} (t, y) + (2n-3)^{\frac{1}{2}} p_{n-2}(t, y) + \ldots ] \notag \\ & (\because 式(20)より、P_n (x) = (2n+1)^{-\frac{1}{2}} p_n (t, y) ) \end{align}$$

これらを用いて、式(12)のODEの右辺を計算する。

まず、第一項を求める。

$$∫0tu(y)(∂∂tpn(t,y))ω(t,y)dy≃∫0tg(t,y)(∂∂tpn(t,y))ω(t,y)dy=∫0t(∑n=0N−1xn(t)pn(t,y)){−(2n+1)12t−1[n(2n+1)−12pn(t,y)+(2n−1)12pn−1(t,y)+(2n−3)12pn−2(t,y)+…]}ω(t,y)dy(∵式(26)を用いた。)=−(2n+1)12t−1[n(2n+1)−12xn(t)+(2n−1)12xn−1(t)+(2n−3)12xn−2(t)+…](∵式(21)の正規直交性を用いた。) \begin{align} &\int^t_0 u(y) (\frac{\partial}{\partial t} p_n (t, y)) \omega (t, y) dy \notag \\ \simeq & \int^t_0 g (t, y) (\frac{\partial}{\partial t} p_n (t, y)) \omega (t, y) dy \notag \\ =&\int^t_0 ( \sum_{n=0}^{N-1} x_n (t) p_n (t, y) ) \{ - (2n+1)^{\frac{1}{2}} t^{-1} [n (2n+1)^{-\frac{1}{2}} p_n (t, y) + (2n-1)^{\frac{1}{2}} p_{n-1} (t, y) + (2n-3)^{\frac{1}{2}} p_{n-2}(t, y) + \ldots ] \} \omega (t, y) dy \notag \\ &(\because 式(26)を用いた。) \notag \\ =& - (2n+1)^{\frac{1}{2}} t^{-1} [n (2n+1)^{-\frac{1}{2}} x_n (t) + (2n-1)^{\frac{1}{2}} x_{n-1} (t) + (2n-3)^{\frac{1}{2}} x_{n-2}(t) + \ldots ] \\ &(\because 式(21)の正規直交性を用いた。) \notag \\ \end{align}$$

つぎに、第二項を求める。

$$∫0tu(y)pn(t,y)(∂∂tω(t,y))dy=∫0tu(y)pn(t,y)(t−1(−ω(t)+δt))dy (∵式(25)を用いた。)=−t−1∫0tu(y)pn(t,y)ω(t)dy+t−1∫0tu(y)pn(t,y)δtdy=−t−1xn(t)+t−1u(t)pn(t,t)=−t−1xn(t)+t−1(2n+1)1/2u(t) \begin{align} &\int^t_0 u(y) p_n (t, y) (\frac{\partial}{\partial t} \omega (t, y) ) dy \notag \\ =&\int^t_0 u(y) p_n (t, y) ( t^{-1} ( - \omega (t) + \delta_t) ) dy \ (\because 式(25)を用いた。) \notag \\ =& -t^{-1} \int^t_0 u(y) p_n (t, y) \omega (t) dy + t^{-1} \int^t_0 u(y) p_n (t, y) \delta_t dy \notag \\ =& -t^{-1} x_n (t) + t^{-1} u(t) p_n (t, t) \notag \\ =& -t^{-1} x_n (t) + t^{-1} (2n+1)^{1/2} u(t) \end{align}$$

式(12)の右辺に式(27), (28)を代入して、以下が成り立つ。

$$x˙n(t)=−(2n+1)12t−1[n(2n+1)−12xn(t)+(2n−1)12xn−1(t)+(2n−3)12xn−2(t)+…]−t−1xn(t)+t−1(2n+1)1/2u(t)=−(2n+1)12t−1[(n+1)(2n+1)−12xn(t)+(2n−1)12xn−1(t)+(2n−3)12xn−2(t)+…]+t−1(2n+1)1/2u(t) \begin{align} \dot{x}_n (t) =& - (2n+1)^{\frac{1}{2}} t^{-1} [n (2n+1)^{-\frac{1}{2}} x_n (t) + (2n-1)^{\frac{1}{2}} x_{n-1} (t) + (2n-3)^{\frac{1}{2}} x_{n-2}(t) + \ldots ] \notag \\ &-t^{-1} x_n (t) + t^{-1} (2n+1)^{1/2} u(t) \notag \\ =& - (2n+1)^{\frac{1}{2}} t^{-1} [(n+1) (2n+1)^{-\frac{1}{2}} x_n (t) + (2n-1)^{\frac{1}{2}} x_{n-1} (t) + (2n-3)^{\frac{1}{2}} x_{n-2}(t) + \ldots ] \notag \\ & + t^{-1} (2n+1)^{1/2} u(t) \end{align}$$

式(29)をベクトル化すると、以下の式(30)が成立する。

$$x˙(t)=−1tAx(t)+1tBu(t)Ank={(2n+1)1/2(2k+1)1/2if n>kn+1if n=k0if n<kBn=(2n+1)12 \begin{align} \dot{x}(t) &= - \frac{1}{t} A x(t) + \frac{1}{t} B u(t) \\ A_{nk} &= \left\{ \begin{array}{ll} (2n+1)^{1/2} (2k+1)^{1/2} & \text{if} \ n>k \\ n+1 & \text{if} \ n=k \\ 0 & \text{if} \ n<k \\ \end{array} \right. \\ B_n &= (2n+1)^{\frac{1}{2}} \\ \end{align}$$

式(31)の行列$AA$は、LSSL、S4、H3、mambaなどの後続研究では、HiPPO行列と呼ばれる。

式(21)の直交性を用いると、式(11)より、

$$xn(t)=∫−∞∞u(y)pn(t,y)ω(t,y)dy=∫0tu(y)pn(t,y)ω(t,y)dy \begin{align} x_n (t) &= \int^{\infty}_{-\infty} u(y) p_n (t, y) \omega(t, y) dy \notag \\ &= \int^t_0 u(y) p_n (t, y) \omega(t, y) dy \end{align}$$

としてみる。

ライプニッツの積分法則

$$∂∂t∫α(t)β(t)f(t,y)dy=∫α(t)β(t)∂∂tf(t,y)dy−α′(t)f(α(t),t)+β′(t)f(β(t),t) \begin{align*} \frac{\partial}{\partial t} \int^{\beta (t)}_{\alpha (t)} f(t, y) dy = \int^{\beta (t)}_{\alpha (t)} \frac{\partial}{\partial t} f(t, y) dy - \alpha' (t) f(\alpha(t), t) + \beta' (t) f( \beta(t), t) \end{align*}$$

について考え、ここで、$β(t)=t\beta (t) = t$、$α(t)=0\alpha (t) = 0$とすると、

$$ddt∫0tf(t,y)dy=f(t,t)ddtt−F(t,0)ddt0+∫0t∂∂tf(t,y)dy=∫0t∂∂tf(t,y)dy+f(t,t) \begin{align*} \frac{d}{dt} \int^t_0 f(t, y) dy &= f(t, t) \frac{d}{dt} t - F(t,0)\frac{d}{dt} 0 + \int^t_0 \frac{\partial}{\partial t} f(t, y) dy \\ &= \int^t_0 \frac{\partial}{\partial t} f(t, y) dy + f(t, t) \end{align*}$$

$f(ty)=u(y)pn(t,y)ω(t,y)f (t y) = u(y) p_n (t, y) \omega(t, y)$として、式(33)の両辺を微分すると、

$$x˙n(t)=ddt∫0tu(y)pn(t,y)ω(t,y)dy=ddt∫0tf(t,y)dy=∫0t∂∂tf(t,y)dy+f(t,t)(∵ライプニッツの積分法則）=∫0tu(y)(∂∂tpn(t,y))ω(t,y)dy+∫0tu(y)pn(t,y)(∂∂tω(t,y))dy+u(t)pn(t,t)ω(t,t) \begin{align*} \dot{x}_n (t) &= \frac{d}{dt} \int^t_0 u(y) p_n (t, y) \omega(t, y) dy \\ &= \frac{d}{dt} \int^t_0 f(t, y) dy \\ &= \int^t_0 \frac{\partial}{\partial t} f(t, y) dy + f(t, t) \quad (\because ライプニッツの積分法則） \\ &= \int^t_0 u(y) (\frac{\partial}{\partial t} p_n (t, y)) \omega (t, y) dy + \int^t_0 u(y) p_n (t, y) (\frac{\partial}{\partial t} \omega (t, y) ) dy \\ & \quad + u(t) p_n(t,t) \omega (t, t) \end{align*}$$

となり、式(12)と比較すると、第三項が増える。

実は、HiPPO-LegSにおいて、この第三項を計算すると、$u(t)pn(t,t)ω(t,t)=(2n+1)12t−1u(t)u(t) p_n(t,t) \omega (t, t) = (2n+1)^{\frac{1}{2}} t^{-1} u(t)$となり、式(32)は$Bn=2(2n+1)12B_n = 2(2n+1)^{\frac{1}{2}}$となると思われるが、実際にHiPPOを用いた実験では式(32)でうまくいってるので、これはおそらく違うのだろうが、それがなぜかわからない。

https://github.com/HazyResearch/hippo-code/issues/11

HiPPO-LegSは時間スケールに依存しない。

式(30)のODEを計算するにあたり、式(6)を求めたように両辺を積分し、GBT（$α=12\alpha=\frac{1}{2}$）を使って離散化して、$xn(t)x_n (t)$の更新式を得る。

$$x(t+Δt)−x(t)=∫tt+Δt−1tAx(t)+1tBu(t)dt≃Δt2(−1tAc(t)+1tBu(t))+(−1t+ΔtAc(t+Δt)+1t+ΔtBu(t+Δt)) \begin{align} x(t+\Delta t) - x(t) &= \int^{t+\Delta t}_t - \frac{1}{t} A x(t) + \frac{1}{t} B u(t) dt \notag \\ &\simeq \frac{\Delta t}{2}{(-\frac{1}{t}Ac(t) + \frac{1}{t}Bu(t))+(-\frac{1}{t+\Delta t}Ac(t+\Delta t) + \frac{1}{t+\Delta t}Bu(t+\Delta t))} \notag \end{align}$$

ゼロ次ホールドを仮定し、$u(t)=u(t+Δt)u(t) = u(t+\Delta t)$であるとすると、

$$x(t+Δt)−x(t)=Δt2(−1tAx(t)+1tBu(t))+Δt2(−1t+ΔtAx(t+Δt)+1t+ΔtBu(t))x(t+Δt)+Δt2(t+Δt)Ax(t+Δt)=x(t)−Δt2tAx(t)+(Δt2t+Δt2(t+Δt))Bu(t)(I+Δt2(t+Δt)A)x(t+Δt)=(I−Δt2tA)x(t)+(Δt2t+Δt2(t+Δt))Bu(t) \begin{align} x(t+\Delta t) - x(t) &= \frac{\Delta t}{2}(-\frac{1}{t}Ax(t) + \frac{1}{t}Bu(t)) + \frac{\Delta t}{2}(-\frac{1}{t+\Delta t}Ax(t+\Delta t) + \frac{1}{t+\Delta t}Bu(t)) \notag \\ x(t+\Delta t) + \frac{\Delta t}{2(t+\Delta t)}Ax(t+\Delta t) &= x(t) - \frac{\Delta t}{2t} Ax(t) + (\frac{\Delta t}{2t} + \frac{\Delta t}{2(t+\Delta t)})Bu(t) \notag \\ (I + \frac{\Delta t}{2(t+\Delta t)}A) x(t+\Delta t) &= (I - \frac{\Delta t}{2t} A)x(t) + (\frac{\Delta t}{2t} + \frac{\Delta t}{2(t+\Delta t)})Bu(t) \end{align}$$

ここで、時間を離散化し、$t=kΔtt=k\Delta t$とし、$xk:=x(kΔt)x_k := x(k\Delta t)$、$uk:=u(kΔt)u_k := u(k \Delta t)$とすれば、

$$(I+Δt2(kΔt+Δt)A)xk+1=(I−Δt2kΔtA)xk+(Δt2kΔt+Δt2(kΔt+Δt))Buk \begin{align*} (I + \frac{\Delta t}{2(k\Delta t+\Delta t)}A) x_{k+1} &= (I - \frac{\Delta t}{2k\Delta t} A)x_k + (\frac{\Delta t}{2k\Delta t} + \frac{\Delta t}{2(k\Delta t+\Delta t)})Bu_k \end{align*}$$

従って、以下が成り立つ。

$$(I+12(k+1)A)xk+1=(I−12kA)xk+(12k+12(k+1))Buk \begin{align} (I + \frac{1}{2(k+1)}A) x_{k+1} &= (I - \frac{1}{2k} A)x_k + (\frac{1}{2k} + \frac{1}{2(k+1)})Bu_k \end{align}$$

式(35)では、時間スケール$Δt\Delta t$がなくなっている。すなわち、HiPPO-LegSは時間スケールに依存しない。これはHiPPO-LegSが持つ特殊な特徴である。

この式に従って$xn(t)x_n (t)$を更新する。入力信号の履歴$u≤tu_{\leq t}$は、各時刻で式(9)に従って以下のように再構成することが出来る。

$$u≤t≃g(t)=∑n=0N−1xn(t)pn(t,y)=∑n=0N−1xn(t)(2n+1)12Pn(2yt−1) \begin{align*} u_{\leq t} \simeq g^{(t)} &= \sum_{n=0}^{N-1} x_n (t) p_n (t, y) \\ &= \sum_{n=0}^{N-1} x_n (t) (2n+1)^{\frac{1}{2}} P_n (\frac{2y}{t} - 1) \end{align*}$$

LSSL(Linear State-Space Layers)

ここまで、HiPPOおよびHiPPO-LegSについて説明してきた。HiPPO-LegSでは、式(30)のODEが導かれ、式(35)のように離散化された。LSSLやその発展であるS4、H3、mambaなどでは、式(1),(2)の線形システムを考え、式(6)、(7)の形に離散化するが、行列$AA$はHiPPO-LegSから導かれ、HiPPO行列と呼ばれる式(31)と同じになるように初期化する。このように初期化することで、$x(t)x(t)$に入力信号$u≤tu_{\leq t}$の履歴を記憶することを可能とし、MNISTベンチマークでの性能を60%から98%に上昇させることが出来る。ここでは、これらの式(6), (7)から導かれるLSSLの性質について見ていく。

畳み込みによる高速化

式(6),(7)を見た時に思ったかもしれないが、式(6)を式(7)に代入することで、以下のように$xx$を用いずに$u≤tu_{\leq t}$のみで出力$yky_k$を計算することが出来る。ただし、$x−1=0x_{-1} = 0$とする。

$$yk=Cxk+Duk=C(Aˉxk−1+Bˉuk)+Duk=CAˉ(Aˉxk−2+Bˉuk−1)+CBˉuk+Duk=C(Aˉ)2(Aˉxk−3+Bˉuk−2)+CAˉBˉuk−1+CBˉuk+Duk=C(Aˉ)3(Aˉxk−4+Bˉuk−3)+C(Aˉ)2Bˉuk−2+CAˉBˉuk−1+CBˉuk+Duk⋮=C(Aˉ)kx0+⋯+C(Aˉ)2Bˉuk−2+CAˉBˉuk−1+CBˉuk+Duk=C(Aˉ)k(Aˉx−1+Bˉu0)+⋯+C(Aˉ)2Bˉuk−2+CAˉBˉuk−1+CBˉuk+Duk=C(Aˉ)kBˉu0+⋯+C(Aˉ)2Bˉuk−2+CAˉBˉuk−1+CBˉuk+Duk(∵x−1=0) \begin{align*} y_k =& C x_k + D u_k \\ =& C (\bar{A} x_{k-1} + \bar{B} u_k) + D u_k \\ =& C \bar{A} (\bar{A} x_{k-2} + \bar{B} u_{k-1}) + C \bar{B} u_k + D u_k \\ =& C (\bar{A})^2 (\bar{A} x_{k-3} + \bar{B} u_{k-2}) + C \bar{A} \bar{B} u_{k-1} + C \bar{B} u_k + D u_k \\ =& C (\bar{A})^3 (\bar{A} x_{k-4} + \bar{B} u_{k-3}) + C (\bar{A})^2 \bar{B} u_{k-2} + C \bar{A} \bar{B} u_{k-1} + C \bar{B} u_k + D u_k \\ \vdots \\ =& C (\bar{A})^k x_0 + \cdots + C (\bar{A})^2 \bar{B} u_{k-2} + C \bar{A} \bar{B} u_{k-1} + C \bar{B} u_k + D u_k \\ =& C (\bar{A})^k (\bar{A} x_{-1} + \bar{B} u_0) + \cdots + C (\bar{A})^2 \bar{B} u_{k-2} + C \bar{A} \bar{B} u_{k-1} + C \bar{B} u_k + D u_k \\ =& C (\bar{A})^k \bar{B} u_0 + \cdots + C (\bar{A})^2 \bar{B} u_{k-2} + C \bar{A} \bar{B} u_{k-1} + C \bar{B} u_k + D u_k \\ &(\because x_{-1} = 0) \end{align*}$$

従って、$Kk(Aˉ,Bˉ,C)=(CBˉ,CAˉBˉ,…,C(Aˉ)kBˉ)\mathcal{K}_k (\bar{A}, \bar{B}, C) = (C \bar{B}, C \bar{A} \bar{B}, \ldots, C (\bar{A})^k \bar{B})$として、以下が成り立つ。

$$yk=(CBˉ,CAˉBˉ,…,C(Aˉ)kBˉ)∗u≤k+Duk=Kk(Aˉ,Bˉ,C)∗u≤k+Duk \begin{align} y_k &= (C \bar{B}, C \bar{A} \bar{B}, \ldots, C (\bar{A})^k \bar{B}) * u_{\leq k} + D u_k \notag \\ &= \mathcal{K}_k (\bar{A}, \bar{B}, C) * u_{\leq k} + D u_k \end{align}$$

ここで、以下のように置いた。

$$Kk(Aˉ,Bˉ,C)=(CAˉiB)i=0,…,k=(CBˉ,CAˉBˉ,…,C(Aˉ)kBˉ)u≤k=(u0,u1,…,uk) \begin{align} \mathcal{K}_k (\bar{A}, \bar{B}, C) &= (C \bar{A}^i B)_{i = 0, \ldots, k} \notag \\ &= (C \bar{B}, C \bar{A} \bar{B}, \ldots, C (\bar{A})^k \bar{B}) \\ u_{\leq k} &= (u_0, u_1, \ldots, u_k) \end{align}$$

また、系列長$LL$の入力信号$u≤L−1=(u0,u1,…,uL−1)u_{\leq L-1} = (u_0, u_1, \ldots, u_{L-1})$が与えられたとき、出力$yL−1y_{L-1}$は、以下のように表される。

$$yL−1=KL−1(Aˉ,Bˉ,C)∗u≤L−1+DuL−1 \begin{align} y_{L-1} = \mathcal{K}_{L-1} (\bar{A}, \bar{B}, C) * u_{\leq L-1} + D u_{L-1} \end{align}$$

このように畳み込み形式で計算することで、全時刻にわたる$y∈RH×Ly \in \mathbb{R}^{H \times L}$を三度のFFTで一度に求めることが出来る。

ただし、LSSLの計算にはボトルネックがある。式(6)の再帰形式では、離散化された状態行列$AA$、つまり$Aˉ\bar{A}$の、行列-ベクトル乗算（matrix-vector multiplication; MVM）、そして式(35)の畳み込み形式では、行列$Aˉ\bar{A}$のべき乗を多く含むKrylov function $KL\mathcal{K}_L$が、それぞれボトルネックとなる。これらをうまく解決し、効率よく計算するアルゴリズムがLSSLでも提案されているが、S4でより効率的な手法が提案されているため、ここでは紹介しない。（S4については、次の記事で説明する。）

LSSLの表現力

ここまで、LSSLが式(6)、(7)の再帰形式、式(39)の畳み込み形式の二つの形式を持つことを見てきた。LSSLでは線形計算のみを行い、RNNやTransformerなどで見られるような非線形関数は用いていない。それでもLSSLは高い表現力を持つことを示す。

畳み込みはLSSLである。

Under Construction

RNNはLSSLである。

Under Construction

Lemma [3, p.6, 3.1]

一次元ゲート付き再帰$xt=(1−σ(z))xt−1+σ(z)ut\bold{x_t = (1 - \sigma(z)) x_{t-1} + \sigma(z) u_t}$（ただし、$σ\bold{\sigma}$はシグモイド関数で$z\bold{z}$は任意の表現）は、1次元の線形ODE$x˙(t)=−x(t)+u(t)\bold{\dot{x}(t) = -x(t) + u(t)}$を$α=1\bold{\alpha=1}$のGBT（つまり、後退オイラー法）で離散化したものとみなすことが出来る。

証明：
$x˙(t)=−x(t)+u(t)\dot{x}(t) = -x(t) + u(t)$を後退オイラー法で離散化すると、以下が成り立つ。

$$xt−1=xt−Δt(−xt+ut) \begin{align*} x_{t-1} = x_t - \Delta t (- x_t + u_t) \end{align*}$$

ここで、$Δt=ez\Delta t = e^z$とすると、

$$xt−1=xt−ez(−xt+ut)=(1+ez)xt−ezut \begin{align*} x_{t-1} = x_t - e^z (- x_t + u_t) = (1 + e^z)x_t - e^z u_t \end{align*}$$

従って、以下が成り立つ。

$$xt=11+ezxt−1+ez1+ezut=e−z1+e−zxt−1+11+e−zut=(1−σ(z))xt−1+σ(z)ut(∵σ(z)=11+e−z) \begin{align} x_t &= \frac{1}{1 + e^z} x_{t-1} + \frac{e^z}{1 + e^z} u_t \notag \\ &= \frac{e^{-z}}{1 + e^{-z}} x_{t-1} + \frac{1}{1 + e^{-z}} u_t \notag \\ &= (1 - \sigma(z)) x_{t-1} + \sigma(z) u_t \quad (\because \sigma(z) = \frac{1}{1 + e^{-z}}) \end{align}$$(証明終)

Lemma [3, p.6, 3.2]

Under Construction

Deep LSSLs

Under Construction

参考文献

[1] Voelker, A., Kajić, I., & Eliasmith, C. (2019). Legendre memory units: Continuous-time representation in recurrent neural networks. Advances in neural information processing systems, 32.
[2] Gu, A., Dao, T., Ermon, S., Rudra, A., & Ré, C. (2020). Hippo: Recurrent memory with optimal polynomial projections. Advances in neural information processing systems, 33, 1474-1487.
[3] Gu, A., Johnson, I., Goel, K., Saab, K., Dao, T., Rudra, A., & Ré, C. (2021). Combining recurrent, convolutional, and continuous-time models with linear state space layers. Advances in neural information processing systems, 34, 572-585.
[4] Gu, A., Goel, K., & Ré, C. (2021). Efficiently modeling long sequences with structured state spaces. arXiv preprint arXiv:2111.00396.
[5] Gu, A., Johnson, I., Timalsina, A., Rudra, A., & Ré, C. (2022). How to train your hippo: State space models with generalized orthogonal basis projections. arXiv preprint arXiv:2206.12037.
[6] 角田良太郎「HiPPO/S4解説」, https://techblog.morphoinc.com/entry/2022/05/24/102648, 2024年2月10日閲覧
[7] Guofeng Zhang, Tongwen Chen, and Xiang Chen. Performance recovery in digital implementation of analogue systems. SIAM journal on control and optimization, 45(6):2207–2223, 2007.
[8] 岡島寛「状態フィードバック制御・状態方程式に基づく制御のまとめ」、制御工学の教科書、2024年2月17日公開、2024年2月21日閲覧