[POD-1] 고전 적합직교분해 (classical POD)

적합직교분해(POD, proper orthogonal decomposition)은 본래 난류 유동(turbulent flow)에서 가장 에너지가 높은 모드를 추출하기 위해서 도입된 수학적인 기법이다.

 

 

\(\mathbf{q}(\mathbf{p},t)\)를 위치벡터 \(\mathbf{p}\)와 시간 \(t\)를 독립변수로 하는 벡터 필드 (예를 들면 유동장에서의 속도 벡터)라고 하자. 이 벡터의 시간 평균을 \(\bar{\mathbf{q}}(\mathbf{p})\)라고 하면 벡터 필드가 평균을 기준으로 변동하는 성분(unsteady component)은 다음과 같이 기저함수(basis function)의 선형 조합으로 나타낼 수 있다.

 

\[ \mathbf{q}(\mathbf{p},t)-\bar{\mathbf{q}}(\mathbf{p}) = \sum_{j=1}^\infty a_j \mathbf{\phi}_j (\mathbf{p},t) \tag{1} \]

 

여기서 \(\mathbf{\phi}_j (\mathbf{p},t)\)는 기저함수이고 \(a_j\)는 기저함수의 중요도를 나타내는 계수다. 이러한 표현 방식을 일반화된 푸리에 급수(generalized Fourier series)라고 한다.

 

 

여기서 한발 더 나아가 계수를 시간의 함수로, 그리고 기저함수를 위치만의 함수로 가정하면 벡터 필드를 시간과 공간을 분리하여 표현할 수도 있다.

 

\[ \mathbf{q}(\mathbf{p},t)-\bar{\mathbf{q}}(\mathbf{p}) = \sum_{j=1}^\infty a_j(t) \mathbf{\phi}_j (\mathbf{p}) \tag{2} \]

 

물론 이런 변수 분리 방식이 모든 문제에 다 적용되는 것은 아니지만, 만약 이러한 분리가 가능하다면 공간적인 특성만 갖는 기저함수를 추출하여 마치 좌표축처럼 어떤 벡터의 운동을 기저함수로 표현할 수 있을 것이다.

기저함수 \(\mathbf{\phi}_j (\mathbf{p})\)가 좌표축의 역할을 하고 \(a_j (t)\)가 축 성분의 역할을 하려면 기저함수가 일정한 영역 내에서 단위 직각함수(orthonormal function)이면 편리할 것이다. 단위 직각함수는 다음과 같은 성질을 갖는 함수를 말한다.

 

\[ \int_V \mathbf{\phi}_i^T (\mathbf{p}) \mathbf{\phi}_j (\mathbf{p}) \ dV = \delta_{ij} \tag{3} \]

 

여기서 \(V\)는 관심 영역이고 \(\delta_{ij}\)는 크로넥커 델타 함수로서 \(i=j\) 이면 \(1\)이고 \(i \ne j\) 이면 \(0\)인 값을 갖는다. 그러면 계수 \(a_j (t)\)는 다음과 같이 편리하게 계산할 수 있다.

 

\[ \int_V \mathbf{\phi}_i^T (\mathbf{p}) \sum_{j=1}^\infty a_j (t) \mathbf{\phi}_j (\mathbf{p}) \ dV = a_i (t) \]

 

벡터 필드가 연속적인 공간상에서 값을 갖는 함수이기 때문에 해석적인 해를 구해야 하나 이는 매우 어렵거나 불가능하고 유체나 구조 해석에서는 공간을 그리드로 분할하여 수치해석적인 방법으로 해를 계산한다.

 

 

관심 영역에 그리드(grid)를 짜고 식 (2)를 \(N\)개의 그리드 포인트에서 표현하면 다음과 같다.

 

\[ \begin{bmatrix} \mathbf{q}(\mathbf{p}_1,t) -\bar{\mathbf{q}}(\mathbf{p}_1 ) \\ \vdots \\ \mathbf{q}(\mathbf{p}_N,t) -\bar{\mathbf{q}}(\mathbf{p}_N ) \end{bmatrix} = \sum_{j=1}^\infty a_j (t) \begin{bmatrix} \mathbf{\phi}_j (\mathbf{p}_1) \\ \vdots \\ \mathbf{\phi}_j (\mathbf{p}_N) \end{bmatrix} \]

 

여기서 \(\mathbf{p}_i\)는 그리드 포인트 \(i\)의 위치벡터를 나타낸다. 위 식을 간단히 다음과 같이 쓰기로 하자.

 

\[ \mathbf{y}(t)= \sum_{j=1}^\infty a_j (t) \mathbf{\psi}_j \tag{4} \]

 

여기서 \(\mathbf{y}(t)\)는 그리드 포인트로 이산화된(discretized) 벡터 필드가 평균을 기준으로 변동하는 성분을 열벡터로 나타낸 것이며 \(\mathbf{\psi}_j\)는 그리드 포인트로 이산화된 기저함수를 열벡터로 나타낸 기저벡터로서 다음과 같이 주어진다.

 

\[ \mathbf{y}(t) = \begin{bmatrix} \mathbf{q}(\mathbf{p}_1,t) -\bar{\mathbf{q}}(\mathbf{p}_1 ) \\ \vdots \\ \mathbf{q}(\mathbf{p}_N,t) -\bar{\mathbf{q}}(\mathbf{p}_N ) \end{bmatrix} \ \in \mathbb{R}^n, \ \ \ \mathbf{\phi}_j = \begin{bmatrix} \mathbf{\phi}_j (\mathbf{p}_1) \\ \vdots \\ \mathbf{\phi}_j (\mathbf{p}_N) \end{bmatrix} \ \in \mathbb{R}^n \]

 

벡터 \(\mathbf{y}(t)\)의 차원 \(n\)은 그리드 포인트 갯수와 벡터 \(\mathbf{q}(\mathbf{p},t)\)의 차원 간의 곱으로 주어진다. 그리드를 촘촘하게 짤수록, 관심 영역의 범위가 넓을수록 \(n\)은 더 커질 것이다. 보통 그 값은 \(n=10^8 \sim 10^{10}\)에 달한다.

기저함수 \(\mathbf{\phi}_j\)는 단위 직각함수로서 식 (3)을 만족해야 하므로 그리드 포인트에서 계산한 기저함수를 이용하면 식 (3)을 다음과 같이 근사화 할 수 있다.

 

\[ \begin{align} \int_V \mathbf{\phi}_i^T (\mathbf{p}) \mathbf{\phi}_j (\mathbf{p}) \ dV &= \delta_{ij} \\ \\ & \approx \Delta V \left( \mathbf{\phi}_i^T (\mathbf{p}_1 ) \mathbf{\phi}_i (\mathbf{p}_1 ) + \cdots + \mathbf{\phi}_i^T (\mathbf{p}_N ) \mathbf{\phi}_i (\mathbf{p}_N ) \right) \\ \\ &= \Delta V \mathbf{\psi}_i^T \mathbf{\psi}_i \end{align} \]

 

여기서 \(\Delta V\)는 일정 간격으로 촘촘히 짜여 있는 그리드의 셀(cell) 부피이다. 위 식은 기저벡터 \(\mathbf{\psi}_j\)도 직각이 됨을 의미하나 그 크기는 셀의 부피에 영향을 받는다는 것을 알 수 있다. 따라서 이 영향을 배제하기 위해서 식 (4)에서 \(\mathbf{\psi}_j\)를 크기를 \(1\)로 고정시킨 직교 단위 벡터 \(\mathbf{w}_j\)로 변경한다.

 

\[ \begin{align} & \mathbf{y}(t) = \sum_{j=1}^\infty a_j (t) \mathbf{w}_j \tag{5} \\ \\ & \mathbf{w}_i^T \mathbf{w}_j= \delta_{ij} \end{align} \]

 

적합직교분해(POD)는 \(n\)차원 공간상에서 \(d\)개의 직교 단위 벡터 \(\mathbf{w}_1, \mathbf{w}_2, ..., \mathbf{w}_d\) 로 구성된 좌표축을 결정하는 문제로서 벡터 \(\mathbf{y}(t)\)를 그 좌표축으로 이루어진 부분 공간에 투사(projection)할 때 투사 오차가 최소가 되도록 좌표축을 계산한다. POD에서는 이 좌표축을 POD 모드(mode)라고 한다.

 

 

먼저 벡터 \(\mathbf{y}(t)\)를 \(d\)개의 좌표축으로 이루어진 부분 공간으로 근사하면 다음과 같다.

 

\[ \begin{align} \mathbf{y}(t) & \approx \sum_{j=1}^d a_j(t) \mathbf{w}_j \\ \\ &= \begin{bmatrix} \mathbf{w}_1 & \mathbf{w}_2 & \cdots & \mathbf{w}_d \end{bmatrix} \begin{bmatrix} a_1(t) \\ a_2(t) \\ \vdots \\ a_d(t) \end{bmatrix} \\ \\ &= W \mathbf{a} (t) \end{align} \]

 

여기서 \(a_j(t)\)는 \(\mathbf{y}(t)\)의 \(\mathbf{w}_j\)축 성분으로서 \(a_j(t)=\mathbf{w}_j^T \mathbf{y}(t)\) 로 주어지며,

 

\[ \mathbf{a}(t)= W^T \mathbf{y}(t) \]

 

이다. \( W \in \mathbb{R}^{n \times d}\)는 서로 직각인 단위 벡터(orthonormal vectors)로 이루어진 행렬로서 다음 식을 만족한다.

 

\[ W^T W=I_d \]

 

투사 오차는 다음과 같으므로

 

\[ \begin{align} \mathbf{e}_{proj} (t) &= \mathbf{y}(t)- \sum_{j=1}^d a_j (t) \mathbf{w}_j \\ \\ &= \mathbf{y}(t)-WW^T \mathbf{y}(t) \end{align} \]

 

시간 영역 \(t=[0, T]\) 에서 POD문제를 다음과 같이 최적화 문제로 정식화 할 수 있다.

 

\[ \begin{align} & \arg\min_W⁡ J = \int_0^T \lVert \mathbf{y}(t)- WW^T \mathbf{y} (t) \rVert_2^2 \ dt \\ \\ & \mbox{subject to } \ \ W^T W=I_d \end{align} \]

 

 

 

함수 \(J\)를 좀더 전개해보자.

 

\[ \begin{align} J &= \int_0^T \lVert \mathbf{y} (t)- WW^T \mathbf{y} (t) \rVert_2^2 \ dt \\ \\ &= \int_0^T \left( \mathbf{y}^T (t) \mathbf{y}(t)- \mathbf{y}^T (t) WW^T \mathbf{y} (t) \right) \ dt \\ \\ &= \int_0^T \mathbf{y}^T (t) \mathbf{y}(t) \ dt - trace \left[ \int_0^T W^T \mathbf{y} (t) \mathbf{y}^T (t) W \ dt \right] \\ \\ &= \int_0^T \mathbf{y}^T (t) \mathbf{y}(t) \ dt - trace \left[ W^T \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) W \right] \end{align} \]

 

그런데 여기서 \(\mathbf{y}(t)\)는 \(W\)에 종속적이지 않으므로, 최소화 문제는 다음과 같이 최대화 문제로 바뀌게 된다.

 

\[ \begin{align} \arg\min_W & \left( \int_0^T \mathbf{y}^T (t) \mathbf{y}(t) \ dt - trace \left[ W^T \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) W \right] \right) \\ \\ &= \arg\min_W \left( - trace \left[ W^T \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) W \right] \right) \\ \\ &= \arg\max_W \left( trace \left[ W^T \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) W \right] \right) \end{align} \]

 

정리하면 POD문제는 다음 최적화 문제로 바꿀 수 있다.

 

\[ \begin{align} & \arg\max_W \left( trace \left[ W^T \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) W \right] \right) \\ \\ & \mbox{subject to } \ \ W^T W=I_d \end{align} \]

 

라그랑지 곱수 \(\lambda_{ij}\)를 이용하여 제약조건이 없는 최적화 문제로 바꾸면 다음과 같다.

 

\[ \max_W L = trace \left[ W^T \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) W \right] + \sum_{i=1}^d \sum_{j=1}^d \lambda_{ij} ( \delta_{ij}-\mathbf{w}_i^T \mathbf{w}_j ) \]

 

최적화의 필요조건에 의하면 다음 미분식을 만족해야 한다.

 

\[ \begin{align} & \frac{\partial L}{\partial \mathbf{w}_i} = 2 \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) \mathbf{w}_i - 2 \lambda_{ii} \mathbf{w}_i - \sum_{j=1,i \ne j}^d \lambda_{ij} \mathbf{w}_j =0 \\ \\ & \frac{\partial L}{\partial \lambda_{ij}} = \delta_{ij} - \mathbf{w}_i^T \mathbf{w}_j =0 \end{align} \]

 

따라서 \(\mathbf{w}_i\)와 라그랑지 곱수 \(\lambda_{ij}\)는 다음 식을 만족한다.

 

\[ \begin{align} & \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) \mathbf{w}_i= \lambda_{ii} \mathbf{w}_i \\ \\ & \lambda_{ij}=0, \ \ \mbox{if} \ i \ne j \end{align} \]

 

즉 라그랑지 곱수 \(\lambda_{ii}\)와 \(\mathbf{w}_i\)는 각각 행렬 \(\left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right)\)의 고유값과 고유벡터가 됨을 알 수 있다.

행렬 \(\left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right)\)는 준정정(semi-positive) 행렬이므로 고유값은 0보다 크거나 같은 실수값이며 고유벡터도 서로 직각인 실수 벡터이다. 따라서 행렬을 다음과 같이 고유값 분해(eigen decomposition)로 표현할 수 있다.

 

\[ \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) = P \Lambda P^T \]

 

여기서 고유값을 크기가 큰 수에서 작은 수로 정렬하면,

 

\[ \begin{align} & \Lambda = diag( \lambda_1, ..., \lambda_d, ..., \lambda_n ), \ \ \ \lambda_1 \ge \cdots \ge \lambda_n \ge 0 \\ \\ & P= \begin{bmatrix} \mathbf{v}_1 & \cdots & \mathbf{v}_d & \cdots & \mathbf{v}_n \end{bmatrix} \end{align} \]

 

POD 모드는 \(\mathbf{v}_i, \ i=1, ..., d\) 를 선택하면 된다.

 

\[ \begin{align} W &= \begin{bmatrix} \mathbf{w}_1 & \mathbf{w}_2 & \cdots & \mathbf{w}_d \end{bmatrix} \\ \\ &= \begin{bmatrix} \mathbf{v}_1 & \mathbf{v}_2 & \cdots & \mathbf{v}_d \end{bmatrix} = P_d \ \in \mathbb{R}^{n \times d} \end{align} \]

 

그러면 투사 오차는 다음과 같이 계산되므로

 

\[ \begin{align} J &= \int_0^T \lVert \mathbf{y} (t)- WW^T \mathbf{y}(t) \rVert_2^2 \ dt \\ \\ &= trace \left[ \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right] - trace \left[ W^T \left( \int_0^T \mathbf{y} (t) \mathbf{y}^T (t) \ dt \right) W \right]\\ \\ &= trace \left[ P \Lambda P^T \right] - trace \left[ P_d^T P \Lambda P^T P_d \right] \\ \\ &= \sum_{i=1}^n \lambda_i -\sum_{i=1}^d \lambda_i \\ \\ &= \sum_{i=d+1}^n \lambda_i \end{align} \]

 

에너지(L2 norm)가 가장 높은 순으로 POD 모드를 추출하게 되는 것이다.

고유값을 이용하면 벡터 필드 데이터의 변동을 나타내는데 필요한 POD 모드 수를 결정할 수 있다. 일반적으로 다음 식을 이용하여 모드 수 \(d\)를 결정한다.

 

\[ \frac{ \sum_{i=1}^d \lambda_i}{ \sum_{i=1}^n \lambda_i } \approx 1 \]

 

고전 POD를 정리하면 다음과 같다.

1. 관심 영역에 그리드(grid)를 구축하고 \(N\)개의 그리드 포인트에서 벡터 필드의 변동 \(\mathbf{q}(\mathbf{p},t)-\bar{\mathbf{q}}(\mathbf{p})\) 를 계산한다.

 

\[ \mathbf{y}(t)= \begin{bmatrix} \mathbf{q}(\mathbf{p}_1,t) - \bar{\mathbf{q}}(\mathbf{p}_1) \\ \vdots \\ \mathbf{q}(\mathbf{p}_N,t) - \bar{\mathbf{q}}(\mathbf{p}_N) \end{bmatrix} \ \in \mathbb{R}^n \]

 

2. 시간 영역 \(t=[0, T]\)에서 \(d\)개의 POD 모드를 다음과 같이 계산한다.

 

\[ \left( \int_0^T \mathbf{y}(t) \mathbf{y}^T (t) \ dt \right) \mathbf{w}_i=\lambda_i \mathbf{w}_i, \ \ i=1, ..., d \]

 

3. POD 모드의 성분은 다음과 같이 계산한다.

 

\[ a_i (t) = \mathbf{w}_i^T \mathbf{y} (t) \]

 

4. 고차원 벡터 필드는 다음과 같이 복원시킬 수 있다.

 

\[ \begin{bmatrix} \hat{\mathbf{q}}(\mathbf{p}_1,t) \\ \vdots \\ \hat{\mathbf{q}}(\mathbf{p}_N,t) \end{bmatrix} = \begin{bmatrix} \bar{\mathbf{q}}(\mathbf{p}_1) \\ \vdots \\ \bar{\mathbf{q}}(\mathbf{p}_N) \end{bmatrix} + \sum_{j=1}^d a_j(t) \mathbf{w}_j \]

 

고전(classical) POD는 공간을 이산화시켰지만 시간은 연속적이다. 따라서 행렬 \( \left( \int_0^T \mathbf{y}(t) \mathbf{y}^T (t) dt \right) \)를 해석적으로 계산해야 하는 어려움이 있다(실제로는 수치해석으로 계산해야 한다). 행렬을 계산한 다음에는 이를 바탕으로 고유값과 고유벡터를 계산해야 하는데 행렬의 사이즈가 \(n \times n\) (여기서 \(n=10^8 \sim 10^{10}\) )에 달하기 때문에 실제적으로 고전 POD를 적용하기는 거의 불가능하다.

고전 POD의 단점을 극복하기 위한 방안으로 스냅샷(snapshot) POD가 개발되었는데 이에 대해서는 다음에 알아보기로 하자.