HiPPO - 3

이전 게시글(https://pasus.tistory.com/363)에서 함수 \(f(x), \ 0 \lt x \le t\) 를 N차원 부분 함수공간으로 투사한 근사 함수 \(g(x), \ 0 \lt x \le t\) 를 다음과 같이 유도하였다.

 

\[ \begin{align} & g(x)= \sum_{n=0}^{N-1} c_n (t) \sqrt{(2n+1)} P_n \left( \frac{2x}{t}-1 \right) \tag{1} \\ \\ & \dot{\mathbf{c}}(t)= - \frac{1}{t} A \mathbf{c}(t)+ \frac{1}{t} Bf(t) \end{align} \]

 

식 (1)에 의하면 HiPPO는 본질적으로 연속시간(continuous-time) 에서 정의된 상미분 방정식(ODE)을 기반으로 데이터의 상태를 업데이트한다. 그러나 센서 신호, 주가 데이터 등의 실제 데이터는 대부분 샘플링(sampling)으로 수집되므로, HiPPO를 실제 시스템에 적용하려면 연속시간 히포 미분방정식(continuous-time HiPPO ODE)을 이산시간(discrete-time) 차분 방정식(difference equation)으로 바꿀 필요가 있다.

 

 

히포 미분방정식은 시간에 따라 변화하는 시변(time-varying) 시스템으로, 시간 \(t\) 에 따라 행렬 \(\frac{1}{t} A\) 와 \(\frac{1}{t} B\) 가 변화한다. 따라서 일반적인 선형 시불변(linear time-invariant, LTI) 시스템과는 달리 이산화 과정에서 주의가 필요하다. 그러나 논문에서 언급된 forward Euler, backward Euler, bilinear 등의 방법은 시변 시스템에도 적용 가능하다. 하지만 ZOH (Zero-Order Hold) 방법의 경우에는 샘플링 구간 동안 행렬을 고정시키든가 또는 평균값을 사용해야 한다(https://pasus.tistory.com/321).

여기서는 forward Euler, backward Euler, bilinear에 대해서만 알아보고자 한다.

미분을 다음과 같이 근사화하는 방법을 오일러 근사(Euler's approximation)법이라고 한다.

 

\[ \begin{align} \dot{\mathbf{c}}(t_k ) & \approx \frac{ \mathbf{c}(t_{k+1} )-\mathbf{c}(t_k )}{\Delta t} \tag{2} \\ \\ \dot{\mathbf{c}}(t_k ) & \approx \frac{ \mathbf{c}(t_k )-\mathbf{c}(t_{k-1} )}{\Delta t} \tag{3} \end{align} \]

 

여기서 \(t=t_k\) 는 \(k\) 번째 샘플링 싯점으로서 \(t_k=k \Delta t\) 이며, \(\Delta t=t_{k+1}-t_k\) 는 샘플링 간격이다. \(k\) 는 시간스텝(time step)을 나타내는 인덱스로서 정수값을 갖는다. 식 (2)는 \(\dot{\mathbf{c}}(t_k )\) 를 근사화할 때 미래시간과 현재시간의 \(\mathbf{c}\) 값을 이용했으므로 forward Euler 근사라고 하고, 식 (3)은 현재시간과 과거시간의 \(\mathbf{c}\) 값을 이용했으므로 backward Euler 근사라고 한다.

식 (2)를 (1)에 대입하면 다음과 같다.

 

\[ \begin{align} \dot{\mathbf{c}}(t_k ) \approx \frac{ \mathbf{c}(t_{k+1} )-\mathbf{c}(t_k )}{\Delta t} = - \frac{1}{t_k} A \mathbf{c}(t_k )+ \frac{1}{t_k} Bf(t_k ) \tag{4} \end{align} \]

 

위 식은 다음과 같이 정리된다.

 

\[ \begin{align} \mathbf{c}(t_{k+1} )= \left( I- \frac{\Delta t}{t_k} A \right) \mathbf{c}(t_k )+ \frac{\Delta t}{t_k} Bf(t_k ) \tag{5} \end{align} \]

 

이산시간 시스템의 기호를 따르기 위하여 다음과 같이 정의하면,

 

\[ \begin{align} \mathbf{c}_{k+1}= \mathbf{c} (t_{k+1} ), \ \ \ \ \mathbf{c}_k= \mathbf{c}(t_k ), \ \ \ \ f_k=f(t_k ), \ \ \ \ \frac{\Delta t}{t_k} = \frac{\Delta t}{k \Delta t}= \frac{1}{k} \tag{6} \end{align} \]

 

다음과 같이 식 (1)에 대한 forward Euler 방식의 이산시간 등가 모델을 얻을 수 있다.

 

\[ \begin{align} \mathbf{c}_{k+1} = \left( I- \frac{1}{k} A \right) \mathbf{c}_k+ \frac{1}{k} Bf_k \tag{7} \end{align} \]

 

이번에는 식 (3)을 (1)에 대입하면 다음과 같다.

 

\[ \begin{align} \dot{\mathbf{c}}(t_k ) \approx \frac{\mathbf{c}(t_k )- \mathbf{c}(t_{k-1} )}{\Delta t}= - \frac{1}{t_k} A \mathbf{c}(t_k )+ \frac{1}{t_k} Bf(t_k ) \tag{8} \end{align} \]

 

위 식을 정리하면 다음과 같다.

 

\[ \begin{align} \mathbf{c}(t_k )= \left( I+ \frac{\Delta t}{t_k} A \right)^{-1} \mathbf{c}(t_{k-1} )+ \left( I+ \frac{\Delta t}{t_k} A \right)^{-1} \frac{\Delta t}{t_k} Bf(t_k ) \tag{9} \end{align} \]

 

이산시간 시스템의 기호를 따르면 다음과 같이 식 (1)에 대한 backward Euler 방식의 이산시간 등가 모델을 얻을 수 있다.

 

\[ \begin{align} \mathbf{c}_k = \left( I+ \frac{1}{k} A \right)^{-1} \mathbf{c}_{k-1} + \left( I+ \frac{1}{k} A \right)^{-1} \frac{1}{k} Bf_k \tag{10} \end{align} \]

 

식 (7)과 식 (10)은 \(f_k\) 가 가해지는 싯점에 차이가 있다. 식 (7)에서는 \(\mathbf{c}_{k+1}\) 를 계산할 때 \(f_k\) 가 적용된 반면 식 (10)에서는 \(\mathbf{c}_k\) 를 계산할 때 \(f_k\) 가 적용되었다. 제어 분야에서는 식 (7)의 형식이 일반적으로 사용되는 반면에, RNN 계열에서는 식 (10)의 형식이 사용된다. 식 (10) 은 \(N \times N\) 행렬의 역행렬을 계산해야 하므로 투사(projection)된 차원이 클 경우 계산상의 문제를 야기할 수 있다.

 

 

Bilinear 변환은 \(t_k\) 와 \(t_{k+1}\) 또는 \(t_k\) 와 \(t_{k-1}\) 사이를 평균화하는 방법이다. 식 (1)의 \(\mathbf{c}(t)\) 항을 샘플링 싯점 \(t=t_k\) 에서 다음과 같이 근사화한다.

 

\[ \begin{align} -\frac{1}{t} A \mathbf{c}(t) & \approx -\frac{1}{2} \left( \frac{1}{t_{k+1}} A \mathbf{c}(t_{k+1}) + \frac{1}{t_k} A \mathbf{c}(t_k ) \right) \tag{11} \\ \\ -\frac{1}{t} A \mathbf{c}(t) & \approx -\frac{1}{2} \left( \frac{1}{t_k} A \mathbf{c}(t_k) + \frac{1}{t_{k-1}} A \mathbf{c}(t_{k-1} ) \right) \tag{12} \end{align} \]

 

식 (11)은 미래시간과 현재시간의 값을 이용한 forward 평균 근사이고, 식 (12)는 현재시간과 과거시간의 값을 이용한 backward 평균 근사다.

식 (11)과 (2)를 식 (1)에 대입하면 다음과 같다.

 

\[ \begin{align} \frac{ \mathbf{c}(t_{k+1} )- \mathbf{c}(t_k )}{\Delta t} = - \frac{1}{2} \left( \frac{1}{t_{k+1}} A\mathbf{c}(t_{k+1} )+ \frac{1}{t_k} A \mathbf{c}(t_k ) \right)+ \frac{1}{t_k} Bf(t_k ) \tag{13} \end{align} \]

 

위 식은 다음과 같이 정리된다.

 

\[ \begin{align} \mathbf{c}(t_{k+1} )= \left( I+ \frac{\Delta t}{2t_{k+1}} A \right)^{-1} \left( I- \frac{\Delta t}{2t_k} A \right) \mathbf{c}(t_k )+ \left( I+ \frac{\Delta t}{2t_{k+1}} A \right)^{-1} \frac{\Delta t}{t_k} Bf(t_k ) \tag{14} \end{align} \]

 

이산시간 시스템의 기호를 따르면 다음과 같이 식 (1)에 대한 forward bilinear 방식의 이산시간 등가 모델을 얻을 수 있다.

 

\[ \begin{align} \mathbf{c}_{k+1} = \left( I+ \frac{1}{2(k+1)} A \right)^{-1} \left( I- \frac{1}{2k} A \right) \mathbf{c}_k+ \left( I+ \frac{1}{2(k+1)} A \right)^{-1} \frac{1}{k} Bf_k \tag{15} \end{align} \]

 

이번에는 식 (12)와 (3)을 식 (1)에 대입하면 다음과 같다.

 

\[ \begin{align} \frac{ \mathbf{c}(t_k )- \mathbf{c}(t_{k-1} )}{\Delta t} = - \frac{1}{2} \left( \frac{1}{t_k} A\mathbf{c}(t_k )+ \frac{1}{t_{k-1}} A \mathbf{c}(t_{k-1} ) \right)+ \frac{1}{t_k} Bf(t_k ) \tag{16} \end{align} \]

 

위 식은 다음과 같이 정리된다.

 

\[ \begin{align} \mathbf{c}(t_k )= \left( I+ \frac{\Delta t}{2t_k} A \right)^{-1} \left( I- \frac{\Delta t}{2t_{k-1}} A \right) \mathbf{c}(t_k )+ \left( I+ \frac{\Delta t}{2t_k} A \right)^{-1} \frac{\Delta t}{t_k} Bf(t_k ) \tag{17} \end{align} \]

 

이산시간 시스템의 기호를 따르면 다음과 같이 식 (1)에 대한 backward bilinear 방식의 이산시간 등가 모델을 얻을 수 있다.

 

\[ \begin{align} \mathbf{c}_k = \left( I+ \frac{1}{2k} A \right)^{-1} \left( I- \frac{1}{2(k-1)} A \right) \mathbf{c}_{k-1}+ \left( I+ \frac{1}{2k} A \right)^{-1} \frac{1}{k} Bf_k \tag{18} \end{align} \]

 

식 (15)와 식 (18)은 \(f_k\) 가 가해지는 싯점에 차이가 있다. 식 (15)에서는 \(\mathbf{c}_{k+1}\) 를 계산할 때 \(f_k\) 가 적용된 반면 식 (18)에서는 \(\mathbf{c}_k\) 를 계산할 때 \(f_k\) 가 적용되었다. 또한 식 (15)와 (18)은 \(N \times N\) 행렬의 역행렬을 계산해야 하므로 투사(projection)된 차원이 클 경우 계산상의 문제를 야기할 수 있다.

계수 차분방정식 (8), (10), (15), (18)에 의하면 수식 상에 샘플링 간격 \(\Delta t\) 가 없으므로 샘플링 간격에 의존하지 않는다는 것을 알 수 있다. 또한 HiPPO는 시간 스케일에 독립적(timescale invariant)이며 어떤 시간 범위에서도 성능을 유지할 수 있다. 시간 스케일에 의존적이면 초당 1회 데이터를 수집하던 센서가 초당 10회 데이터를 수집하면 성능이 저하될 수도 있다.

HiPPO가 시간 스케일에 독립적이라는 것은 다음과 같이 증명할 수 있다. 함수 \(f(x)\) 에서 시간 스케일이 변경된 함수를 \(\tilde{f}(x)=f(\alpha x)\) 라고 하자. 그러면 계수 적분방정식에 의해서 \(\tilde{f}(x)\) 의 계수 \(\tilde{c}_n (t)\) 는 다음과 같이 계산할 수 있다.

 

\[ \begin{align} \tilde{c}_n (t) &= \int_0^t \tilde{f}(x) \sqrt{(2n+1) } P_n \left( \frac{2x}{t}-1 \right) \frac{1}{t} \ dx \tag{19} \\ \\ &= \int_0^t f(\alpha x) \sqrt{(2n+1) } P_n \left( \frac{2 \alpha x}{\alpha t}-1 \right) \frac{1}{\alpha t} \alpha \ dx \\ \\ &= \int_0^{\alpha t} f(y) \sqrt{(2n+1) } P_n \left( \frac{2y}{\alpha t}-1 \right) \frac{1}{\alpha t} \ dy \\ \\ &=c_n (\alpha t) \end{align} \]

 

여기서 \(y= \alpha x\) 로 변수를 변환했다. 따라서 함수의 시간 스케일을 바꾸면 동일한 스케일로 계수도 바뀐다는 것을 알 수 있다.

HiPPO는 시간에 따라 변화하는 데이터를 직교 다항식 기반으로 최적화하여 식 (8), (10), (15) 또는 (18)을 이용하여 계수 \(\mathbf{c}_k= \mathbf{c}(t_k)\) 의 형태로 함축적으로 저장한다. 그리고 원래 함수 \(f(x)\) 는 \(\mathbf{c}(t_k)\) 를 이용하여 다음과 같은 근사 함수 \(g(x)\) 로 복원할 수 있다.

 

\[ \begin{align} g(x_i )= \sum_{n=0}^{N-1} c_n (t_k ) \sqrt{(2n+1) } P_n \left( \frac{2x_i}{t_k} -1 \right) , \ \ \ \ x_i \le t_k \tag{20} \end{align} \]

 

다음은 \(x \in [0, \ 150]\) 구간에서 \(\Delta t=0.1\), \(N=10\) 을 사용하여 \(f(x)=\cos \left(\frac{x}{20} \right) \sin \left( \frac{x}{5} \right )\) 을 HiPPO 로 근사한 것이다. 아래 그림은 \(t_k=75\) 에서 근사한 것이다.

 

 

다음 그림은 \(t_k=150\) 에서 근사한 것이다.

 

 

다음 그림은 \(N=20\) 을 사용하여 \(t_k=150\) 에서 근사한 것이다.

 

 

 

 

다음은 해당 매트랩 코드다.

 

% Comparison among HiPPO ODE, forward Euler, backward Euler, 
% and Bilinear forwad and backward methods
%
% (c) st.watermelon

clear; 

% HiPPO setup
N = 10; 
tf = 150;
dt = 0.1; % time interval
t_step = dt:dt:tf;

% HiPPO A, B
[A, B] = hippo_AB(N);

% input signal
f = @(t) cos(t/20).*sin(t/5);

ft_orig = f(t_step);
ft = ft_orig + 0.1 * randn(1,length(t_step));

% HiPPO ODE 
hippo_ode = @(t, c) (-A / t) * c + (B / t) * f(t);

% initial condition
c0 = zeros(N, 1);

% 0. HiPPO ODE
[time, c_ode] = ode45(hippo_ode, t_step, c0);


% 1. forward Euler
c_euler_fwd = zeros(N, 1);
c_euler_fwd_history = zeros(length(t_step), N); 
for k = 1:length(t_step)-1
    c_euler_fwd = (eye(N) - A/k) * c_euler_fwd + (B/k) * ft(k); 
    c_euler_fwd_history(k+1, :) = c_euler_fwd'; 
end

% 2. backward Euler
c_euler_bwd = zeros(N, 1);
c_euler_bwd_history = zeros(length(t_step), N);
for k = 2:length(t_step)
    M = eye(N) + A/k;
    invM = inv(M);
    c_euler_bwd = invM * c_euler_bwd + invM * (B/k) * ft(k); 
    c_euler_bwd_history(k, :) = c_euler_bwd'; 
end

% 3. bilinear forward 
c_bi_fwd = zeros(N, 1);
c_bi_fwd_history = zeros(length(t_step), N);
for k = 1:length(t_step)-1
    M = eye(N) + A/(2*(k+0));
    invM = inv(M);
    c_bi_fwd = invM * (eye(N) - A/(2*k)) * c_bi_fwd + invM * (B/k) * ft(k); 
    c_bi_fwd_history(k+1, :) = c_bi_fwd'; 
end

% 4. bilinear backward 
c_bi_bwd = zeros(N, 1);
c_bi_bwd_history = zeros(length(t_step), N);
for k = 2:length(t_step)
    M = eye(N) + A/(2*k);
    invM = inv(M);
    c_bi_bwd = invM * (eye(N) - A/(2*(k-0))) * c_bi_bwd + invM * (B/k) * ft(k); 
    c_bi_bwd_history(k, :) = c_bi_bwd'; 
end


%%
% fitting

cut = round(length(t_step)*1/2);
g = hippo_fit(c_bi_fwd_history(1:cut,:), dt);


% plotting
figure, plot(t_step, ft, '.', 'color', [.7 .7 .7]), hold on
plot(t_step, ft_orig, 'r')
plot(t_step(1:cut), g,'b', 'LineWidth', 1.2)
legend('noisy f(t)','f(t)','g(t)')
xlabel('Time');
title(['N = ' num2str(N)]);

%%

function g = hippo_fit(c, dt)
%
% g = hippo_fit(c, dt)
% HiPPO fitting
% input
%   c(t): coefficients
%   dt: sampling time interval
% output
%   g: fitted function
%
% (c) st.watermelon

[k, N] = size(c);

g = [];
for ii=1:k
    x = ii*dt;
    sum = 0;
    for n=0:N-1
        p=LegendrePoly(n);
        Pn=polyval(p,2*x/(k*dt)-1);
        sum = sum + c(k, n+1)*sqrt(2*n+1)*Pn;
    end
    g = [g; sum];
end

end


%%
function lpoly = LegendrePoly(n)
%
% lpoly = LegendrePoly(n)
% N-th order Legendre polynomials
%
% input: order of Legendre polynomials
% output: polynomial coefficients
%
% (c) st.watermelon

for ii=0:n
    
    if ii==0
        p{ii+1}=1;
    elseif ii==1
        p{ii+1}=[1 0];
    else
        p{ii+1}=(2*ii-1)/ii*[p{ii} 0]-(ii-1)/ii*[0 0 p{ii-1}];
    end
    
end

lpoly=p{n+1};

end