SINDy

1. 개요

SINDy Sparse Identification of Nonlinear Dynamics
입력	상태 시계열 x(t) (와 그 시간미분)
출력	항이 몇 개뿐인 미분방정식 dx/dt = f(x)
핵심 가정	f 가 후보 함수 라이브러리의 희소 선형결합
기본 풀이	순차 임계 최소제곱 (STLSQ)
제안	브런턴·프록터·쿠츠, 2016
구현	PySINDy

데이터를 신경망에 넣으면 잘 맞는 블랙박스가 나온다. SINDy 에 넣으면 교과서에 적을 수 있는 방정식이 나온다. 운이 좋으면.

SINDy(Sparse Identification of Nonlinear Dynamics, 비선형 동역학의 희소 식별)는 측정한 상태 시계열로부터, 미리 정해 둔 후보 함수(다항식·삼각함수 등) 라이브러리의 몇 안 되는 항만으로 이루어진 지배 미분방정식을 회귀로 찾아내는 방법이다. 2016년 스티븐 브런턴·조슈아 프록터·네이선 쿠츠가 PNAS 에 발표했다.¹

핵심 관찰은 단순하다. 물리 법칙은 대개 항이 적다. 로렌츠 방정식은 세 식에 항이 일곱 개뿐이고,² 나비에-스토크스도 몇 개 항의 합이다. 그러니 “가능한 항을 잔뜩 깔아 놓고, 그중 대부분의 계수가 0인 해”를 찾으면 모형 선택과 계수 추정이 한 번에 된다. 이는 라쏘·압축센싱의 희소 회귀를 동역학 방정식에 적용한 것이고, 2009년 슈미트·립슨의 기호 회귀(유전 프로그래밍으로 식 자체를 탐색)를 훨씬 싼 선형대수 문제로 바꿨다는 데 의의가 있다.

2. 회귀 문제로 쓰기

시각 $t_1, \dots, t_m$ 의 상태를 행으로 쌓은 $X \in \mathbb R^{m\times n}$ 과 그 시간미분 $\dot X$ 를 준비한다. 후보 함수 $p$ 개를 열로 늘어놓은 라이브러리

$$\Theta(X) = \big[\ \mathbf 1\ \ X\ \ X^{P_2}\ \ X^{P_3}\ \cdots\ \sin X\ \cdots\ \big] \in \mathbb R^{m\times p}$$

를 만든다($X^{P_2}$ 는 $x^2, xy, xz, \ldots$ 같은 2차 단항식 전부). 그러면 미지수는 계수 행렬 $\Xi \in \mathbb R^{p\times n}$ 하나이고

$$\dot X \approx \Theta(X)\,\Xi$$

를 $\Xi$ 의 각 열이 희소하도록 푼다. 열 $k$ 가 곧 $\dot x_k = \Theta(\mathbf x)^{\mathsf T}\boldsymbol\xi_k$ 의 계수다. 상태변수 하나당 독립된 선형 회귀이니 병렬화도 공짜다.

2.1. STLSQ — 순차 임계 최소제곱

원 논문의 기본 풀이는 놀랄 만큼 소박하다.

1. 최소제곱으로 $\Xi$ 를 푼다.
2. 절댓값이 임계값 $\lambda$ 보다 작은 계수를 0으로 만든다.
3. 남은 열만으로 최소제곱을 다시 푼다. 2–3을 지지 집합이 안 바뀔 때까지 반복한다.

$\ell_1$ 벌점을 쓰는 라쏘보다 단순하고, 살아남은 계수에 수축 편향이 없다는 장점이 있다. 장-셰퍼(2019)는 이 반복이 유한 번 안에 멈추고, 적절한 조건에서 $\ell_0$ 벌점 문제의 국소 최소점으로 간다는 것을 보였다. 실무에서는 $\lambda$ 와 함께 능형회귀 벌점을 살짝 섞는 변형(PySINDy 의 기본값)이 흔하다.

3. 검산 — 로렌츠 계 되찾기

$\sigma=10$, $\rho=28$, $\beta=8/3$ 의 로렌츠 계를 $(-8, 7, 27)$ 에서 4차 룽게-쿠타법으로 $\Delta t = 0.002$, 10000 스텝($t \in [0, 20]$) 적분했다. 미분은 중심차분으로 근사하고, 라이브러리는 $\{1, x, y, z, x^2, xy, xz, y^2, yz, z^2\}$ 열 개, $\lambda = 0.1$ 로 STLSQ 를 돌렸다.

잡음 없음. 30개 후보 계수 중 정확히 7개가 살아남았다.

$$\dot x = -9.9992\,x + 9.9992\,y,\quad \dot y = 27.992\,x - 0.9985\,y - 0.9998\,xz,\quad \dot z = -2.6664\,z + 0.9999\,xy$$

참값 $(-10, 10;\ 28, -1, -1;\ -8/3, 1)$ 과 넷째 자리까지 맞는다. 남은 오차는 중심차분의 $O(\Delta t^2)$ 이다.

임계값 없이(그냥 최소제곱). 상태에 표준편차 0.1의 잡음을 얹고 $\lambda = 0$ 으로 풀면, 큰 계수는 대충 맞지만 30개 전부가 0이 아니다. $\dot x$ 식에 상수 $-0.37$, $z$ 계수 $0.049$ 같은 가짜 항이 붙는다. 예측은 그럭저럭 되지만 “방정식을 발견했다”고는 말할 수 없다. 희소화 단계가 하는 일이 바로 이 쓰레기 청소다.

잡음을 키우면. 같은 잡음 0.1에 $\lambda=0.1$ 이면 일곱 항이 그대로 복원된다. 큰 계수는 0.2% 안쪽이고, 가장 작은 $-y$ 계수만 2% 남짓 벗어난다. 잡음을 1.0으로 키우자(상태 표준편차 약 8의 12% 정도) 사정이 달라졌다.

잡음	$\lambda$	결과
1.0	0.1	세 식 모두 가짜 항 다수. $\dot x$ 계수가 $3.49x + 1.75y$ 로 엉망
1.0	0.5	구조는 거의 복원($\dot x$ 2항·$\dot z$ 2항), 계수는 $-8.77x + 8.88y$ 처럼 약 12% 작게 편향, $\dot y$ 에서 $-y$ 항이 탈락

임계값을 올리면 가짜 항은 지워지지만 작은 진짜 항도 같이 지워진다. $\dot y$ 의 $-y$ 는 계수가 1이지만 $28x$ 옆에서는 상대적으로 작아서 희생됐다. 계수가 일제히 0쪽으로 줄어든 건 설명변수 $X$ 자체에 잡음이 낀 탓의 감쇠 편향(errors-in-variables)이다 — 총최소제곱이 따로 존재하는 이유와 같다.

4. 진짜 적은 미분이다

위 실험에서 상태 잡음 1.0은 중심차분을 거치면 미분 잡음 $1.0/(\sqrt2\,\Delta t) \approx 350$ 이 된다. 로렌츠 계의 $\dot x$ 크기가 수십~백 정도니 미분 신호 대 잡음비가 1보다 작다. 그런데도 구조가 반쯤 살아남은 건 만 개의 샘플이 평균을 내 줬기 때문이고, 샘플이 적었다면 그것도 없었다. 수치미분은 잡음을 $1/\Delta t$ 로 증폭하고, 그 잡음이 회귀의 좌변 전체에 들어간다.

그래서 SINDy 변형의 상당수는 미분을 어떻게 피하거나 다스리느냐의 이야기다.

• 평활 미분. 사비츠키-골레이 필터나 전변분 잡음제거형 정규화 미분으로 $\dot X$ 를 먼저 다듬는다. 가장 흔하지만 평활 폭이 또 하나의 손잡이가 된다.
• 적분형. 양변을 시간 적분해 $X(t) - X(t_0) = \int \Theta\,dt\ \Xi$ 로 바꾼다. 미분 대신 적분이 잡음을 평균 낸다.
• 약형식 SINDy(Messenger & Bortz, 2021). 방정식에 매끄러운 시험함수를 곱하고 부분적분해 미분을 시험함수 쪽으로 옮긴다. 유한요소법의 약형식과 같은 아이디어다. 잡음에 대한 강인성이 자릿수 단위로 좋아진다고 보고됐다.
• 앙상블·부트스트랩. 데이터를 재표집해 여러 번 SINDy 를 돌리고, 대부분의 표본에서 살아남는 항만 채택한다. 계수의 불확실성 추정이 덤으로 따라온다.

5. 확장 — PDE·제어·축소모형

• PDE-FIND(Rudy 등, 2017). 라이브러리에 $u_x, u_{xx}, uu_x, \ldots$ 같은 공간미분 항을 넣어 편미분방정식을 찾는다. 시공간 격자 데이터에서 버거스·KdV·나비에-스토크스(와도 형태)를 되찾은 것이 대표 예다. 공간 미분까지 수치로 해야 하니 잡음 문제는 더 심하다.
• SINDYc. 제어 입력 $\mathbf u$ 를 라이브러리에 넣어 $\dot{\mathbf x} = \mathbf f(\mathbf x, \mathbf u)$ 를 식별한다. 찾은 희소 모형은 모형예측제어(MPC)에 넣기 좋다 — 평가가 싸고 해석이 된다.
• 축소모형과 결합. 고차원 유동장을 적합직교분해 모드 몇 개로 줄이고, 그 모드 계수의 동역학을 SINDy 로 찾는다. 갤러킨 투영 없이 데이터만으로 축소차수모델을 만드는 길이다. 실린더 후류의 몇 모드 모형이 원 논문부터 시연 예제였다.
• 좌표 발견. 오토인코더로 잠재 좌표를 학습하면서 그 좌표에서의 동역학이 희소하도록 함께 훈련하는 방식도 있다. “어떤 좌표에서 보면 방정식이 단순해지는가”라는 질문은 쿠프만 연산자·동적 모드 분해가 선형성으로 답하려는 질문과 같은 뿌리다.

6. 한계

• 라이브러리 밖의 항은 절대 못 찾는다. $1/(1+x^2)$ 이 진짜 항인데 다항식만 깔았다면, SINDy 는 그걸 다항식 여러 개로 흉내 내며 희소성을 잃는다. 라이브러리 설계에 사람의 물리 지식이 들어간다.
• 열끼리 강한 공선성. 데이터가 끌개 위에 몰려 있으면 $x^2$ 과 $y^2$ 같은 열이 거의 평행해 어느 쪽을 남길지가 불안정해진다. 여러 초기조건·과도 구간을 섞어야 한다.
• $\lambda$ 선택. 위 실험처럼 결과가 $\lambda$ 에 민감하다. 교차검증이나 아카이케 정보기준으로 파레토 곡선(오차 대 항 개수)의 꺾임점을 고르는 게 국룰이다.
• 찾은 모형이 맞다는 보장은 없다. 짧은 구간에서 잘 맞는 희소 모형이 장시간 적분하면 발산하거나 다른 끌개로 갈 수 있다. 식별 뒤 반드시 적분해서 원 데이터와 비교해야 한다. 이건 모든 데이터 기반 모형의 검증 및 확인 숙제다.³

7. 관련 문서

• 라쏘 · 압축센싱 · 능형회귀 · 총최소제곱 · 교차검증
• 동적 모드 분해 · 쿠프만 연산자 · 확장 DMD · 적합직교분해 · 축소차수모델
• 로렌츠 방정식 · 시스템 식별 · 수치미분 · 전변분 잡음제거 · 기호 회귀
• 물리 정보 신경망 · 편미분방정식

8. Footnotes

1. 읽을 때는 보통 「신디」라고 한다. 약어를 사람 이름처럼 들리게 짓는 건 이 바닥 논문 마케팅의 국룰이다. ↩

2. 그 일곱 항으로 나비 모양 끌개와 카오스가 다 나온다. 항이 적다고 거동이 단순하다는 뜻은 아니다 — 그래서 희소 회귀가 먹힌다. ↩

3. “데이터로 방정식을 발견했다”는 논문 대부분은 답을 아는 계를 되찾은 것이다. 답을 모르는 계에서 찾은 식이 맞는지는 결국 실험이 정한다. 갈릴레이가 들으면 웃을 이야기다. ↩