[수치해석] Operation cost / Condition number(Stiffness)

normal-engineer.tistory.com/83

이 글과 이어지는 내용입니다.

 

 

1. Solution 계산에 필요한 operation cost

 

앞에서 modeling - transform to multiple ODE - transform to algebraic equation 까지 진행한 다음에는

역행렬을 계산해야합니다.

이 글에서는 matrix의 형태에 따라 역행렬을 구하는 operation cost가 달라진다는 것을 보이려고 합니다.

 

 

참고

$\sum_{k=1}^{n}1=n$

$\sum_{k=1}^{n}k=\frac{n(n+1)}{2}$

$\sum_{k=1}^{n}k^{2}=\frac{n(n+1)(2n+1)}{6}$

$\sum_{k=1}^{n}k^{3}=\left ( \frac{n(n+1)}{2}\right )^{2}$

 

 

1) Diagonal matrix

$\begin{bmatrix}7&0&0\\0&5&0\\0&0&1\end{bmatrix}\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\end{bmatrix}=\begin{bmatrix}8\\3\\0\end{bmatrix}$

 

$n\times n $diagonal matrix

$\text{Operation cost }~\theta (n^{1})$

 

division을 n번 나눴으므로 total operation count n번.

 

 

2) Upper Triangular matrix

$\begin{bmatrix}3&5&6\\0&2&5\\0&0&7\end{bmatrix} \begin{bmatrix}x_{1}\\x_{2}\\x_{3}\end{bmatrix}=\begin{bmatrix}6\\3\\8\end{bmatrix}$

 

back substitution을 사용합니다.

 

$7x_{3}=8 \rightarrow$ 1 division

$2x_{2}+5x_{3}=3 \rightarrow$ 1 substitution(multiplication), 1 substraction, 1 division 

$3x_{1}+5x_{2}+6x_{3}=8 \rightarrow$ 2 multiplication, 2 substraction, 1 division

 

매 row마다 division이 1번씩 존재하고, row가 올라갈 때마다 multiplication과 substraction(일종의 addition)이 1회씩 늘어납니다.

$n\times n\text{ matrix }$의 operation cost : $\sum_{k=1}^{n}2(k-1)+1=n(n+1)-n=n^{2}$

$\text{Operation cost }~\theta (n^{2})$

 

 

3) General (Sparse) matrix $\rightarrow$ Gauss elimination

$A\underline{x} = \underline{b}$

 

$\begin{bmatrix}a_{11} & a_{12} & \cdots & a_{1n}\\ a_{21} & \ddots &  & a_{2n}\\ \vdots &  &  & \vdots\\ a_{n1} &  &  & a_{nn}\end{bmatrix} \begin{bmatrix}x_{1}\\ x_{2}\\ \vdots \\ x_{n}\end{bmatrix}=\begin{bmatrix}b_{1}\\b_{2}\\ \vdots \\b_{n}\end{bmatrix}$

 

위의 1),2)와 달리 모든 항에 임의의 숫자가 있는 경우, 목표는 general matrix를 upper triangular matrix로 만들어서 2)번을 적용하는 것입니다.

 

알고리즘

(1). if $a_{11}=0 \overset{\text{yes}}{\longrightarrow} \text{row change}$ : 만약 $a_{11}=0$일 경우 $a_{i1}\neq0$ 인 row와 바꿔줍니다. (permutation)

$a_{11}\neq 0$이면 $l_{2}=\frac{a_{21}}{a_{11}}$으로 1 division을 시행합니다.

 

(2). first row에 $l_{2}$를 곱해줍니다. 총 n번의 multiplication($a_{11}$을 제외한 나머지 element와 source element인 $b_{1}$까지)

 

(3). second row에 $l_{2}$를 곱한 first row를 빼줍니다. 총 n번의 substraction. 어차피 첫번째 항은 사라질 것이라 보기 때문에 굳이 곱셈과 뺄셈을 진행하지 않아도 됩니다.

 

(4). 마찬가지로 third, fourth, fifth... row의 첫번째 항을 제거하기 위해 1),2),3) 과정을 계속 거칩니다.

(5). 첫번째 항이 모두 제거되었다면, second row를 이용해 third row부터 final row까지 두번째 항을 제거합니다.

(6). (n-1)th row까지 진행하면 끝. (nth row는 위의 operation을 통해 마지막 element밖에 없음)

 

따라서 총 들어가는 operation count

Total division : $\sum_{k=1}^{n-1}(n-k) \sim \theta (n^{2})$

Total multiplication : $\sum_{k=1}^{n-1} (k+1)k \sim \theta(n^{3})$

Total addition : $\sum_{k=1}^{n-1} (k+1)k \sim \theta(n^{3})$

 

$\text{Total Operation cost }~\theta (n^{3})$

 

이렇게 만들어놓고 다시 upper triangular matrix에 대해 solution을 구해야하므로 $n^{2}$만큼의 operation이 더 필요합니다.

order가 3인 경우는 지나치게 operation cost가 크므로 아예 사용을 배제한다고 합니다.

 

물론 아주 정확한 방법이기는 하지만 시간이 너무 오래 걸린다면 효율적이지 못합니다. 따라서 특정 method를 통해 근사화된 답을 구함으로써 order를 2, 더 줄여서 1까지로 만들 수 있습니다.

 

 

4) Tridiagonal matrix

이번에는 tridiagonal matrix의 operation cost에 대해서 알아보겠습니다.

 

$\begin{pmatrix}b_{1} & c_{1} & 0 & 0 & \cdots & 0 \\ a_{2} & b_{2} & c_{2} & 0 & \cdots & 0\\ 0 & a_{3} & b_{3} & c_{3} & \cdots & 0\\ \vdots &  &  &  & \ddots & \vdots\\ 0 &  &  &  & a_{n} & b_{n} \end{pmatrix}$

 

(1). Eliminate $a_{2}$

 

1 division : $l=a_{2}/b_{1}$

2 multiplication : 첫번째 row에 l을 곱합니다. $c_{1}\rightarrow c_{1}\times l \text{  and source element multiplied by }l$

2 addition(substraction) : $b_{2}-c_{1}\times l \text{  and source element substraction}$

 

$a_{2}$는 당연히 제거되므로 계산에 포함시키지 않았습니다.

$a_{2}$를 제거하면 first row의 과정은 끝나고, second row로 가서 위 과정을 진행합니다. 따라서 모든 row에 대해서 똑같이 1D, 2M, 2A operation이 필요합니다.

 

(2). Make upper triangular matrix : total operation은 (n-1)D+2(n-1)M+2(n-1)A $\sim \theta (n^{1})$

 

(1)번 과정을 거치면 tridiagonal matrix는 bidiagonal matrix가 됩니다. (아래의 element들은 임의의 숫자이므로, 위의 triangular matrix에 나온 항과는 별개입니다.)

(참고 : bidiagonal matrix는 tridiagonal matrix에서 a항, 또는 c항이 전부 0인 경우)

 

$\text{bidiagonal matrix : }\begin{pmatrix}b_{1} & c_{1} & 0 & 0 & \cdots & 0 \\ 0 & b_{2} & c_{2} & 0 & \cdots & 0\\ 0 & 0 & b_{3} & c_{3} & \cdots & 0\\ \vdots &  &  &  & \ddots & \vdots\\ 0 &  &  &  & 0 & b_{n}\end{pmatrix}$

 

 

이와 같은 bidiagonal matrix(이자 upper triangular matrix)로 만든 다음에 backward substitution을 합니다.

위의 예시는 $\theta(n^{2})$이었지만 지금은 다릅니다.

 

(3). Backward substitution : 각 row 당 1 division, 1 multiplication, 1 substraction 총 3번 operation이 필요하다. (마지막 row 제외하고)

total operation of backward substitution  $\sum_{k=2}^{n} 3 + 1 = 3(n-1)+1 \sim \theta(n^{1})$

 

따라서 (1),(2),(3)을 고려하면 tridiagonal matrix의 operation cost는 다음과 같습니다.

$\text{Total Operation cost }~\theta (n^{1})$

 

 

tridiagonal matrix일 때 훨씬 계산이 덜 필요합니다.

 

수치해석에서는 이렇게 operation cost가 적은 matrix를 얻기 위해 노력합니다.

CFD 난류 해석에서 10만개에서 몇 억개의 n이 필요하다는 것을 생각했을 때 order가 중요함을 쉽게 납득할 수 있습니다.

 

 

 

질문 : permutation은 operaction cost에 고려하는가?
답변 : 단순한 row change는 포인터만 바꿔주면 되는 문제라서 cost에 거의 영향을 미치지 않는다. 하지만 logical operation(위 matrix가 invertible인지 아닌지 밝혀내는 과정)은 cost가 매우 크다. 또한 permutation을 수행해서 미리 row를 정렬하는 과정을 해놓을 수도 있는데 이것을 pre-conditioning이라고 한다. (cost가 많이 들 수 있음에도) pre-conditioning을 하는 이유는 A matrix가 잘 정립되면 그 후부터는 b(source term)만 바꿔가면서 solution을 구하기 때문이다.

---

2. Condition number

 

 

마지막으로 condition number(stiffness)에 대해서 이야기하면서 introduction을 마칩니다.

 

수치해석에서 고려하는 parameter 중 하나가 condition number이자 stiffness입니다.

 

이 stiffness는 solution의 안정성을 의미합니다.

일반적으로 $Ax=b$에서 A matrix를 구하고 나면 b라는 source term을 바꿔가면서 solution x를 구합니다.

 

그런데 아래의 예시를 보면,

 

$\left\{\begin{matrix}x_{1}+x_{2}=2\\ x_{1}+1.0001x_{2}=2\end{matrix}\right.\Rightarrow \begin{matrix}x_{1}=2\\ x_{2}=0\end{matrix}$

 

$\left\{\begin{matrix}x_{1}+x_{2}=2\\ x_{1}+1.0001x_{2}=2.0001\end{matrix}\right.\Rightarrow \begin{matrix}x_{1}=1\\ x_{2}=1\end{matrix}$

 

작은 source term의 변화에도 $x_{1},x_{2}$값이 크게 변한 것을 확인할 수 있습니다. matrix에 따라 작은 변화에도 크게 solution에 영향을 미친다면 문제가 발생할 수 있습니다.

 

컴퓨터의 유효숫자는 32bit(single precision)와 64bit(double precision)에 따라서 나타낼 수 있는 maximum 유효숫자 자리가 다릅니다. 따라서 어떤 컴퓨터에서는 2.0000000과 2.0000001을 같은 숫자로 저장할 수 있고 위와 같이 작은 변화에도 민감하다면 부정확한 solution을 얻게 되는 것입니다.

 

 

그래서 b의 미소변화량에 대한 x의 미소변화량을 구해보려고 합니다.

 

$Ax = b$일 때

$A(x+\delta x)=b+\delta b$ : 이 식에서 $\delta b$에 비해 $\delta x$가 너무 크면 stiffness가 크다고 볼 수 있습니다. (less stiff is better)

 

$Ax-b+A\delta x = \delta b$

$A\delta x = \delta b$ $(\because Ax=b)$

$\delta x = A^{-1}\delta b$

 

 

Cauchy-schwarz inquality

(1). $\left \| \delta x \right \|=\left \| A^{-1}\delta b \right \|\leq \left \| A^{-1} \right \| \left \| \delta b \right \|$

(2). $\left \| b \right \| \leq \left \| A \right \| \left \| x \right \|$

 

$\Rightarrow \frac{\left \| \delta x \right \|}{\left \| x \right \|} \leq \left \| A^{-1} \right \|\left \| A \right \|\frac{\left \| \delta b \right \|}{\left \| b \right \|}$

 

위 부등식을 보면 $\frac{\left \| \delta x \right \|}{\left \| x \right \|}$는 x의 변화율이고 $\frac{\left \| \delta b \right \|}{\left \| b \right \|}$는 b의 변화율입니다.

 

두 식을 연결해주는 것이

$\left \| A^{-1} \right \|\left \| A \right \|=K(A) \text{ : A matrix의 condition number}$

이 condition number가 큰 matrix는 stiff합니다.

 

주로 $K(A) \approx \frac{\left | \lambda_{max} \right |}{\left | \lambda_{min} \right |}$

두 eigenvalue maximum, minimum 값의 비율로 볼 수 있습니다.

 

 

---

다음에는 interpolation에 대해 이야기합니다.