[Complexity] 리아프노프 지수 - 로지스틱 맵

 시계열 데이터의 Complexity를 정량화 하는 방법에 리아프노프 지수(Lyapunov exponents)가 있다. Lyapunov 지수는 간단하면서도 카오스적 성질을 지닌 Logistic map의 예제를 통해 그 의미를 가늠해볼 수 있다.

 시계열 데이터에서 Complexity란 단순히 진동하거나 증가/감소 하는 성질이 아닌 복잡하고도 예측 불가능한 성질이다. 이러한 성질은 카오스라고도 불리우며, 위상공간에서 기이한 끌개(strange attractor)의 운동을 보이는데, '초기 조건에의 민감성'을 갖는다.

 '초기 조건에의 민감성'은 나비효과로 널리 알려 있는데, 아주 작은 초기 상태의 차이가 시간이 지날 수록 점점 거대한 차이로 나타나는 특성으로, 계의 예측을 어렵게 혹은 불가능하게 만드는 카오스 상태의 본질이다.

 리아프노프 지수는 초기 조건에의 민감성을 말 그대로 식으로 나타내어 정량화 한 값이다. 위상 공간에서 작은 초기값의 차이가 시간이 지남에 따라 지수적으로 증가할 때, 지수 함수의 계수를 나타낸다.

 

 x는 위상공간에서  두 궤적의 거리를 나타내며, 궤적의 거리는 지수 λ가 0보다 작을 때는 수렴하고, 0일 때는 유지, 0보다 클 때는 지수적으로 증가하여 예측하기 힘든 성질을 야기한다. 이러한 위상공간을 지닌 계는 '카오스'라고 불리우며, 1차원 위상공간의 간단한 예(로지스틱 맵)을 통해 확인할 수 있다.

로지스틱 사상 (Logistic map)

---

 로지스틱 맵은 현 상태가 다음 상태를 결정하는 이산 시간 동역학계로, 다음과 같은 식으로 나타낼 수 있다.

x는 상태이며, r에 따라 x의 거동이 변화한다.

 X의 초기 값 X_0을 0.1로 두고, step별 변화를 나타낸 그림은 다음과 같다.

r = 2 에서는 x는 한 값으로 수속, 3에서는 주기적 진동, 4에서는 카오스적 행태를 보인다.

 이번에는 초기 값 0.1은 파란색, 0.1+10^-8은 빨간색으로 변화를 살펴본 결과이다.

r = 4, 즉, 카오스 상태에서는 초기의 아주 작은 차이 10의 -8승으로 인해 두 궤적에 큰 변화가 나타난다.

소스코드)

import matplotlib.pyplot as plt


def logistic_map_state(r, step, x0):
    x = x0
    x_list = [x]
    for n in range(step):
        x_n = r * x * (1 - x)
        x = x_n
        x_list += [x]
    return x_list

fig = plt.figure(figsize=(8, 3))
subplot = fig.add_subplot(1, 3, 1)
r = 2
x_list = logistic_map_state(r=r, step=50, x0=0.1)
x_list_d = logistic_map_state(r=r, step=50, x0=0.1+1e-8)
plt.plot(x_list)
plt.plot(x_list_d, c='red')
plt.ylim([-0.1, 1.1])
plt.title('r=2')
plt.ylabel('x')
plt.xlabel('step')
subplot = fig.add_subplot(1, 3, 2)
r = 3
x_list = logistic_map_state(r=r, step=50, x0=0.1)
x_list_d = logistic_map_state(r=r, step=50, x0=0.1+1e-8)
plt.plot(x_list)
plt.plot(x_list_d, c='red')
plt.ylim([-0.1, 1.1])
plt.title('r=3')
plt.xlabel('step')
subplot = fig.add_subplot(1, 3, 3)
r = 4
x_list = logistic_map_state(r=r, step=50, x0=0.1)
x_list_d = logistic_map_state(r=r, step=50, x0=0.1+1e-8)
plt.plot(x_list)
plt.plot(x_list_d, c='red')
plt.ylim([-0.1, 1.1])
plt.title('r=4')
plt.xlabel('step')
plt.tight_layout()

 동역학계의 상태를 표시하는 방법에는 Bifurcation diagram이라는 방법이 있다. step에 따른 status의 변화가 수속하는지, 진동하는지, 카오스적 행태를 보이는지 확인할 수 있는 직관적인 방법인데, Logistic map에서는 계수 r에 따른 status의 변화를 다음과 같이 나타낼 수 있다.

r이 2.5에서 약 3 부근 사이에서는 x는 하나의 값으로 수속한다. r이 3 부근에서 x의 값에 분기가 일어나며, x는 두 값 사이에서 진동한다. x의 상태는 r~3.5 부근에서 또 한번의 분기하여 4개의 값 사이에서 진동하고, 조금 더 높은 r에서는 8개의 값 사이에서 진동하다가 더 이상 주기성을 알 수 없는 카오스적 상태에 돌입하는 것을 직관적으로 확인할 수 있다.

소스코드)

import numpy as np
import matplotlib.pyplot as plt


def logistic(r, x):
    return r * x * (1 - x)

# 관측할 r과 각 r의 상태 x 생성
n = 10000
r = np.linspace(2.5, 4.0, n)
x = 0.1 * np.ones(n)

# step 수 및 diagram에 나타낼 state 개수(last) 설정
step = 1000
last = 100

fig = plt.figure()
for n in range(step):
    x = logistic(r, x)
    # 각 r에 대한 마지막 100개의 상태 x 플롯
    if n >= (step - last):
        plt.plot(r, x, marker=',', ls='', c='black', alpha=0.2)
plt.xlim(2.5, 4)
plt.title("Bifurcation diagram")
plt.xlabel('r')
plt.ylabel('x')
plt.tight_layout()

---

리아프노프 지수 (Lyapunov exponent)

---

 다시, 리아프노프 지수로 돌아와보자. 리아프노프 지수는 이러한 초기조건 민감성을 더 정량적으로 확인할 수 있는 방법이다. 우리는 제일 처음 제기하였던 식에서 λ를 구해야 한다.

 로지스틱 사상에서 변이는 step별로 나타나므로, t는 n으로 대체할 수 있다. 이 때 λ는 다음과 같이 구할 수 있다.

 다음, 로지스틱 사상에서 초기 값의 차이가 시간이 지날 수록 어떻게 성장하는지 살펴보자. 로지스틱 사상은 다음식으로 나타낼 수 있다.

로지스틱 사상을 함수 f로 두고, 스탭 0에서 스탭 1로의 변화를 살펴보자. 

x_0은 x_1로, x_0+dx_0은 x_1+dx_1로 맵핑된다. 이 때, dx_1을 윗식으로부터 전개하면 다음과 같이 나타낼 수 있다.

 이 과정을 스탭 n으로 확장하면, 다음과 같다.

 리아프노프 지수에서는 두 궤적의 차이의 '크기'에 관심이 있기 때문에 궤적 간 거리에 대한 식은 다음과 같다.

 이 식을 위에서 구한 λ를 구하는 식에 대입하면, 다음과 같다.

 이 때, n을 무한대(아주 큰 수)로 두고 구한 λ를 maximal Lyapunov exponent라고 부르는데, 변덕스러운 fluctuation을 잡고, 큰 흐름 속에서의 궤적의 변화를 정량화하기 위함이다.

 이를 로지스틱 사상에 대입하면 다음과 같고, 매 스탭 구하여 더해준 후 평균 값을 취하면 된다.

이렇게 구한 각 r 별 Lyapunov exponent를 Bifurcation diagram과 같이 나타내면 다음과 같다.

 수렴구간, 주기구간에서 리아프노프 지수는 0 이하의 값을, 카오스 구간에서 0 이상의 값을 갖는 것을 알 수 있다. 이렇게, 계가 지닌 카오스적 성질은 Lyapunov exponent를 통해 가늠해 보는 것이 가능하다.

소스코드)

import numpy as np
import matplotlib.pyplot as plt


def logistic(r, x):
    return r * x * (1 - x)


# 관측할 r과 각 r의 상태 x 생성
n = 10000
r = np.linspace(2.5, 4.0, n)
x = 0.1 * np.ones(n)

# step 수 및 diagram에 나타낼 state 개수(last) 설정
step = 1000
last = 100

# r에 따른 Lyapunov exponents 설정
lyapunov = np.zeros(n)

fig = plt.figure()
subplot = fig.add_subplot(2, 1, 1)
for n in range(step):
    x = logistic(r, x)
    # Lyapunov exponents 계산
    lyapunov += np.log(abs(r - 2 * r * x))
    # 각 r에 대한 마지막 100개의 상태 x 플롯
    if n >= (step - last):
        plt.plot(r, x, marker=',', ls='', c='black', alpha=0.2)
lyapunov = lyapunov / step
plt.xlim(2.5, 4)
plt.title("Bifurcation diagram")
plt.xlabel('r')
plt.ylabel('x')
subplot = fig.add_subplot(2, 1, 2)
plt.plot(r, lyapunov, marker=',', ls='', c='black', alpha=0.2)
plt.axhline(y=0, c='red')
plt.xlim(2.5, 4)
plt.title("Lyapunov exponent")
plt.tight_layout()

 추가로 Maximal Lyapunov exponent(MLE)가 무엇을 뜻하는지, 구간 별 Lyapunov 지수와 함께 살펴보자.

r = 2.5일 경우, r=3일 경우, r=4일 경우의 각 스탭 별 f'(x_n)값을 검은 점으로, 그 평균 값인 MLE를 빨간 선으로 나타내었다.

 값이 수렴하는 영역인 r=2.6에서는 깔끔하게 Lyapunov지수가 수렴하지만, 주기영역인 r=3.3과 카오스영역의 r=4에서는 그렇지 아니하다. 고로, 위의 값들의 평균을 취함으로써, fluctuation을 제거한 전체적인 trend를 살펴보는 방법으로써 MLE가 도입된 것을 알 수 있다.

 Logistic map의 경우에서는 f(x)의 미분 값을 구할 수 있는 경우였지만, 실제 Data에서는 model을 알 수 없는 경우가 많으므로, 미분 값을 대입하는 것이 아닌 시간 변화와 궤적 간 거리의 변화를 이용해서 구간 별 λ를 구하고, 안정된 영역에서 이 값의 평균을 구하는 방식으로 리아프노프 지수를 구한다.

---

참고자료

• https://ocw.hokudai.ac.jp/wp-content/uploads/2016/01/ChaosFractal-2011-Note-04.pdf
• http://hatano-lab.iis.u-tokyo.ac.jp/thesis/soturon2000/soturon_fujibuchi.pdf
• http://chaosbook.org/chapters/Lyapunov.pdf
• http://fy.chalmers.se/~f99krgu/dynsys/DynSysLecture10.pdf
• http://www.fulviofrisone.com/attachments/article/412/Chaos%20From%20Simple%20Models%20to%20Complex%20Systems.pdf
• http://blog.daum.net/neuronics/7567593
• https://ipython-books.github.io/121-plotting-the-bifurcation-diagram-of-a-chaotic-dynamical-system/