Next Generation Reservoir Computing 코드 실습 (MATLAB)

오늘은 MATLAB으로 Next Generation Reservoir Computing (NG-RC)를 구현해보려고 한다.

 

 

 

코드의 전체적인 구성은 다음과 같다.

1. Lorenz attractor 데이터셋 생성
2. Feature vector 구성
3. Ridge regression을 통해 가중치 계산
4. 시계열 예측

이번 실습에서는 Fig 2의 결과를 재현해보고자 한다.

즉, \( x,\,y,\,z \) 를 모두 넣어서 전체 시스템을 forecasting하는 것을 목표로 한다.

 

 

1. Lorenz attractor 데이터셋 생성 코드

clear; clc; close all;

%% 1. 데이터 생성 (Ground Truth: Lorenz System)
sigma = 10; rho = 28; beta = 8/3;
lorenz = @(t,s) [sigma*(s(2)-s(1)); s(1)*(rho-s(3))-s(2); s(1)*s(2)-beta*s(3)];

dt = 0.025; 
[t, raw_data] = ode45(lorenz, 0:dt:100, [1; 1; 1]);

% 데이터 전처리 (Z-score Normalization)
data = zscore(raw_data)'; % (Dimension x Time) 형태로 전치: 3 x N
[dim, total_len] = size(data);

% train/test 데이터 분할
train_len = 2000;
test_len = length(data(1,:))-train_len;
train_data = data(:, 1:train_len);
test_data = data(:, train_len+1:train_len+test_len);

 

예쁜 나비

 

 

 

NG-RC 원 논문을 보면 \( k \)와 \( s \) 두개의 parameter가 등장하는데, 각각 delay 길이와 stride를 의미한다.

%% 2. Feature vector 구성

% parameter setting
k = 2;  % Delay taps (현재와 바로 전 과거 1개만 봄)
s = 1;  % Stride
alpha = 1e-4; % Ridge regression parameter (Tikhonov)
bias = 1.0; % bias term

% Warm-up: k*s 만큼의 과거 데이터가 있어야 feature vector를 만들 수 있음
warmup = k*s;

버리는 데이터 길이인 warmup은 논문에서처럼 \( sk \) 로 설정했다.

 

 

 

%% 2. Feature vector 구성

O_train = [];
Y_target = []; % 다음 스텝 예측 목표 (X_t+1)

for i = warmup : train_len - 1
    
    % --- [Step A] Linear Part (Eq. 5) ---
    % 현재(t)와 과거(t-k)의 상태 벡터를 가져와서 일렬로 펼칩니다.
    % 현재로부터 k개를 가져오는 것이기 때문에 for loop의 인덱스는 0부터 시작
    vec_linear = [];
    for delay = 0:(k-1)
        vec_linear = [vec_linear; train_data(:, i - delay*s)]; 
    end
    
    % --- [Step B] Nonlinear Part (Eq. 6) ---
    % linear vector끼리의 outer product
    full_prod = vec_linear * vec_linear'; % Outer Product 
    
    % upper triangle 추출
    idx = triu(true(size(full_prod))); 
    vec_nonlinear = full_prod(idx);
    
    % Total Feature Vector ---
    % Bias + Linear + Nonlinear
    O_vec = [bias; vec_linear; vec_nonlinear];
    
    O_train = [O_train, O_vec];      % 입력 데이터 행렬
    Y_target = [Y_target, train_data(:, i+1)]; % 정답(타겟): 바로 다음 스텝
end

 

Step A에서는 논문의 5번 수식을 구현하였다.

 

 

Step B에서는 논문의 6번 수식을 구현( \( p=2\) 로 설정) & 모든 term을 합쳐주었다.

 

우리의 경우 \( k=2 \) 이고 \( d=3 \)이기 때문에 total feature vector가 1+6+21 = 28 차원이 나온다.

 

 

%% 3. Ridge regression을 통해 가중치 행렬 (W_out) 계산

feat_dim = size(O_train, 1);
W_out = (Y_target * O_train') *  inv(O_train * O_train' + alpha * eye(feat_dim));

 

 

 

 

% Lorenz attractor prediction

Y_pred = zeros(dim, test_len);

% train data의 마지막 부분(초기값)을 가져옵니다.
current_history = train_data(:, train_len - warmup : train_len);

for t = 1:test_len
    % 1. 현재 history buffer로 feature vector 생성
    vec_linear = [];
    for delay = 0:(k-1)
        % history의 끝에서부터 과거로 가져옴
        vec_linear = [vec_linear; current_history(:, end - delay*s)];
    end
    
    full_prod = vec_linear * vec_linear';
    vec_nonlinear = full_prod(triu(true(size(full_prod))));
    O_vec = [1; vec_linear; vec_nonlinear];
    
    % 2. 다음 스텝 예측 (행렬 곱)
    next_step = W_out * O_vec;
    Y_pred(:, t) = next_step;
    
    % 3. History 업데이트
    % 가장 오래된 데이터 버리고, 방금 예측한 값을 끝에 추가 (Closed-loop)
    current_history = [current_history(:, 2:end), next_step];
end

 

 

실제 값과 예측값 비교

 

%% 시각화

figure('Position', [100, 100, 1000, 500]);
plot(test_data(1, :), 'k', 'LineWidth', 1.5); hold on;
plot(Y_pred(1, :), 'r--', 'LineWidth', 1.5);
title('Time Series Prediction (x-component)');
legend('Ground Truth', 'NG-RC Prediction');
xlabel('Time Steps'); ylabel('x'); grid on;

논문의 그림과 비슷하게, 특정 Lyapunov time까지는 아주 잘 맞추는 것을 볼 수 있다.

 

 

 

<전체 코드>

 

clear; clc; close all;

%% 1. 데이터 생성 (Ground Truth: Lorenz System)
sigma = 10; rho = 28; beta = 8/3;
lorenz = @(t,s) [sigma*(s(2)-s(1)); s(1)*(rho-s(3))-s(2); s(1)*s(2)-beta*s(3)];

dt = 0.025; 
[t, raw_data] = ode45(lorenz, 0:dt:100, [1; 1; 1]);

% 데이터 전처리 (Z-score Normalization)
data = zscore(raw_data)'; % (Dimension x Time) 형태로 전치: 3 x N
[dim, total_len] = size(data);

% train/test 데이터 분할
train_len = 2000;
test_len = length(data(1,:))-train_len;
train_data = data(:, 1:train_len);
test_data = data(:, train_len+1:train_len+test_len);

%% 2. Feature vector 구성

% parameter setting
k = 2;  % Delay taps (현재와 바로 전 과거 1개만 봄)
s = 1;  % Stride
alpha = 1e-4; % Ridge regression parameter (Tikhonov)
bias = 1.0; % bias term

% Warm-up: k*s 만큼의 과거 데이터가 있어야 feature vector를 만들 수 있음
warmup = k*s;

O_train = [];
Y_target = []; % 다음 스텝 예측 목표 (X_t+1)

for i = warmup : train_len - 1
    
    % --- [Step A] Linear Part (Eq. 5) ---
    % 현재(t)와 과거(t-k)의 상태 벡터를 가져와서 일렬로 펼칩니다.
    % 현재로부터 k개를 가져오는 것이기 때문에 for loop의 인덱스는 0부터 시작
    vec_linear = [];
    for delay = 0:(k-1)
        vec_linear = [vec_linear; train_data(:, i - delay*s)]; 
    end
    
    % --- [Step B] Nonlinear Part (Eq. 6) ---
    % linear vector끼리의 outer product
    full_prod = vec_linear * vec_linear'; % Outer Product 
    
    % upper triangle 추출
    idx = triu(true(size(full_prod))); 
    vec_nonlinear = full_prod(idx);
    
    % Total Feature Vector ---
    % Bias + Linear + Nonlinear
    O_vec = [bias; vec_linear; vec_nonlinear];
    
    O_train = [O_train, O_vec];      % 입력 데이터 행렬
    Y_target = [Y_target, train_data(:, i+1)]; % 정답(타겟): 바로 다음 스텝
end

%% 3. Ridge regression을 통해 가중치 행렬 (W_out) 계산

feat_dim = size(O_train, 1);
W_out = (Y_target * O_train') *  inv(O_train * O_train' + alpha * eye(feat_dim));


% Lorenz attractor prediction

Y_pred = zeros(dim, test_len);

% train data의 마지막 부분(초기값)을 가져옵니다.
current_history = train_data(:, train_len - warmup : train_len);

for t = 1:test_len
    % 1. 현재 history buffer로 feature vector 생성
    vec_linear = [];
    for delay = 0:(k-1)
        % history의 끝에서부터 과거로 가져옴
        vec_linear = [vec_linear; current_history(:, end - delay*s)];
    end
    
    full_prod = vec_linear * vec_linear';
    vec_nonlinear = full_prod(triu(true(size(full_prod))));
    O_vec = [1; vec_linear; vec_nonlinear];
    
    % 2. 다음 스텝 예측 (행렬 곱)
    next_step = W_out * O_vec;
    Y_pred(:, t) = next_step;
    
    % 3. History 업데이트
    % 가장 오래된 데이터 버리고, 방금 예측한 값을 끝에 추가 (Closed-loop)
    current_history = [current_history(:, 2:end), next_step];
end


%% 시각화

figure('Position', [100, 100, 1000, 500]);
plot(test_data(1, :), 'k', 'LineWidth', 1.5); hold on;
plot(Y_pred(1, :), 'r--', 'LineWidth', 1.5);
title('Time Series Prediction (x-component)');
legend('Ground Truth', 'NG-RC Prediction');
xlabel('Time Steps'); ylabel('x'); grid on;