선형시스템 알고리즘의 연산 복잡도 비교

선형대수학에서 가장 기본이 되는 과제 중 하나는 선형시스템이라고 불리는 식, 즉 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 효율적으로 풀어내는 것이다. 여기서 $\mathbf{A}$는 $n\times n$ 크기의 정사각행렬, $\mathbf{x}$는 $n\times 1$ 크기의 미지수 벡터, 그리고 $\mathbf{b}$는 $n\times 1$ 크기의 주어진 벡터로 가정한다. 고전적인 관점에서 선형시스템 풀이는 행렬을 직접 다루는 “직접해(direct method)”와 필요한 연산을 반복적으로 수행하여 근사해를 구하는 “반복해(iterative method)”로 크게 나눌 수 있다. 이 두 방법은 목적은 같지만 알고리즘의 구조와 연산 복잡도가 서로 다르며, 부동소수점 연산이 얼마나 수행되는지에 따라 계산 시간과 메모리 사용량, 그리고 결과 오차 등이 달라진다.

크기가 $n$인 선형시스템을 해결하기 위해서는 일반적으로 $O(n^p)$ 형태의 연산이 필요하다고 알려져 있다. $p$의 값은 알고리즘의 차수(order)를 말하며, 주어진 알고리즘이 어떤 방식으로 문제의 크기에 비례하는지(또는 그 이상의 비례식을 갖는지)를 나타낸다. 본 장에서는 가장 많이 쓰이는 직접해 알고리즘들과 대표적인 반복해 알고리즘들을 예로 들어, 그 연산 복잡도를 비교하고 상황별로 어떻게 달리 적용되는지 자세히 살펴본다.

직접해 알고리즘의 개요

직접해 알고리즘은 이론적으로 유한 단계 안에(부동소수점 연산에서의 반올림 오차는 무시한다고 가정할 때) 정확한 해를 구하는 것이 가능하다는 장점을 갖는다. 대표적인 예로 가우스 소거법(Gaussian Elimination)을 들 수 있으며, 이를 변형한 $LU$ 분해나 $QR$ 분해, 그리고 대칭행렬에 특화된 조셉스키(Cholesky) 분해 등이 있다. 연산 복잡도를 구체적으로 살펴보기 위해 가장 기초가 되는 가우스 소거법부터 살펴보자.

가우스 소거법은 다음 단계로 크게 요약할 수 있다. 첫째, 적절한 피벗을 선택해가며 $\mathbf{A}$를 삼각화(triangularization)한다. 둘째, 삼각화된 행렬에서 후진대입(back-substitution) 또는 전진대입(forward-substitution)을 통해 $\mathbf{x}$를 구한다. 여기서 삼각화 과정에서 생기는 연산량이 전체 알고리즘 연산의 대부분을 차지한다. 특정 $n$에 대하여 가우스 소거법의 연산량을 부동소수점 연산 횟수로 엄밀히 계산해 보면 $O(n^3)$이 된다.

예를 들어 $LU$ 분해 역시 가우스 소거법을 살짝 변형한 것으로, $\mathbf{A} = \mathbf{L}\mathbf{U}$ 꼴로 분해한 뒤 $\mathbf{L}\mathbf{y} = \mathbf{b}$, $\mathbf{U}\mathbf{x} = \mathbf{y}$의 두 단계로 선형시스템을 푸는 방법이다. 이때 $LU$ 분해 자체는 여전히 $O(n^3)$ 정도의 연산을 요구하며, 분해가 완료된 이후에는 $\mathbf{b}$가 달라져도 $\mathbf{x}$를 구하는 데 필요한 연산량은 $O(n^2)$ 정도로 작아진다. 이러한 특징 때문에 $\mathbf{A}$는 변하지 않고 $\mathbf{b}$만 여러 번 바뀌는 시스템에서 매우 유용하다.

직접해 알고리즘을 구현할 때 피벗팅(pivoting)은 필수적으로 고려되어야 한다. $\mathbf{A}$의 특정 행이나 열의 원소가 매우 작아서 수치적으로 불안정하면, 그 위치를 바꾸어 연산 안정성을 높이는 것이 일반적이다. 이를 부분 피벗팅(partial pivoting) 또는 전체 피벗팅(full pivoting)이라 부르는데, 경우에 따라 연산 횟수는 약간 증가할 수 있지만, 수치 안정성 측면에서 매우 중요한 역할을 한다.

가우스 소거법의 연산량 상세

가우스 소거법의 연산 과정을 단계별로 구체적으로 살펴보면 다음과 같은 형태로 연산 횟수를 근사할 수 있다. 먼저 첫 번째 열을 피벗으로 삼아 모든 행에 대해 소거 연산을 수행하면, 약 $n-1$개의 행을 대상으로 각각 $n-1$번 정도의 곱셈 및 뺄셈이 일어난다. 그 다음 두 번째 열을 피벗으로 삼아 소거할 때에는 약 $n-2$개의 행에 대해 $n-2$번의 곱셈 및 뺄셈이 발생한다. 이런 식으로 반복되는 과정을 모두 합산하면 결국

n(n1)(n2)6n36\frac{n(n-1)(n-2)}{6}\approx \frac{n^3}{6}

정도의 부동소수점 연산이 필요함을 알 수 있다. 물론 계수가 $\frac{1}{6}$인 것은 정확한 연산 횟수를 계산할 때 등장하는 것이고, 빅오 표기법 관점에서는 $O(n^3)$로 나타낸다.

후진대입이나 전진대입 과정은 각각 $O(n^2)$ 정도의 연산만 소요되므로, 전체 알고리즘의 지배적인 항은 삼각화 과정에 의해 결정된다. 따라서 가우스 소거법의 전체 연산 복잡도는 $O(n^3)$이다. 이러한 결론은 $LU$ 분해나 $QR$ 분해 또는 조셉스키 분해 등에서도 동일하게 유지되거나, 적어도 $n^3$ 차수의 연산량을 요구한다.

행렬 특성에 따른 분해 알고리즘의 변화

만일 $\mathbf{A}$가 특별히 희소(sparse) 구조를 갖거나 대칭 성질을 갖는다면, 추가적인 최적화를 통해 연산 복잡도를 낮출 수 있다. 예를 들어 대칭이고 양의 정부호(positive definite)인 $\mathbf{A}$에 대해서는 조셉스키 분해가 가능한데, 이 경우에도 $O(n^3)$ 복잡도를 갖지만, 분해 과정에서 필요한 연산 횟수의 상수 계수가 일반적인 $LU$ 분해보다 작아진다. 또한 대각 행렬이나 삼각 행렬과 같이 더욱 간단한 구조를 가진다면, $O(n)$ 또는 $O(n^2)$ 단위에서 계산이 끝나는 경우도 있다.

하지만 일반적인 경우, 특별한 구조를 갖지 않는 임의의 밀집(dense) 행렬 $\mathbf{A}$에 대해서는 $O(n^3)$ 복잡도를 갖는 직접해 알고리즘 이상으로 빨리(예: $O(n^2)$나 $O(n\log n)$) 해결할 수 있는 대표적인 방법은 존재하지 않는다. 따라서 문제의 크기가 매우 클 때에는 메모리나 시간 측면에서 직접해 알고리즘의 비용이 상당히 커지며, 이런 상황에서 반복해 알고리즘이 유리할 수 있다.

반복해 알고리즘의 소개

반복해 알고리즘은 초기 추정값 $\mathbf{x}^{(0)}$에서 시작하여, 특정 규칙에 따라 $\mathbf{x}^{(k+1)}$을 갱신해 나가면서 원하는 수준의 오차 이내로 해를 근사하는 방법이다. 예를 들어 야코비(Jacobi) 방법, 가우스-자이델(Gauss-Seidel) 방법, 그리고 이보다 더 발전된 사후조건자(preconditioner)를 사용하는 방법들이 많다. 반복해 알고리즘은 한 번의 반복마다 드는 연산량은 직접해보다 훨씬 작을 수 있으나, 필요한 반복 횟수에 따라 총 연산량이 달라진다. 조건수가 큰 행렬을 다룰 경우에는 수렴 속도가 매우 느려져 오히려 $O(n^3)$보다 큰 비용이 들기도 한다.

이 장의 다음 부분에서는 대표적인 반복해 알고리즘들의 수렴 특성과 연산 복잡도를 살펴보고, 대규모 시스템에서 이 방법들을 어떻게 활용해야 하는지 구체적으로 살펴본다. 특히 대규모 희소 행렬에서 반복해 알고리즘이 어떻게 효율을 낼 수 있는지, 그리고 사후조건자를 어떠한 방식으로 설계해야 하는지 등을 심도 있게 다룰 것이다.

반복해 알고리즘의 기본 원리

반복해 알고리즘(iterative method)은 초기 추정값 $\mathbf{x}^{(0)}$에서 시작하여, 특정 규칙에 따라 해를 조금씩 갱신해 가며 목표 오차 범위 안에서 해를 근사하는 방식이다. 수렴이 빠르다면 직접해에 비해 연산량이 대폭 절감될 수 있지만, 행렬의 조건수(condition number)가 좋지 않으면 수렴에 많은 반복이 필요하여 전체 연산이 늘어날 수 있다. 이 장에서는 반복해 알고리즘이 어떤 식으로 설계되는지, 그리고 그 연산 복잡도를 어떻게 파악하는지 살펴본다.

제일 간단한 케이스인 야코비(Jacobi) 방법과 가우스-자이델(Gauss-Seidel) 방법을 통해 반복해 알고리즘의 기본 아이디어를 소개하고, 이어서 좀 더 고차원적인 기법인 SOR(Successive Over-Relaxation) 방법, 대규모 선형시스템에 적합한 사후조건자(preconditioner) 기법, 그리고 Krylov 서브스페이스(Krylov subspace) 기법 중 대표적인 CG(Conjugate Gradient), GMRES(Generalized Minimal Residual) 등에 대해서도 개략적으로 짚어볼 것이다.

야코비(Jacobi) 방법과 가우스-자이델(Gauss-Seidel) 방법

야코비 방법은 $\mathbf{A} = \mathbf{D} + \mathbf{L} + \mathbf{U}$ 형태로 분할했을 때, 주대각 성분만 따로 떼어낸 $\mathbf{D}$를 이용하여 해를 갱신하는 전형적인 알고리즘이다. 구체적으로 $i$번째 성분 $x_i^{(k+1)}$는

xi(k+1)=1aii(bijiaijxj(k))x_i^{(k+1)} = \frac{1}{a_{ii}} \Bigl( b_i - \sum_{j \neq i} a_{ij} x_j^{(k)} \Bigr)

와 같이 업데이트된다. 모든 성분을 동시에 업데이트한다는 점이 야코비 방법의 특징으로, 병렬화가 용이하다는 장점이 있다. 하지만 수렴 속도는 비교적 느린 편이며, 행렬 $\mathbf{A}$의 대각 성분이 충분히 커서(즉 대각우세(diagonally dominant) 등) 안정된 수렴을 보장할 때만 사용된다.

가우스-자이델 방법은 야코비 방법과 달리, 한 개의 성분을 갱신하면서 곧바로 그 값을 다음 성분의 업데이트에 사용한다. 즉

xi(k+1)=1aii(bij=1i1aijxj(k+1)j=i+1naijxj(k))x_i^{(k+1)} = \frac{1}{a_{ii}} \Bigl( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)} \Bigr)

와 같은 식으로 갱신되어, 전체 순회 한 번을 마치면 $\mathbf{x}^{(k+1)}$가 완성된다. 야코비에 비해 일반적으로 더 빠른 수렴 성능을 보이지만, 병렬화 측면에서는 덜 유리할 수 있다.

야코비나 가우스-자이델 방법 모두, 한 번의 반복에서 필요한 연산량은 대략 $O(n^2)$ 정도에 해당한다. 왜냐하면 각 성분마다 약 $n$ 정도의 합이나 곱이 이루어지기 때문이다. 따라서 이들이 실제로 효율적인지는 “필요한 반복 횟수”에 크게 의존한다. 반복 횟수가 $K$라면 총 연산량은 $O(K n^2)$가 된다. 이 $K$가 매트릭스의 조건수나 행렬 구조에 의존하여 달라지므로, 실제로 얼마나 빨리 수렴하는지는 사전 지식이나 추가적인 분석이 필요하다.

SOR와 사후조건 기법

가우스-자이델 방법은 야코비 방법에 비해 빠른 수렴을 보이긴 하지만, 여전히 문제의 특성(예: 행렬의 스펙트럼 반경 등)에 따라 매우 느리게 수렴할 수도 있다. 이를 개선하기 위한 대표적인 변형 기법 중 하나가 SOR(Successive Over-Relaxation) 방법이다. SOR에서는 가우스-자이델의 업데이트 식에 과보정(over-relaxation) 파라미터 $\omega$를 곱해, 다음과 같이 업데이트를 진행한다.

xi(k+1)=xi(k)+ω(x~i(k+1)xi(k))x_i^{(k+1)} = x_i^{(k)} + \omega \Bigl(\tilde{x}_i^{(k+1)} - x_i^{(k)}\Bigr)

여기서 $\tilde{x}_i^{(k+1)}$는 가우스-자이델 방법에 의해 계산된 값이고, $0 < \omega < 2$의 범위에서 적절한 $\omega$를 선택함으로써 수렴 속도를 개선할 수 있다. 이론적으로 최적의 $\omega$는 행렬 $\mathbf{A}$의 고유값 분포에 의해 결정되며, 일반적인 문제에서는 실험적 혹은 경험적으로 결정하는 경우가 많다.

사후조건자(preconditioner)를 사용하는 기법은 반복 과정에서 행렬 $\mathbf{A}$ 대신에 $\mathbf{M}^{-1}\mathbf{A}$ (또는 그 유사 형태)를 사용하는 방법으로, $\mathbf{A}$의 조건수를 인공적으로 개선하여 수렴을 빠르게 만드는 전략이다. 이때 $\mathbf{M}$을 사후조건자라고 부르는데, $\mathbf{M}$은 쉽게 역을 구할 수 있으면서도 $\mathbf{M}^{-1}\mathbf{A}$가 수렴 특성에서 큰 이득을 얻을 수 있도록 설계된다. 사후조건자는 간단한 대각행렬을 쓰는 방법부터, 대규모 문제에 적합한 복잡한 기법까지 매우 다양하다.

사후조건 기법의 효과를 정량화하려면 행렬 $\mathbf{M}^{-1}\mathbf{A}$의 스펙트럼 반경이나 조건수 등을 살펴봐야 하며, 이를 계산하는 일 자체도 만만치 않다. 하지만 적절한 사후조건자를 잘 선택했을 때, 반복 횟수를 대폭 줄여주는 효과가 있기 때문에, 대규모 희소 행렬 시스템에서는 사실상 사후조건 기법이 필수적으로 사용된다. 다만 $\mathbf{M}$을 구하고 $\mathbf{M}^{-1}$ 곱셈을 수행하는 비용도 무시할 수 없으므로, 전체 복잡도를 감안하여 균형점을 찾아야 한다.

Krylov 서브스페이스 방법의 연산 복잡도

직접해가 $O(n^3)$ 복잡도를 보이는 반면, 반복해 방법 중 일부는 큰 차원 문제에서 효과적인 수렴을 보여 실무적으로도 자주 사용된다. 그중에서도 Krylov 서브스페이스 방법들이 널리 알려져 있다. Krylov 서브스페이스란, 초기 잔차(residual) $\mathbf{r}^{(0)} = \mathbf{b} - \mathbf{A}\mathbf{x}^{(0)}$에서 시작하여, $\mathbf{A}^k\mathbf{r}^{(0)}$ ($k=0, 1, 2,\dots$)에 의해 생성되는 부분공간을 의미한다. 이 부분공간 내에서 특정 조건을 만족하는 $\mathbf{x}$를 찾는 방식이 Krylov 방법의 핵심이다.

대표적인 Krylov 방법 중 하나는 CG(Conjugate Gradient) 방법이다. CG는 $\mathbf{A}$가 대칭 양의 정부호(symmetric positive definite)라는 가정 아래, 적절히 정의된 에너지 노름(energy norm)을 최소화하는 해를 Krylov 공간에서 효율적으로 찾는 알고리즘이다. 한 번 반복할 때 필요한 주요 연산은 $\mathbf{A}\mathbf{p}$ 형태의 곱셈(여기서 $\mathbf{p}$는 탐색 방향 벡터)과 몇 개의 벡터 연산이므로, 대략 $O(n^2)$ 정도로 평가할 수 있다. 반복 횟수는 행렬 $\mathbf{A}$의 조건수 $\kappa$에 의해 크게 좌우되는데, 균등(equivalently well-conditioned)하다면 적은 횟수로 수렴한다.

GMRES(Generalized Minimal Residual)는 일반적인(대칭이 아닐 수도 있는) 선형시스템에 적용 가능한 Krylov 서브스페이스 방법이다. 매 반복마다 $\mathbf{A}$ 곱셈 연산과 잔차를 최소화하기 위한 직교화(orthogonalization) 과정을 거쳐야 하므로, 구현 방식에 따라서는 추가 메모리와 연산이 필요하다. 반복 횟수 역시 문제 크기가 커짐에 따라 늘어날 수 있으나, 사후조건자와 결합하거나 재시작 기법(restarting)을 사용하면 실용적인 계산 시간을 확보할 수 있다.

이와 같은 Krylov 방법의 공통된 특징은 “매 반복마다 $O(n^2)$ 정도의 연산”을 필요로 하며, “반복 횟수 $K$가 행렬의 조건수나 스펙트럼 분포에 따라 크게 변한다”는 것이다. 따라서 최적화된 사후조건자를 적용하여 $K$를 매우 작게 만들 수 있다면, 결과적으로 $O(K n^2)$의 연산으로 $n^3$에 준하는 직접해보다 효율적인 해결책을 얻을 수 있다. 반면 사후조건자가 부실하거나, 문제 구조가 삐뚤어져 있어 $K$가 $O(n)$ 이상으로 커져버린다면, 결국 $O(n^3)$에 도달하거나 이를 초과할 수 있다는 점에 주의해야 한다.

중규모 vs. 대규모 문제에서의 적용

중규모 문제(예: $n$이 몇 천에서 몇 만 정도)에서는 아직 $O(n^3)$가 불가능하지 않을 수도 있다. 예를 들어, 현대적인 선형대수 라이브러리를 사용하면 병렬 연산과 고성능 연산장치(GPU, 다중코어 CPU 등)의 도움으로, $n$이 수천에서 수만 정도인 문제까지도 직접해 알고리즘(예: $LU$ 분해, $QR$ 분해 등)로 해결할 수 있다. 반면 그 이상의 초대형 문제(예: $n$이 수십만, 수백만 이상)에서는 메모리와 계산 시간이 너무 커서 직접해 알고리즘을 사실상 적용하기 어렵다.

이런 초대형 문제에서 반복해 알고리즘은 메모리 접근 패턴이 간단하고, 행렬-벡터 곱셈에 특화된 희소 행렬 처리 기법을 활용하면 극적으로 빨라질 수 있다. 또한 일부 행렬 요소만 접근하게 되어 캐시 이용(또는 분산 메모리 구조)에서 효율을 보이는 경우가 많다. 반면 반복 횟수를 최소화하기 위해서는 반드시 사후조건자나 적절한 초기값 선택이 필요하므로, 문제의 특성에 맞는 알고리즘 설계가 필수적이다.

병렬화와 대규모 문제 해결

현대 컴퓨팅 환경은 다중코어 CPU와 GPU, 클러스터, 클라우드 기반의 대규모 분산 메모리 시스템을 활용할 수 있으므로, 선형시스템 알고리즘에도 병렬화(parallelization)가 매우 중요한 요소가 되었다. 본질적으로 $O(n^3)$의 연산량을 갖는 직접해 알고리즘이라고 하더라도, 큰 $n$에서 병렬 연산이 효율적으로 이루어진다면, 실질적인 계산 시간은 상당히 단축될 수 있다. 예를 들어 LU 분해의 삼각화 과정은 피벗팅이 없다면(또는 블록 단위 피벗을 활용한다면) 병렬화가 비교적 용이하며, 대규모 고성능 선형대수 라이브러리들이 이를 잘 활용하고 있다. GPU 가속을 사용하는 경우, 벡터-행렬 연산이나 BLAS 레벨-3 연산(예: 행렬 곱셈 등)을 효율적으로 처리하는 방식으로 $n$이 수만, 수십만을 넘어서는 문제까지도 직접해가 적용되는 사례가 늘어나고 있다.

반복해 알고리즘도 마찬가지로 병렬화 측면에서 다양한 개선이 이루어지고 있다. 야코비(Jacobi) 방법은 각 성분을 독립적으로 갱신한다는 특성상 병렬화가 가장 쉽지만, 수렴 속도가 느리다는 약점이 있다. 가우스-자이델(Gauss-Seidel)이나 SOR, Krylov 서브스페이스 계열의 알고리즘은 순차적으로 의존하는 연산이 포함되므로 병렬 효율이 떨어질 수 있지만, 대규모 희소 행렬에 대해선 행렬-벡터 곱셈을 병렬화하여 전체 성능을 끌어올리는 전략이 가능하다. 또한 사후조건(preconditioner)을 블록 단위로 적용하거나, 영역분할(domain decomposition) 기법을 사용하여 부분문제를 동시에 풀고 결과를 합치는 병렬 알고리즘도 활발히 연구 및 활용된다.

영역분할(domain decomposition) 및 블록 기반 기법

초대형 문제를 풀 때 자주 쓰이는 방법 중 하나가 영역분할(domain decomposition) 기법이다. 예를 들어 물리공간에서 정의되는 편미분방정식을 유한요소법이나 유한차분법으로 이산화하면 거대한 선형시스템 $\mathbf{A}\mathbf{x} = \mathbf{b}$가 형성되는데, 이를 여러 하위 영역으로 분할한 뒤, 각 영역에 해당하는 부분행렬(submatrix)을 독립적으로 풀고 경계(interface)에서만 적절한 조건을 만족하도록 연결하는 방식이다. 각 영역 내에서는 직접해나 반복해를 어떤 식으로든 사용할 수 있으며, 영역 간의 경계 정보를 갱신하면서 전체 문제의 해로 수렴한다.

블록 기반 기법도 유사한 사고방식을 취한다. 행렬 $\mathbf{A}$를 적당한 크기의 블록으로 나누고, 각 블록을 효율적으로 처리한 뒤 그 결과를 합성하는 방식이다. 이를테면 블록 Jacobi나 블록 Gauss-Seidel, 블록 SOR 등을 정의할 수 있는데, 각 블록 내에서의 연산이 독립적이므로 병렬화가 가능하다. 한 블록은 작은 차원의 부분문제를 이루며, 이 부분문제를 직접해로 처리하거나 간단한 사후조건자에 기반한 반복해를 적용할 수도 있다.

이런 영역분할과 블록 기법 모두, 이론적 연산 복잡도는 문제 전체 크기 $n$에 대해 여전히 $O(n^3)$ 혹은 반복 횟수가 포함된 $O(Kn^2)$ 형태가 될 수 있다. 하지만 실제로 하드웨어 자원을 최대한 효율적으로 활용함으로써, 계산 시간을 대폭 줄일 수 있는 장점이 있다. 따라서 대규모 병렬 컴퓨팅 환경에서 효율적 스케일링을 달성하기 위해서는, 단일 알고리즘에 집착하기보다는 분산 메모리 구조와 캐시 구조를 고려한 알고리즘 설계가 필수적이다.

반올림 오차와 수치 안정성

선형대수 알고리즘의 연산 복잡도뿐만 아니라, 수치 안정성(numerical stability)에 대한 고려도 매우 중요하다. 부동소수점 연산을 무한 정밀도로 할 수 없다면, 모든 계산 단계에서 반올림 오차(round-off error)가 누적된다. 직접해를 예로 들면, 가우스 소거법에서 피벗팅(pivoting)을 적절히 수행하지 않으면, 작은 피벗 요소에 의해 큰 오차가 증폭될 가능성이 있다. 이러한 문제는 $LU$ 분해, $QR$ 분해, 조셉스키(Cholesky) 분해 등에서도 각각의 방식에 맞는 조건이 맞아야 안정적으로 수행된다.

반복해 방법에서도 오차 축적 문제는 무시할 수 없다. 반복 횟수가 늘어날수록 작은 반올림 오차가 누적되어 해가 어긋나는 현상이 있을 수 있다. 그러나 대부분의 현대적 Krylov 방법은 잔차(residual) 기반의 갱신을 수행하므로, 직교화 과정을 주기적으로 하거나(예: GMRES) 스칼라 곱 구간을 주의 깊게 처리하는 방식으로 안정성을 개선한다. 또한 사후조건 기법에서 $\mathbf{M}$이나 $\mathbf{M}^{-1}$이 특정 구조를 갖는다면, 이 부분에서도 반올림 오차를 최소화하기 위한 전략이 필요하다.

수치 안정성을 정량적으로 평가하려면, 보통 $\kappa(\mathbf{A})$로 표기하는 행렬 $\mathbf{A}$의 조건수를 통해 접근한다. $\kappa(\mathbf{A})$가 크면(즉 조건이 나쁘면), 작은 반올림 오차가 해에 크게 반영될 수 있으므로, 알고리즘 선택뿐 아니라 사후조건자 설계나 피벗팅 전략을 꼼꼼히 살펴야 한다. 동시에 부동소수점 포맷(배정밀도, 단정밀도 등)에 따라서도 실제 계산에서의 오류 양상이 달라질 수 있음을 유념해야 한다.

---과 연결 고찰

본 장에서는 선형시스템 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 해결하기 위한 여러 가지 알고리즘의 연산 복잡도를 중심으로 살펴보았다. 일반적인 밀집 행렬에 대한 직접해는 $O(n^3)$ 복잡도를 갖고, 반복해는 $O(Kn^2)$ 정도로 볼 수 있으나, 실제로는 조건수와 구조적 특성, 그리고 병렬화 여부에 따라 더 복잡한 양상이 나타난다. 크기가 중간 정도인 문제에서는 여전히 직접해가 선호되는 경우가 많지만, 매우 큰 스케일의 문제에서는 반복해 방법과 사후조건자가 사실상 필수적이다. 또한 하드웨어 자원을 적극 활용할 수 있는 병렬화 기법이나 영역분할 기법도 현대 선형대수계산에서 핵심적인 위치를 차지하고 있다.

앞으로 남은 부분에서는 실제 예제 문제에 대해 여러 알고리즘을 구현해 보고, 연산 복잡도와 계산 시간, 그리고 수치 안정성 측면에서 어떤 차이가 발생하는지 살펴볼 것이다. 또한 희소(sparse) 행렬이 포함된 대표 문제들의 사례 연구를 통해, 알고리즘 선택이 결과와 자원 효율에 얼마나 큰 영향을 미치는지 구체적으로 제시할 것이다.

다중격자(Multigrid) 기법

대규모 선형시스템을 풀 때 대표적으로 고려되는 강력한 방법 중 하나가 다중격자(multigrid) 기법이다. 이는 편미분방정식(PDE)을 이산화하여 생기는 매우 큰 스케일의 희소행렬 시스템에 특히 잘 맞는데, 특정 격자(mesh) 해상도에서 발생하는 오차를 보다 큰 스케일(더 성긴 격자)에서 효율적으로 보정할 수 있다는 아이디어에 기반한다.

다중격자 기법에서 전체 격자는 여러 단계로 나뉘는데, 이를 “격자 계층(mesh hierarchy)”라고 부른다. 가장 섬세한 격자를 ‘최상위 격자’라 하고, 그보다 공간 해상도가 낮은 격자(더 적은 수의 자유도)를 ‘하위 격자’로 차례차례 내려간다. 예를 들어 2차원 문제에서 $h$ 간격을 갖는 격자보다 요소가 훨씬 적은 $2h, 4h$ 간격의 격자들을 순차적으로 구성할 수 있다. 이 과정을 통해 문제 해상도가 작은 격자 공간에서는 상대적으로 큰 스케일의 잔차(residual)를 빠르게 제거할 수 있고, 다시 해상도가 높은 격자에서는 세밀한 오차를 정밀하게 줄이는 방식으로 전체 문제를 빠르게 수렴시킨다.

알고리즘 개념 흐름

다중격자 방법을 구성하는 핵심 루틴은 보통 “스무딩(smoothing)”, “제약(restriction)”, 그리고 “보간(prolongation)”으로 요약된다.

스무딩 단계는, 현재 해상도의 격자에서 잔차를 최대한 줄이기 위해 가우스-자이델(Gauss-Seidel)이나 SOR와 같은 간단한 반복해를 여러 번 적용한다. 이 과정을 통해 고주파(high-frequency) 성분의 오차가 급격히 완화되며, 저주파(low-frequency) 성분이 주로 남게 된다.

제약 단계는, 스무딩 후 남은 잔차(저주파 성분이 주도)를 더 성긴 격자로 “이동”하여, 저해상도 문제에서 이 오차를 보정할 기회를 얻는 것이다. 행렬 $\mathbf{A}$도 그에 맞춰 “응축(coarse grid operator)” 형태로 축소되며, 이는 주로 격자의 구성 방식을 바탕으로 수학적으로 정의된다.

보간 단계는, 하위 격자에서 구해진 보정값을 다시 원래(또는 상위)의 해상도 격자로 옮겨주는 과정이다. 저해상도에서 빠르게 줄어든 저주파 오차 정보를 상위 격자의 해에 반영함으로써, 고해상도 문제에서 전체적인 오차가 일시에 큰 폭으로 개선될 수 있다.

이후 상위 격자에서 다시 스무딩을 수행하고, 필요하다면 더욱 하위 격자로 내려가거나 올라오는 과정을 반복한다. 이 다단계 사이클을 V-사이클, W-사이클, F-사이클 등 다양한 패턴으로 구성할 수 있는데, 문제 특성과 계산 자원에 따라 알맞은 구성과 사이클 횟수를 설정한다.

연산 복잡도 분석

다중격자 기법의 이상적인 이론치에 따르면, 2차원(또는 3차원)에서 주어진 PDE로부터 유도된 희소 행렬 $\mathbf{A}$가 적절한 부조건(엘립틱 문제, 규칙적인 격자 구조, 매끄러운 해 등)을 만족할 경우, 문제 크기 $n$에 대해 $O(n)$ 혹은 $O(n\log n)$ 정도의 연산으로 수렴해내는 것이 가능하다고 알려져 있다. 이는 “하위 격자에서의 연산”이 빠르고, 전체 격자에 걸친 잔차의 저주파 성분과 고주파 성분을 단계적으로 제거하여 반복 횟수 $K$가 상대적으로 거의 상수(또는 매우 완만하게 증가)로 유지된다는 특징 때문이다.

보다 구체적으로, 각 격자 레벨마다 수행되는 연산량은 해당 격자에서의 자유도 수 $n_\ell$에 비례한다. 만일 전체 자유도 $n$이 여러 격자 레벨에 걸쳐 비슷한 비율로 나누어진다면, 한 번의 다중격자 사이클(V-사이클 등)을 수행하는 데 드는 총 연산량을 대략

n+n1++n1n_\ell + n_{\ell-1} + \dots + n_1

형태로 합산하게 된다. 격자를 2배씩 줄일 수 있다고 가정하면, $n + \frac{n}{4} + \frac{n}{16} + \dots$ 식으로 기하급수적 감소가 일어나므로, 최종 합은 상수배의 $n$에 수렴한다. 여기에 사이클 횟수가 문제에 따라 상수 혹은 로그 스케일로 증가하므로, 실질적으로 $O(n)$ 또는 $O(n \log n)$의 복잡도로 해결이 가능해진다.

물론 실제 문제에서는 격자 구조가 불규칙하거나, 경계 조건이 복잡하거나, PDE가 엘립틱 형태가 아닌 혼합형 편미분방정식이거나, 심지어 비선형성이 존재할 수도 있다. 이 경우 다중격자 기법이 이론적인 $O(n)$ 또는 $O(n\log n)$ 성능을 보장하지 못할 수도 있지만, 여전히 수많은 실용적 상황에서 매우 빠르고 우수한 성능을 발휘한다.

병렬성과 결합된 응용

다중격자 방법은 자연스럽게 병렬화가 잘 되는 측면이 많다. 각 격자 레벨에서 스무딩 연산은 가우스-자이델이나 SOR와 같은 순차 의존도가 있는 방법을 그대로 사용하기보다는, 레드-블랙(적-흑) 가우스-자이델 같은 병렬 변형을 쓰거나, 야코비 방법 계열을 병렬로 활용할 수 있다. 또한 대규모 클러스터나 GPU 환경에서는 제약·보간 과정에서도 영역분할(domain decomposition) 기법과 결합하여 각 부분영역에서 독립적으로 다중격자를 적용한 뒤, 경계만 상호 교환하는 전략을 펼칠 수도 있다.

다중격자의 성공적인 활용을 위해서는, 격자 생성부터 제약·보간 연산자 설계, 그리고 스무딩 기법 선택까지 통합적으로 고려해야 한다. 이 때문에 수치해석과 계산과학 분야에서 다중격자 기법은 많은 이론적 연구와 함께 실제 엔지니어링 문제(예: 유체역학, 구조해석, 전자기장 해석 등)에서 표준적인 솔버(solver)로 자리매김해 있다.

다중격자 구현 스케치 (예: Python)

다음은 다중격자 알고리즘의 골격을 단순화하여 표현한 코드 예시(Python)이다. 실제 문제에서는 격자 자료구조와 $\mathbf{A}$, 제약·보간 연산자를 체계적으로 구성해야 하며, 스무딩 연산 또한 병렬화된 방식을 사용할 수 있다.

위 코드는 다중격자 V-사이클 알고리즘의 개념적 뼈대를 보여주며, 실제로는 문제 특성에 따라 여러 격자 레벨을 계층적으로 준비하고, 각 레벨마다 스무딩 파라미터와 연산자를 달리 설정하며, 적절한 반복 종료 조건을 설정해야 한다. 또한 $A$, $R$, $P$가 모두 희소 행렬 구조로 관리될 때에는, NumPy가 아닌 별도의 희소행렬 라이브러리(예: SciPy의 sparse 모듈)를 사용해 연산을 가볍게 해주어야 실질적인 효율을 얻을 수 있다.

다중격자 기법은 이처럼 PDE 기반 선형시스템에 특화되어 있으며, 때로는 비정형 격자(unstructured grid)나 적응형 격자(adaptive mesh refinement) 구조에서도 성공적으로 적용될 수 있다. 실제 산업 및 연구 현장에서는 Krylov 방법과 다중격자를 결합한 유연한 솔버가 많이 쓰이는데, 대규모 병렬 환경에서 고성능 계산을 해야 할 때 매우 중요한 도구로 인정받는다.

응용과 실제 벤치마크

실제 현장에서 어떤 선형시스템 솔버를 채택할지 결정하는 과정은 단순히 이론적 연산 복잡도를 비교하는 것만으로 끝나지 않는다. 문제의 구조(밀집행렬 또는 희소행렬, 대칭성 여부, 양의 정부호 여부 등), 문제 규모($n$), 사용 가능한 메모리, 그리고 하드웨어 특성(GPU 가속, 멀티코어, 분산 메모리 시스템 등)을 모두 고려해, 궁극적으로 “실제 계산 시간”과 “허용 가능한 오차 범위”를 종합적으로 검토해야 한다.

연구 및 실무에선 특정 문제를 두고 다양한 솔버를 벤치마크(benchmark)하여 비교한다. 예를 들어, 유한요소법으로 이산화된 구조해석 문제나 열전달 문제, 전자기장 해석 문제 등을 놓고, $n$의 크기를 여러 단계로 늘리면서 다음을 측정한다.

어떤 솔버(예: 직접해, CG, GMRES, 다중격자, 사후조건 Krylov 등)를 썼을 때,

  • 전체 계산 시간(초),

  • 반복 횟수(반복해 알고리즘인 경우),

  • 메모리 사용량(GB),

  • 최종 해의 정확도(오차 규준으로 측정)

  • 병렬화 스케일링(코어 수나 GPU 장치 수를 늘렸을 때 성능이 어떻게 변하는지)

와 같은 요인들을 수집하여, 각 솔버가 어느 구간에서 가장 유리한지 판단한다. 어떤 솔버는 작은 $n$에서 매우 빠르지만, $n$이 커지면 오히려 다른 솔버가 앞설 수도 있다. 또한 동일한 크기의 문제라도 계수 행렬의 구조(조건수, 희소 패턴, 주대각 우세 등)에 따라 결과가 크게 달라질 수 있다.

예를 들어 병렬화가 잘 되어 있는 LU 분해 루틴은 중간 규모에서는 최고의 성능을 보이다가, 초대규모로 가면 메모리 제약 때문에 “사전에 분해” 과정조차 힘들어질 수 있다. 이때 반복해 솔버는 분해 과정을 건너뛰고 곧바로 행렬-벡터 곱셈 연산 위주로 구성되어, 메모리 사용량이 적고 병렬 스케일링도 더 좋을 수 있다. 하지만 조건수가 나빠서 반복 횟수가 폭증하면 다시 계산 시간이 길어질 수 있다. 결국, 어떤 솔버가 최선인가 하는 문제는 “최소한의 튜닝이 이루어졌을 때 실제로 가장 빨리 문제를 풀 수 있는가”라는 관점에서 접근해야 한다.

반복해 솔버를 사용할 때에는 사후조건(preconditioner)의 선택이 엄청난 변수다. 사후조건자에 따라 실제 반복 횟수가 10배 이상 차이 나기도 하며, 이 사후조건자를 효율적으로 구성하는데 걸리는 오버헤드가 전체 계산에서 무시 못 할 비중을 차지할 수도 있다. Krylov 방법과 다중격자를 혼합한 알고리즘도 존재하는데, 다중격자로 저주파 오차를 제거하고 Krylov으로 고주파 성분을 세밀하게 완화하는 식으로 구성해 큰 시너지 효과를 내기도 한다. 그러나 이러한 고급 기법일수록 구현 복잡도와 파라미터 튜닝의 부담이 높아진다는 단점이 있다.

병렬 클러스터 환경에서의 분산 메모리 벤치마크 사례를 살펴보면, 수천 개 이상의 코어를 동원한 영역분할(Domain Decomposition) + Krylov + 사후조건 결합 기법이 매우 효율적인 성능을 발휘하는 경우가 많다. 행렬을 여러 서브도메인(subdomain)으로 나누고, 각 서브도메인에서 독립적으로 (블록) 직접해 분해를 진행하거나 간단한 반복해를 수행한 뒤, 전체 해를 한데 모아서 수정하는 절차를 반복한다. 이때 통신 비용과 각 노드의 캐시 사용 패턴 등을 최적화하여 높은 병렬 효율을 도출하는 것이 목표다. 예를 들어 특정 초대형 해석 문제에서, 다중 노드를 이용한 영역분할 + 다중격자 + Krylov 결합으로 천만 차원 이상의 문제를 수십 초 이내에 풀어내는 실제 사례도 보고되어 있다.

결국 이러한 벤치마크와 응용 사례들은 “이론적 복잡도 + 실제 구현상의 이점 + 병렬/메모리 구조 + 문제 특성”이라는 네 가지 축의 균형을 맞출 때, 가장 실용적인 결론이 도출됨을 보여준다. 본 장의 뒤에서는, 실제 예제 문제를 몇 가지 선정해 가우스 소거법, LU 분해, CG, GMRES, 다중격자 등을 비교하는 사례 연구를 진행할 것이다. 문제 크기, 행렬 구조, 하드웨어 스펙 등을 바꿔가며 측정하는 과정을 통해, 각 알고리즘이 언제 강점을 지니고 언제 약점을 보이는지 경험적으로 살펴보자.

고성능 계산 라이브러리와 구현 이슈

직접해든 반복해든, 대규모 선형시스템을 효율적으로 풀기 위해서는 고성능 선형대수 라이브러리를 적극적으로 활용하는 것이 매우 중요하다. 이러한 라이브러리는 보통 다음과 같은 특징을 갖는다.

  • 행렬 연산(예: 행렬-행렬 곱, 행렬-벡터 곱)을 고도의 최적화 기법과 캐시 구조에 맞춰 설계

  • BLAS(Basic Linear Algebra Subprograms)와 LAPACK(LA Package) 인터페이스 지원

  • CPU 멀티코어, 벡터화, GPU 가속 등 다양한 하드웨어를 고려한 병렬 루틴 제공

  • 희소(sparse) 연산을 위한 별도 루틴(MKL Sparse, CUDA Sparse, etc.) 마련

예컨대 오픈소스 진영에서는 OpenBLAS, ATLAS, BLIS 등이 있고, 상용 라이브러리로는 Intel MKL, NVIDIA cuBLAS, IBM ESSL 등이 널리 활용된다. 반복해 알고리즘의 경우, 행렬-벡터 곱셈이 전체 연산의 핵심이므로, 희소 행렬을 어떻게 효율적으로 저장하고 접근하는지가 성능에 큰 영향을 미친다. 대표적으로 CSR(Compressed Sparse Row) 혹은 CSC(Compressed Sparse Column) 포맷을 써서 실제 연산에 필요한 데이터만 빠르게 순회할 수 있도록 한다.

병렬화 관점에서는 OpenMP나 MPI를 통한 다중코어 및 분산 메모리 접근이 중요한데, 예를 들어 분산 메모리 환경에서 LU 분해를 적용하려면 행렬을 여러 노드에 나누어 저장하고, 각 단계마다 피벗팅과 부분행렬 갱신을 병렬화/통신화해야 한다. 반복해 알고리즘은 주로 행렬-벡터 곱셈을 반복 수행하기 때문에, 희소 구조 접근과 벡터 연산을 고루 병렬화/통신화한다. 이때 사후조건자 적용(예: 블록 대각, 영역분할)은 각 노드가 자기 영역에 해당하는 부분문제를 독립적으로 처리하고, 경계 정보를 교환하여 전체 문제의 해에 수렴하도록 설계한다.

혼합 정밀도(Mixed Precision) 기법

최근에는 혼합 정밀도(mixed precision) 기법이 GPU 등 고성능 하드웨어에서 매우 주목받고 있다. 이는 계산 과정 일부(또는 전부)에 단정밀도(single precision, 32비트) 또는 절반정밀도(half precision, 16비트)를 사용하여 연산을 가속하고, 오차가 크게 축적되지 않도록 필요한 단계에서만 배정밀도(double precision, 64비트) 연산을 수행함으로써 계산 시간을 크게 줄이는 방식이다.

예를 들어 반복해 알고리즘에서, 잔차가 아직 클 때에는 단정밀도 연산만으로도 큰 문제가 없을 수 있다. 오차가 어느 정도 줄어든 후에, 최종 해를 정교하게 구하기 위해 배정밀도 연산으로 전환한다. 이렇게 하면 배정밀도 연산만 계속 수행하는 경우에 비해 연산량이 줄고, 특히 GPU에서는 낮은 정밀도일수록 처리량(throughput)이 훨씬 커지는 특성이 있어 전체 계산 시간을 단축할 수 있다.

혼합 정밀도 전략을 적용할 때에는, 행렬의 조건수가 너무 좋지 않으면 낮은 정밀도 구간에서 오차가 쉽게 퍼질 수 있으므로, 반복의 각 단계나 사후조건자 구성에서 “언제 어떻게 정밀도를 전환할지”를 신중히 결정해야 한다. Krylov 계열의 반복해를 단정밀도에서 돌리다가, 최종 국면에 배정밀도로 스위칭하거나, 특정 스무딩이나 사후조건자 계산만 배정밀도로 처리하는 식의 선택지가 존재한다. 이 분야 역시 하드웨어 스펙과 문제 특성에 따라 최적 해법이 달라지므로, 실험적 튜닝이 필요하다.

에너지 최적화와 친환경 컴퓨팅

슈퍼컴퓨터나 대규모 데이터센터 환경에서 문제를 풀 때, “전력 소모를 최소화하거나(그린 컴퓨팅)” “동일 전력 내에서 최대 성능을 내는(성능-전력 효율)” 것이 또 다른 주요 이슈다. 같은 $O(n^3)$ 직접해 알고리즘이라도, 라이브러리에 따라 전력 사용 패턴이 크게 달라질 수 있고, 반복해 알고리즘도 병렬화 구조에 따라 전력 대비 성능 지표가 달라진다.

특히 행렬-벡터 곱셈이나 벡터 연산 위주로 구성된 방법이 메모리 대역폭과 GPU 연산 능력을 극대화하여 효율적으로 쓸 수 있다면, 전력 소비량을 상대적으로 낮추면서 높은 성능을 낼 수도 있다. 예컨대 GPU 클러스터에서 희소 행렬 연산을 잘 최적화해놓으면, 같은 문제를 CPU 클러스터 대비 훨씬 적은 전력으로 해결할 수 있다는 연구 결과도 다수 존재한다. 반면 너무 많은 통신이 수반되거나, 알고리즘 구조상 CPU 쪽에서 자주 동기화(sync)가 일어나면 오히려 에너지 효율이 떨어질 수 있다.

이러한 “에너지-성능” 트레이드오프 관점은 앞으로 더 중요해질 것으로 예측되며, 알고리즘 연구에서도 전력 소비 모델을 고려한 복잡도 해석(energy complexity)이나, 동적인 클록 주파수 제어(DVFS, Dynamic Voltage and Frequency Scaling)까지 포함한 최적화가 진행 중이다.

전문가 시스템과 자동 튜닝

라이브러리 단에서 BLAS/LAPACK 호출을 자동으로 최적화(tuning)하는 기능이 발전해 왔듯이, 반복해 알고리즘과 사후조건자 선택에서도 사용자 개입 없이 “문제 구조를 자동 분석해, 최적 솔버를 추천·세팅”해 주는 전문가 시스템이 연구되고 있다. 예를 들어 어떤 소프트웨어 프레임워크는

  • 행렬의 희소도

  • 대칭 여부

  • 양의 정부호 여부

  • 대각 우세성, 대각 블록 구조

  • 부분 도메인별 크기와 연결성

  • 기존에 알려진 문제군(case) 데이터베이스

등을 자동으로 스캔한 뒤, 가장 적절한 사후조건자 유형(블록 대각, incomplete LU, 다중격자, etc.)과 적당한 알고리즘(CG, BiCGSTAB, GMRES, MINRES 등)을 제안한다. 그 다음 몇 번의 테스트 런을 거쳐 반복 횟수와 계산 시간을 체크한 뒤, 실제 해석에 투입하는 식이다. 이런 접근은 여전히 초기에 세밀한 설정 및 사용자의 이해가 필요하지만, 장기적으로는 “자동 튜닝(auto-tuning)” 기술이 더욱 발전해, 초보 사용자도 대규모 문제를 비교적 간단히 해결하게 될 가능성이 높다.

케이스 스터디를 통한 비교

실제 사례 분석을 위해, 다음 장에서는 물리적 모델을 기반으로 한 구체적 예시를 두세 가지 들고, 여러 솔버를 각각 적용하는 실험을 진행할 것이다. 예컨대:

  • 2D 열전달 문제를 유한차분으로 이산화하여 얻은 포아송(Poisson) 타입 행렬

  • 3D 탄성해석 문제를 유한요소로 이산화하여 얻은 희소 행렬

  • 임의의 밀집행렬을 포함한 회귀분석 문제(작은·중간 규모)

등의 예시를 놓고, 가우스 소거법, $LU$ 분해, CG, GMRES(+사후조건자), 다중격자 등을 각각 구현·실행해 본다. 문제 크기($n$)를 서서히 늘리고, 병렬 코어 수나 GPU 사용 여부, 그리고 라이브러리 종류를 바꿔가며, 총 계산 시간과 반복 횟수, 메모리 사용량, 에너지 소비량까지 측정해본다. 이를 통해 이론적 복잡도와 실제 성능 사이에서 어떤 차이가 발생하고, 어떤 요인이 그 차이를 결정짓는지 좀 더 체감할 수 있을 것이다.

사례 연구 1: 2D 열전달(포아송) 문제

편미분방정식(PDE)에 기반한 선형시스템의 대표적 예시로 2차원 포아송(Poisson) 방정식을 들 수 있다. 열전달(steady-state heat conduction)이나 전기 퍼텐셜, 공기역학의 일부 문제 등이 사실상 다음 형태의 PDE로 귀결된다.

2u(x,y)=f(x,y)-\nabla^2 u(x, y) = f(x, y)

여기서 $u(x, y)$는 해석하고자 하는 물리량(온도, 퍼텐셜 등), $f(x, y)$는 문제 영역 내의 소스 항, 그리고 $\nabla^2$는 라플라시안(Laplacian) 연산자다. 경계조건(Dirichlet, Neumann 등)에 따라 문제 설정이 달라지며, 간단히 디리클레 경계를 예로 들면,

u(x,y)=g(x,y)onΩu(x, y) = g(x, y)\quad \text{on}\quad \partial\Omega

등의 조건이 붙는다. 이를 유한차분법이나 유한요소법으로 이산화하면, $n$차원의 희소 선형시스템

Au=f\mathbf{A}\mathbf{u} = \mathbf{f}

이 만들어진다. 행렬 $\mathbf{A}$는 대체로 대칭(symmetric), 양의 정부호(positive definite), 그리고 특정 패턴의 희소성(sparsity)을 갖는 특성을 지닌다. 다음은 가장 단순화된 5-점차분(5-point stencil) 유한차분을 직사격형 격자에 적용해 얻는 경우의 일반 형태다.

5-점차분에 의한 간단한 행렬 구조

2D 도메인 $\Omega$를 $h$ 간격으로 등분하여 내부 격자점을 $(i, j)$ ($i=1,\dots,N_x$, $j=1,\dots,N_y$)로 설정하면, $u_{i,j}\approx u(x_i, y_j)$에 대한 근사를 구하게 된다. 라플라시안은

2u(xi,yj)ui+1,j2ui,j+ui1,jh2ui,j+12ui,j+ui,j1h2-\nabla^2 u(x_i, y_j) \approx -\frac{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}}{h^2} - \frac{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}}{h^2}

과 같이 근사된다. 이에 따라 한 점 $(i, j)$에 대응되는 방정식을 적어 보면, 대략

ui+1,j+ui1,j+ui,j+1+ui,j14ui,j=h2fi,ju_{i+1,j} + u_{i-1,j} + u_{i,j+1} + u_{i,j-1} - 4u_{i,j} = -h^2 f_{i,j}

형태가 되어, 5-점차분 방식임을 알 수 있다. 이렇게 전체 격자점에 대해 식을 쓰면, 희소 행렬 $\mathbf{A}$는 대각 원소와 소수의 인접 원소들만 값이 있는 “밴드 구조”를 갖게 된다. 정확히 말해서 각 행에 최대 5개(본인 및 상하좌우 이웃)의 비영(非零) 원소가 나타난다.

Python 예시 코드로 살펴보는 행렬 생성

아래 코드는 Python에서 SciPy의 희소행렬(sparse matrix) 기능을 간단히 활용해 2D 포아송 문제를 이산화하고, 이를 여러 솔버로 풀 수 있도록 준비하는 골격 예시다. 실제로는 경계 조건이나 소스 항 $f_{i,j}$ 처리 등을 좀 더 정교하게 구현해야 한다.

위 코드에서 generate_poisson_2d 함수는 대칭 양의 정부호 희소행렬 A와 우변 벡터 b를 반환한다. 이후 spsolve는 내부적으로 희소행렬 $LU$ 분해를 수행하는 것으로 알려져 있으며, cg는 SciPy에서 제공하는 단순 CG(Conjugate Gradient) 반복 솔버다. 이 예시만으로는 조건수나 수렴 속도, 실제 연산 시간을 충분히 보여주지 않으므로, 더 큰 문제 크기나 다른 사후조건자, 다중격자 라이브러리 등을 적용해가며 성능 차이를 관찰할 수 있다.

이론적 복잡도 vs. 실제 계산 시간

5-점차분 포아송 문제의 경우, 내부 격자점이 대략 $N_x \times N_y$개이므로, $n = N_x N_y$가 전체 미지수 차원이 된다. 희소행렬 $\mathbf{A}$의 각 행에는 최대 5개의 비영 원소만 있으므로, 전체적으로 $O(n)$ 개의 비영 원소를 가진다고 볼 수 있다. 따라서 CG나 GMRES 같은 반복해에서 한 번의 행렬-벡터 곱셈이 $O(n)$에 수행될 수 있고, 다중격자를 포함하면 대부분 최종 해결까지 $O(n)$ 혹은 $O(n \log n)$ 수준의 연산량으로 끝낼 수 있다. 이론적으로는 매우 효율적인 성능을 기대할 수 있다.

하지만 실제 계산 시간은 다음과 같은 요인에 의해 달라진다.

직접해(spsolve)를 사용할 경우, 희소 $LU$ 분해 루틴이 내부적으로 어떤 알고리즘을 쓰느냐에 따라 성능이 달라진다. 일반적으로 스파스 $LU$ 분해에서 “fill-in(분해 과정에서 0이 아닌 요소가 새로 생기는 현상)”이 발생하면, $O(n^{1.5})$ 혹은 $O(n^2)$ 수준의 메모리와 연산을 요구하기도 한다.

반복해(CG 등)에서는 한 번 반복당 $O(n)$ 연산이 가능하다 하더라도, 실제 필요한 반복 횟수 $K$가 조건수나 경계조건 등에 따라 달라진다. 사후조건자를 사용하지 않으면 $K$가 커져 전체 $O(Kn)$이 $O(n^{1.5})$ 정도로 악화되는 경우도 종종 있다.

병렬환경에서 행렬-벡터 곱셈을 어떻게 분할하고, 통신비용을 얼마나 줄이느냐도 중요하다. $\mathbf{A}$가 행기반(block row)으로 분할되어 각 프로세스가 자기 부분을 처리한다면, 프로세스 간 경계 행 정보만 교환하면 되므로, 높은 병렬 효율을 얻을 수 있다. 그러나 누락된 사후조건자나 부적절한 병렬 스케줄링은 수렴 속도를 떨어뜨리거나 통신 병목을 야기한다.

후속 연구와 사례 확장

2D 포아송 예시는 “가장 기본적이고 단순한 구조”를 갖고 있기 때문에, 수치해석 이론과 실제 계산이 비교적 잘 부합하는 편이다. 그러나 실제 산업 및 학술 연구 현장에서는 더 복잡한 물리, 불규칙한 도메인, 비정형 격자, 비선형 경계 조건 등이 합쳐져서, 행렬 $\mathbf{A}$의 희소 패턴이 예측하기 어려워지고, condition number도 훨씬 커질 수 있다. 그럴 때는 반복해 솔버를 위한 효과적인 사후조건자나 다중격자 기법을 적극적으로 도입해야 한다.

예를 들어 3D 문제로 넘어가면, 행렬 크기가 $N_x \times N_y \times N_z$로 급격히 증가하고, 실질적으로 $n = O(N^3)$인 수준이 되므로, 직접해는 매우 부담스러워진다. 이때 반복해 솔버(특히 다중격자)는 이론적으로 $O(n)$ 또는 $O(n \log n)$의 연산량으로 해결할 수 있어, 대규모 3D 문제의 표준적 접근법으로 받아들여지고 있다.

다음 사례 연구에서는 조금 더 복잡한 편미분방정식(예: 탄성 역학, 유체 역학)이나, 밀집 구조를 갖는 문제에 대해서도 비슷한 비교를 진행하면서, 실제 실험 결과를 어떻게 해석해야 하는지 살펴보기로 하자.

사례 연구 2: 3D 탄성해석 문제

2D 포아송 문제는 비교적 단순한 PDE 사례이지만, 실제 공학 현장에서는 훨씬 복잡한 물리계가 주어지는 경우가 많다. 예를 들어 탄성해석(Elasticity)은 구조물의 변형과 응력/변형률 분포를 구하는 대표적 분야다. 3차원 도메인에서 물체가 외력을 받았을 때 발생하는 변위(Displacement)와 응력(Stress)을 계산하기 위해서는, 라메(Lamé) 계수나 탄성 계수, 포아송 비 등 물성을 고려한 탄성 방정식을 풀어야 한다. 이를 유한요소법(Finite Element Method, FEM)으로 이산화하면, 대규모 희소 선형시스템

Ku=f\mathbf{K}\mathbf{u} = \mathbf{f}

이 생성된다. 여기서 $\mathbf{K}$는 보통 ‘강성행렬(stiffness matrix)’이라 불리고, $\mathbf{u}$는 구조물의 변위 벡터, $\mathbf{f}$는 외력 벡터에 해당한다.

유한요소법으로 구성된 강성행렬의 특성

3D 탄성문제에서, 요소(element)가 헥사(Hexa), 사면체(Tetra), 사각형 기반의 복합 요소 등 다양하게 설정될 수 있다. 요소의 형태와 개수에 따라 유한요소 해석 결과가 달라지는데, 일반적으로 다음 특징을 가진 행렬이 만들어진다.

  • 매우 큰 차원. 실제 산업적 문제에서는 요소 개수가 수백만에서 수천만까지도 올라가므로, 미지수($n$)의 규모도 비례해서 수백만 이상이 될 수 있다.

  • 희소 구조. 각 요소는 제한된 개수의 자유도(예: 헥사 요소 한 개당 8개의 정점 x 3방향 변위 = 24 자유도)를 연결하므로, 각 행에 등장하는 비영 원소가 상대적으로 제한적이다. 전체적으로 $O(n)$ 정도의 비영 원소를 갖는 희소행렬이 형성되는 경우가 많다.

  • 대칭 양의 정부호(symmetric positive definite). 탄성해석에서 $\mathbf{K}$는 일반적으로 이 조건을 만족한다(비선형 접선강성이나 접촉문제 등을 가정하지 않는다면). 이는 CG(Conjugate Gradient) 계열의 반복해 사용에 유리한 요인이다.

  • 조건수. 구조 해석에서 $\mathbf{K}$가 상당히 고르지 못한(ill-conditioned) 양상을 보이기도 한다. 예를 들어 매우 얇은 부재나 서로 큰 강성 차이가 있는 재질을 동시에 고려하면, 어떤 부분은 강성이 매우 크고 어떤 부분은 매우 작아서, 수치상으로 불안정한 시스템이 될 수 있다. 이때 사후조건자(preconditioner) 없이 CG를 쓰면 반복 횟수가 크게 늘어난다.

직접해 vs. 반복해의 선택 고려

중간 규모(예: 수십만 자유도) 문제에서는, 적절하게 구현된 희소 $LU$ 분해 루틴(또는 희소 Cholesky 분해)이 여전히 효과적일 수 있다. 그러나 수백만~수천만 규모로 넘어가면, 분해 과정에서 발생하는 fill-in 때문에 $O(n^{1.5}\sim n^2)$ 이상의 연산량과 대규모 메모리를 요구하게 된다. 실제로 이런 초대형 문제에서 직접해를 시도하면, 고성능 클러스터를 동원해도 어마어마한 계산 시간이 들거나, 메모리 한계로 인해 실패할 수 있다.

따라서 대규모 3D 탄성해석에는 반복해 알고리즘이 사실상 표준으로 자리잡는다. 특히 CG, MINRES, PCG(Preconditioned CG), FETI(Finite Element Tearing and Interconnecting) 방법 등 다양한 Krylov 서브스페이스 기법이 적용된다. FETI 계열 기법은 영역분할(domain decomposition)과 연결되어, 큰 구조물을 여러 부분 도메인으로 나누고 경계에서 조건을 맞춰주면서 각 부분문제를 병렬로 해결하는 방식이다.

사후조건자 선택도 매우 중요한데, 블록 대각 사후조건자나 Incomplete Cholesky(IC) 계열, 다중격자(Multigrid) 방식 등이 활용된다. 경우에 따라서는 AMG(Algebraic Multigrid)처럼 격자 생성 과정을 자동화한 다중격자 기법을 쓰기도 한다.

3D 문제에서의 연산 복잡도와 병렬화

이론적으로 $\mathbf{K}$가 대칭 양의 정부호라면, CG 반복 한 번당 $O(n)$ 연산(희소 행렬-벡터 곱셈, 내적, 벡터합 등)으로 구성될 수 있다. 그러나 반복 횟수 $K$가 조건수에 의해 좌우되므로, 최종 복잡도는 $O(Kn)$이다. 사후조건자가 잘 작동하면 $K$를 거의 상수 수준(또는 다소 완만하게 증가)으로 억제할 수 있으므로, 이상적으로 $O(n)$ 또는 $O(n\log n)$ 수준으로 전체 문제를 해결할 수 있다.

하지만 실제로는 병렬 환경에서 각 반복마다 경계나 전역 데이터를 통합·교환하는 통신 비용이 추가된다. CPU 클러스터의 MPI 통신, GPU 간에 걸친 NVLink나 Infiniband, 혹은 하이브리드 아키텍처에서 노드 간 통신은 병렬 효율을 결정짓는 주요 요인이 된다. 영역분할법을 적용하면 각 서브도메인의 자유도 수가 감소해 부분문제를 빠르게 푸는 대신, 서브도메인 경계에서의 통신이 반복적으로 일어난다. 최적의 영역분할이 이뤄지지 않으면 통신 부하가 커지거나, 각 노드가 맡은 연산량이 불균등해져 전체 계산 시간이 늘어난다.

3D 탄성해석에서 흔히 사용하는 다중격자(Geometric Multigrid) 또는 AMG(Algebraic Multigrid) 기법 역시 큰 힘을 발휘한다. 하위 해상도의 격자(또는 알gebraic 준거)를 구성하여, 고주파 오차를 상위 레벨에서 제거하고 저주파 오차는 하위 레벨에서 제거함으로써 반복 횟수를 크게 줄일 수 있다. 특히 양의 정부호 대칭행렬에서는 MG나 AMG가 CG와 자연스럽게 결합되어 PCG(Preconditioned CG) 방식으로 사용되는 사례가 많다.

수치 안정성과 접촉/비선형성

단순 선형 탄성해석이 아니라 접촉(contact) 문제나 큰 변형(large deformation)을 포함한 비선형 탄성(non-linear elasticity) 문제로 넘어가면, 한 번의 뉴턴(Newton) 반복에서 접선강성행렬(tangent stiffness matrix)이 달라지는 구조가 된다. 이때도 선형시스템 풀이는 필수이며, 전체 해석 과정은 다음과 같은 형태로 이뤄질 수 있다.

  1. 현재 변형 상태에서 접선강성행렬 $\mathbf{K}(\mathbf{u}^{(k)})$ 생성

  2. 선형시스템 $\mathbf{K}(\mathbf{u}^{(k)}) \Delta \mathbf{u} = \mathbf{r}(\mathbf{u}^{(k)})$ 풀기

  3. $\mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} + \Delta \mathbf{u}$ 갱신

  4. 수렴 시까지 반복

여기서 2번 단계가 곧 선형해석이며, 문제 크기가 일정히 크므로 다시 $O(n)$ 혹은 $O(Kn)$의 반복해 솔버를 효율적으로 써야 한다. 매 스텝마다 $\mathbf{K}(\mathbf{u}^{(k)})$가 바뀌므로, $LU$ 분해 같은 직접해를 재활용하기 어려워지는 반면, 반복해 알고리즘은 큰 부하 없이 새로운 행렬에 적용 가능하다는 장점이 있다. 다만 비선형 문제는 해가 점진적으로 바뀌므로, 사후조건자도 유연하게 조정하거나 재구축해야 할 수 있고, 이는 추가적인 비용이 될 수 있다.

접촉 문제(Contact mechanics)에서는 경계 조건이 “비활성”에서 “활성”으로 바뀌는 식으로, 해석 도중에 행렬 구조가 변화하기도 한다. 이 경우, 특정 부분 영역의 행렬이 새로 커플링되면서 $\mathbf{K}$의 희소 패턴마저 변할 수 있다. 이런 환경에서 직접해는 계속해서 분해를 갱신해야 하므로 매우 비효율적일 수 있고, 반복해는 구조 변화를 온전히 수용하면서 계속해서 업데이트를 진행할 수 있다는 이점이 있다.

간단한 코드 스케치와 확장

3D 탄성해석 전체를 Python이나 Octave 코드로 예시하기는 매우 방대하므로, 간단히 “어떤 식으로 반복해를 적용할 수 있는지”만 스케치해보면 다음과 같은 흐름을 갖는다.

generate_3d_elastic_stiffness 함수 내부는 유한요소 라이브러리(예: Fenics, Deal.II, CalculiX 등)를 통해 자동으로 구성될 수 있다. AMGPreconditioner는 사용 라이브러리에 따라 구체적인 구현이 달라지며, SciPy나 PyAMG, Trilinos, PETSc 같은 라이브러리를 활용하면 다양한 고성능 사후조건자나 Krylov 솔버를 활용할 수 있다.

이 과정을 병렬 클러스터 환경으로 확장하려면, MPI나 OpenMP 지원이 있는 라이브러리를 써야 하며, PETSc, Trilinos, hypre 등이 대표적인 예다. 실제로 대규모 산업 문제에서는 수백~수천 개 이상의 코어를 동원하여 FETI, BDDC(Balancing Domain Decomposition by Constraints), 혹은 AMG+Krylov 솔버를 조합해 문제를 해결한다.

전체적 시사점

3D 탄성해석 문제는 2D 포아송 문제보다 훨씬 복잡하고, 실제 문제가 대규모로 확장되면서 “반복해 + 사후조건자 + 병렬화” 구도가 확연히 중요해진다. 특히 구조에 따라 해석 시간을 $10$배 이상 단축하거나 수렴이 불가능해지는 경우도 흔하므로, 문제 특성과 솔버/사후조건자 선택이 매우 민감하다. 그런 이유로 수치해석과 고성능 계산 분야에서 끊임없이 새로운 기법들이 연구되고 있으며, 실제 상업용 소프트웨어(예: Abaqus, ANSYS, LS-DYNA, COMSOL) 내부에도 여러 솔버가 탑재되어 사용자의 문제 특성에 따라 적절히 선택·적용되고 있다.

다음 사례에서는 이러한 편미분방정식 기반 문제가 아닌, 밀집(dense) 행렬이 자연스럽게 발생하는 또 다른 예를 들어본다. 예컨대 통계학적 회귀분석이나 기계학습, 또는 특정 커널 기법에서 $n\times n$ 완전 밀집행렬을 다루는 상황을 살펴볼 것이며, 이때에는 직접해의 $O(n^3)$ vs. 반복해의 $O(Kn^2)$ 비교가 다른 양상으로 전개된다는 점을 짚어볼 것이다.

사례 연구 3: 밀집행렬에서의 선형시스템

지금까지는 PDE(편미분방정식) 기반의 해석 문제를 다루었으며, 이때 얻어지는 계수 행렬은 대체로 희소(sparse) 구조를 갖는다는 특징이 있었다. 반면, 통계학이나 기계학습, 최적화 문제 등에서는 완전 밀집(dense) 행렬이 자연스럽게 형성되기도 한다. 이 경우, 행렬 크기가 $n\times n$이라면 모든 원소가 0이 아니거나(또는 거의 대부분이 0이 아님) 특별한 희소 구조 최적화를 기대하기 어렵다. 대표적으로 대규모 회귀분석(예: 최소제곱법), 커널 방법(kernel methods), 또는 전치행렬과 정방행렬 간의 고차원 곱셈이 얽힌 문제에서 이러한 상황이 발생한다.

예시 1: 최소제곱법(OLS)에서의 밀집행렬

전통적인 다중선형회귀(Multiple Linear Regression) 모델을 생각해보자. 독립변수 벡터가 $x \in \mathbb{R}^p$, 종속변수가 스칼라인 $y$로 주어졌을 때, 표본 데이터가 $n$개 있으면, 이를 행렬 $\mathbf{X}\in\mathbb{R}^{n\times p}$와 벡터 $\mathbf{y}\in\mathbb{R}^n$로 정리할 수 있다. 일반적인 최소제곱(OLS) 문제는

minβRpyXβ22\min_{\beta \in \mathbb{R}^p} \, \|\mathbf{y} - \mathbf{X}\beta\|_2^2

의 형태가 되며, 통상적으로 $p \le n$이라 해도, 보통 $\mathbf{X}$는 밀집행렬로 간주하거나, 데이터 크기가 아주 커서 희소 행렬 기법을 적용하기 애매한 경우가 있다.

가장 단순한 직접해는 정규방정식(normal equation)

XXβ=Xy\mathbf{X}^\top \mathbf{X}\,\beta = \mathbf{X}^\top \mathbf{y}

을 풀어 $\beta = (\mathbf{X}^\top \mathbf{X})^{-1} \mathbf{X}^\top \mathbf{y}$를 구하는 것이다. 여기서 $\mathbf{X}^\top \mathbf{X}$는 $p\times p$ 행렬이지만, 만약 $p$가 수십만, 수백만에 이를 정도로 매우 크다면, 전체가 밀집 행렬이므로 직관적인 $O(p^3)$ 연산이 부담스러워진다. 더구나 수치 안정성 관점에서도, $\mathbf{X}^\top \mathbf{X}$가 잘못된 조건수를 가질 수 있어 직접 역행렬을 구하는 방식은 별로 권장되지 않는다.

이 문제를 더 효율적으로 풀기 위해 QR 분해나 SVD 분해를 적용하는데, 이들 역시 전형적으로 $O(np^2)$, $O(p^3)$ 혹은 이와 유사한 차수의 연산이 필요하다(단, $n \gg p$의 경우 $O(np^2)$가 지배적일 수 있다). 만약 $p$가 커서 이 방식조차 버겁다면, 반복해 기법(예: LSQR, LSMR 등)을 활용하거나, 랜덤화된 저랭크 근사(randomized low-rank approximation) 기법 등을 통해 연산을 줄이기도 한다.

예시 2: 커널 방법에서의 밀집행렬

기계학습의 커널 기법(kernel method) 중 대표적인 서포트벡터머신(SVM)이나 가우시안 프로세스(Gaussian Process) 회귀 등을 생각해보면, 훈련 표본이 $n$개 있을 때, 커널행렬(kernel matrix) $\mathbf{K}\in\mathbb{R}^{n\times n}$을 구성해야 하는 상황이 온다. 보통 커널함수 $\kappa(x_i, x_j)$를 모든 쌍 $(i, j)$에 대해 계산해 $\mathbf{K}_{ij} = \kappa(x_i, x_j)$로 설정한다. 이 행렬은 대칭 양의 정부호(symmetric positive definite)지만, 일반적으로 완전 밀집행렬에 해당한다.

또한 가우시안 프로세스 회귀를 수행할 때에는 $\mathbf{K} + \sigma^2 \mathbf{I}$ (노이즈 항 포함)을 역행렬 혹은 분해 형태로 다루게 되는데, $n$이 수천~수만 단위일 때조차 $O(n^3)$ 연산과 $O(n^2)$ 메모리가 필요하므로, 현실적으로 매우 부담스럽다. 이러한 이유로 “커널법은 $n$이 크면 사실상 쓰기 어렵다”는 말이 종종 나온다.

물론, 근사 기법이 다양하게 개발되어 있다. 예: 랜덤 피처(RF, Random Fourier Features)나 Nystrom 근사, 혹은 특정 구조(저랭크 근사, 쿼드 트리 분할, FMM(Fast Multipole Method) 유사 기법 등)를 이용해, $O(n^2)$ 혹은 $O(n\log n)$ 정도로 줄이는 시도들이 계속되고 있다. 하지만 이런 근사 없이 정면으로 커널행렬을 분해하거나 역행렬을 구하려면, $n$이 수만만 되어도 $O(n^3)$ 복잡도가 크나큰 장벽이 된다.

직접해와 반복해의 활용도

밀집행렬 문제에서는 “행렬-벡터 곱셈이 $O(n^2)$을 요한다”는 점이 희소행렬과의 결정적 차이점이다. 크기가 $n\times n$인 밀집행렬 $\mathbf{A}$에 대해 $\mathbf{A}\mathbf{x}$를 구하면, 곱셈 연산이 대략 $n^2$ 스케일로 발생한다. 따라서 반복해 알고리즘이 한 번 반복할 때마다 $O(n^2)$을 소모한다. 반복 횟수가 $K$라면 전체는 $O(K n^2)$ 복잡도다.

직접해로 $\mathbf{A}$를 분해($LU$, $QR$, $SVD$ 등)하는 경우, 그 분해 과정이 $O(n^3)$에 달하지만, 분해가 끝나면 $\mathbf{b}$가 바뀌어도 $O(n^2)$만에 해 $\mathbf{x}$를 구할 수 있다. 문제에 따라 “분해를 한 번 해두고 우변만 여러 번 바뀌는” 시나리오라면, 직접해가 훨씬 실용적일 수 있다. 또한 반복해는 매번 $O(n^2)$ 연산을 반복해야 하므로, $K$가 수십 회만 되어도 금세 $O(n^3)$에 육박하게 되며, 조건수가 나쁘면 $K$가 더 늘어난다.

이 때문에 밀집행렬 문제에서 “한 번만 해를 풀면 되는” 상황이라면, 오히려 잘 최적화된 직접해가 앞설 수 있다. BLAS 레벨-3 연산(행렬-행렬 연산 등)을 통해 $O(n^3)$ 알고리즘이 GPU에서 매우 빠르게 수행될 수도 있다. 예컨대 $n$이 수천에서 수만 정도인 경우, 고성능 GPU에서는 행렬 분해가 의외로 빠르게 끝나기도 한다. 반면 반복해는 빅오 표기로만 보면 $O(n^2)$ per iteration으로 좋아 보일 수 있겠지만, 실제로는 $K$에 따라 총연산이 $O(n^3)$ 이상이 될 위험이 있다.

대규모 밀집행렬에서의 근사 기법

밀집행렬의 크기가 정말 커서(수십만 이상) $O(n^3)$나 $O(n^2)$ 계산 자체가 불가능한 수준이라면, 근사 기법을 고려해야 한다. 예를 들어,

  • 저랭크 근사(low-rank approximation): $\mathbf{A} \approx \mathbf{U}\mathbf{V}^\top$ 꼴($\mathbf{U}$, $\mathbf{V}$는 $n\times k$), 여기서 $k \ll n$

  • H-행렬(H-matrix) 계열: 계층적(block) 구조를 이용해 빠른 곱셈 및 분해

  • FMM(Fast Multipole Method) 유사 아이디어: 물리적 거리나 상호작용 구조를 활용

  • 랜덤화 기법: $\mathbf{A}$를 랜덤 투영으로 줄여, 근사 문제를 빠르게 해결

등이 존재한다. 이러한 기법들은 $O(n^2)$에서 더 나아가 $O(n \log n)$ 또는 $O(n)$에 가까운 계산량으로 행렬 연산을 처리할 수 있는 가능성을 열어준다. 단, 이때엔 근사 오차가 발생하므로, “어느 정도 정확도를 희생하더라도 빠른 계산을 택하겠다”라는 응용 목표가 뚜렷해야 한다.

정리하면, 밀집행렬 문제에서는

  • “하나의 문제를 매우 정확하게 풀어야 하고, $n$이 중간 수준(수천~수만)”이라면 직접해($O(n^3)$)가 병렬/가속화와 결합되어 빠르고 안정적인 해법이 될 수 있다.

  • “우변 벡터가 여러 번 바뀌면서 계속 해를 구해야 한다”면, 한 번 분해로 여러 번 해가 가능하므로 직접해가 유리하다.

  • “$n$이 엄청나게 커서 $O(n^3)$도 수행이 불가능”하다면, 반복해를 쓰기보다는 근사 기법(저랭크, FMM, etc.) 또는 랜덤화 기법을 결합하여 $O(n^2)$ 이하의 알고리즘을 고안해야 할 것이다.

  • “밀집 + 반복해” 방식은 $K$가 충분히 작으면 $O(Kn^2)$가 가능하지만, 조건수가 나쁘면 $K$가 커져버리므로 사실상 직접해와 큰 차이가 없거나 더 느려지는 경우가 흔하다.

실험 예시와 시사점

통계/머신러닝 측면에서, $n$ 또는 $p$가 수천수만 범위라면, NumPy/SciPy, MATLAB, R 등에서 LAPACK 기반의 $O(n^3)$ 분해가 GPU나 병렬 라이브러리로 비교적 빠르게 동작한다. 예컨대 $n = 10^4$ 정도라면, 수 초수십 초 안에 분해가 끝나기도 한다(고성능 GPU 사용 시). 조건수가 어느 정도 양호하다면 수치 안정성도 문제없다.

하지만 $n = 10^5$ 이상을 밀집행렬로 다룬다면, 아무리 고성능 GPU가 있어도 $O(10^{15})$ 스케일의 연산량은 사실상 절망적이며, 메모리도 TB 단위를 넘어설 수 있다. 그래서 근본적으로 데이터 차원을 축소하거나, 커널 저랭크 근사 기법을 적용하는 등 발상의 전환이 필요해진다. 이때 반복해를 무턱대고 시도하기보다는, “행렬의 특정 구조”나 “물리적/통계적 근사”를 적극 활용하여 계산량을 획기적으로 줄이는 전략이 권장된다.

밀집행렬 문제는 PDE 기반 문제와 달리, 희소 구조 최적화가 힘들다는 점에서, 이론적 복잡도와 실제 구현의 부담이 확연히 달라진다. 결국 “중규모 밀집행렬”에선 직접해가 강세를 보이지만, “초대형 밀집행렬”은 반복해와 근사기법, 혹은 문제 자체의 구조적 특성을 통한 차원 축소를 모색해야 한다는 점이 중요한 시사점이다.

메모리 구조와 캐싱 효과

수많은 선형대수 알고리즘에서 이론적 연산 복잡도가 동일하더라도, 실제 계산 시간은 메모리 접근 패턴이나 캐시 사용 효율에 따라 크게 달라진다. 현대 CPU나 GPU는 연산 속도(FLOPS) 자체도 중요하지만, 동시에 메모리 대역폭(memory bandwidth)과 캐시 계층 구조를 어떻게 활용하느냐가 알고리즘의 실제 성능에 지대한 영향을 미친다. 이를테면 다음과 같은 요인들이 핵심이다.

  • 연속적 배열 접근: 벡터나 행렬 원소를 연속적인 메모리 구역에 배치하면, CPU 캐시나 GPU 텍스처/공유 메모리 등의 효율이 극대화되어 성능 향상에 유리하다.

  • 블록 연산(blocked operations): 행렬 곱셈이나 분해를 할 때, 큰 연산을 여러 작은 블록으로 나누어 처리하면 캐시에 더 잘 적재되어 반복 사용될 수 있다. BLAS-3 연산(예: GEMM, GEMV)에서 이러한 블로킹(blocking) 기법이 광범위하게 사용된다.

  • 희소 행렬 포맷: CSR(Compressed Sparse Row), CSC(Compressed Sparse Column), ELLPACK, COO 등 다양한 희소 행렬 저장 포맷이 있으며, 각 포맷은 특정 유형의 희소 패턴이나 하드웨어 구조에 특화되어 있다. 제대로 된 포맷을 선택하지 않으면, 실제 계산이 이론적으로는 $O(n)$이라도 계속 랜덤 접근(random access)이 발생해 캐시 미스(cache miss)가 빈번하게 일어날 수 있다.

  • NUMA(Non-Uniform Memory Access) 구조: 멀티소켓 CPU 시스템에서, 각 소켓의 로컬 메모리에 접근하는 것과 다른 소켓 메모리에 접근하는 것의 지연 시간이 다를 수 있다. 선형시스템 알고리즘이 여러 소켓에 걸쳐 동작할 때, 데이터가 어느 소켓에 배치되었는지에 따라 성능 차이가 커질 수 있다. MPI나 OpenMP를 이용해 병렬 처리를 할 때 이러한 NUMA 환경을 잘 고려해야 한다.

특히 대규모 반복해 알고리즘에서 매 반복마다 행렬-벡터 곱셈을 수행할 때, 희소 구조와 포맷에 따라 캐시 적중률이 크게 좌우된다. 이를 개선하기 위해선, 도메인 분할과 함께 로컬 도메인 내의 행렬을 연속된 메모리에 저장하거나, 행렬을 블록-CSR 형태로 변환하는 방법 등을 활용한다. GPU 가속 시에도 공유 메모리(shared memory)나 텍스처 캐시 활용 정도가 알고리즘 성능을 좌우하므로, 단순히 이론적 연산 횟수만 볼 것이 아니라 “메모리 접근 최적화 전략”을 비롯한 구현 세부 요소가 중요하다.

분산 메모리와 대규모 병렬 환경

오늘날 초고성능 계산(HPC)을 수행하는 슈퍼컴퓨터나 클러스터는 대규모 노드가 네트워크로 연결된 분산 메모리 환경을 이룬다. 한 노드는 여러 개의 CPU 코어와 GPU를 포함할 수 있으며, 각 노드 간에는 MPI(Message Passing Interface) 등의 통신 프로토콜을 통해 데이터를 주고받는다. 이런 분산 메모리 구조에서 선형시스템을 풀 때는 다음과 같은 추가 이슈가 등장한다.

  • 행렬 분할(matrix partitioning): 전체 행렬을 어떻게 노드별로 나눌지 결정해야 한다. 예: 행 단위로 나누어 $A$를 블록 행렬로 만들고, 각 노드가 자기 행에 해당하는 부분을 저장하고 연산을 담당한다. 반복해 알고리즘에서는 행렬-벡터 곱셈 시 각 노드가 부분결과를 계산한 뒤, 경계(다른 노드에 저장된 요소) 정보를 교환해야 한다.

  • 통신 오버헤드: 노드가 늘어날수록, 통신량과 동기화 비용이 커진다. 간단한 CG 반복에서조차 내적(dot product)을 위해 전역적으로 모든 프로세스가 합산(reduction) 동작을 하게 되고, 이것이 병목이 될 수도 있다. 영역분할법이나 다중격자에서도, 각 부분 도메인 경계에서 정보 교환이 자주 일어난다.

  • 동적 부하 분산(load balancing): 행렬이나 격자 영역을 균등하게 나눴다고 해서, 반드시 각 노드가 같은 연산량을 갖는 것은 아니다. 실제 연산 시간은 캐시 히트율, 데이터 밀도, 네트워크 지연 등에 영향을 받으므로, 어떤 노드는 일찍 연산을 마치고 대기하는 반면, 다른 노드는 늦게 끝날 수 있다. 이를 최소화하기 위한 부하 분산 전략이 필요하다.

결국 “이론적 연산 복잡도가 $O(n^3)$ 또는 $O(n^2)$라 하더라도, 실제로는 분산 메모리 구조에서 통신 비용과 부하 분산이 얼마나 최적화되었는지가 성능을 결정짓는다”는 사실이 중요하다. 일반적으로 크기가 $n$인 문제를 $p$개의 노드에 분산하면, 순수 계산량은 $\frac{1}{p}$로 줄어들지만 통신과 동기화가 추가되므로, 이상적인 선형 스케일링(Linear Scaling)을 달성하기는 쉽지 않다. 고성능 라이브러리(PETSc, Trilinos, hypre 등)는 이 문제를 해결하기 위한 다양한 병렬 알고리즘 및 자동 영역분할 기능을 제공하고 있다.

동적 정밀도와 적응형 전략

앞서 소개한 혼합 정밀도(mixed precision) 개념을 확장하면, 알고리즘이 진행되면서 정밀도와 사후조건의 수준을 동적으로 바꾸는 적응형(adaptive) 전략도 고려할 수 있다. 예를 들면, 초기에 오차가 큰 구간에서는 단정밀도(single precision)로 빠르게 큰 폭의 오차를 줄이고, 어느 정도 수렴이 이뤄지면 배정밀도(double precision) 연산으로 전환해 최종 해를 정교하게 다듬는 식이다. 혹은 사후조건자 역시 초기에 대략적인(가벼운) 사후조건을 적용하다가, 수렴이 더딘 특정 단계에서만 좀 더 정교한 사후조건자를 구축·적용해 반복 횟수를 줄이도록 설계할 수도 있다.

이와 같이 “계산 자원을 적재적소에 배분”하는 방식은, 최근 HPC나 딥러닝 분야에서 큰 관심을 끌고 있다. 딥러닝 최적화에도 혼합 정밀도 기법이 폭넓게 쓰이는데, 선형시스템 문제에서도 유사한 아이디어를 적용하면 상당한 계산 시간 단축과 에너지 절감을 기대할 수 있다. 다만, 이러한 적응형 전략은 구현 난이도가 높으며, 실시간으로 상태를 모니터링하고 결정을 내리는 로직이 필요하므로, 상대적으로 복잡한 소프트웨어 구조가 요구된다.

소프트웨어 생태계와 알고리즘 선택

최적의 알고리즘을 택하더라도, 실제 현장에서 구현과 관리가 용이하지 않다면 실용화가 어려울 수 있다. 따라서 많은 수치해석 소프트웨어는 알고리즘 라이브러리화와 함께, 사용자 친화적인 인터페이스를 제공한다. 예를 들어:

  • 수치해석 라이브러리: BLAS/LAPACK, ScaLAPACK, PETSc, Trilinos, MUMPS, SuperLU, MKL, cuBLAS, cuSOLVER 등

  • FEM/CFD 플랫폼: deal.II, FEniCS, OpenFOAM, FreeFEM 등 (유한요소·유체해석 특화)

  • 상업용 패키지: MATLAB, Abaqus, ANSYS, COMSOL, LS-DYNA, NASTRAN 등

이들은 내부적으로 다양한 선형시스템 알고리즘을 포함하고, 문제 유형에 맞게 “자동 선택”하거나, 사용자 매뉴얼에서 “권장 알고리즘”을 제시한다. 예컨대, PETSc나 Trilinos는 CG/BiCGSTAB/GMRES/TFQMR 등 Krylov 계열과, 여러 사후조건자를 메뉴처럼 제공하며, 사용자가 옵션을 통해 조합하여 시도할 수 있도록 해준다. FEM 툴도 계산 백엔드를 PETSc나 Trilinos와 연결해 다중격자나 영역분할 기법을 투명하게 적용한다.

이 생태계에서 핵심은 “다양한 선택지 중, 어느 것을 언제·어떻게 쓰면 되는지 직관을 얻는 것”이다. 알고리즘 전문가나 응용 분야 엔지니어는, 특정 문제의 규모($n$), 구조(희소/밀집), 조건수, 하드웨어 자원, 정밀도 요구사항 등을 종합해 맞춤형 솔루션을 설계한다. 반면 초보 사용자는 “자동 설정” 기능에 의존할 수밖에 없지만, 그래도 최소한 어떤 알고리즘이 어떤 상황에서 쓰이는지 개념적으로 파악해 두면, 문제 해결 속도와 결과 품질을 훨씬 높일 수 있다.

추가적 예제와 확장 방향

본 장에서는 다양하고 대표적인 선형시스템 알고리즘의 연산 복잡도와 적용 사례들을 살펴보았다. 다음으로는, 실제 프로그램이나 오픈소스 라이브러리를 이용해 간단한 예제를 실행해 보고, 문제 크기와 조건을 바꿨을 때 계산 시간과 오차가 어떻게 변하는지 구체적으로 관찰할 것이다. 예컨대,

  • 희소계수행렬 문제에서 LU 분해와 CG+사후조건을 각각 적용해 성능 차이를 측정

  • 3D 탄성해석 문제를 domain decomposition + Krylov 방법으로 해결하며 코어 수를 늘려 병렬 스케일링 지표를 관찰

  • 밀집행렬 차원을 서서히 키워가며 직접해와 반복해, 혹은 저랭크 근사 기법의 실행 시간 비교

등을 실제로 수행해, 교과서적인 이론 복잡도와 현장에서 관찰되는 성능의 괴리가 어떻게 나타나는지 확인해볼 것이다. 이를 통해 알고리즘 선택이 어떤 맥락에서 이뤄지는지, 그리고 병렬 하드웨어와 수치 안정성, 에너지 효율 등 다면적 요소가 어떠한 식으로 복합 작용하는지를 좀 더 실감하게 될 것이다.

간단한 실험 예시: Python 기반 비교

이제까지 제시한 이론과 구현 이슈들을 실제로 실험해 볼 간단한 예시를 마련해 보자. 문제 규모가 그리 크지 않은 상황에서, Python의 NumPy/SciPy나 추가 라이브러리를 활용해 서로 다른 솔버를 적용하고, 연산 시간과 반복 횟수, 오차 등을 비교하는 방식이다. 아래에서는 세 가지 유형의 계수행렬을 간단히 생성하여, 실제 성능 차이를 체감해 보는 예시 코드를 제시한다. (실제 산업적 규모와는 거리가 있지만, 알고리즘 성질을 파악하는 학습용으로 유용하다.)

위 예시는 단순 무작위 행렬을 이용하기 때문에, 실제 물리 문제나 구조적 특성은 반영되지 않았다. 그러나 다음과 같은 시사점을 확인할 수 있다.

직접해(np.linalg.solvespsolve)는 밀집행렬일 때 $O(n^3)$ 연산이 발생하므로, $n$이 커지면 시간이 급격히 늘어난다. 하지만 반복해(CG)를 적용하려면 매 반복마다 행렬-벡터 곱셈이 $O(n^2)$을 소모하고, 조건수가 좋지 않으면 많은 반복($K$)이 필요할 수 있어 결과적으로 예상보다 느릴 수 있다.

희소행렬(test_solver_sparse)의 경우, spsolve는 내부적으로 희소 $LU$ 분해를 시도하므로, $n$이 매우 커지고 fill-in이 심하면 메모리·시간 부담이 커진다. CG는 희소 행렬-벡터 곱셈이 $O(n)$ 수준으로 수행 가능하나, 여전히 반복 횟수에 따라 실제 총연산이 정해지므로, 상황에 따라 어느 한쪽이 절대적으로 빠르다고 단정하기 어렵다.

더욱이 위 코드에서는 사후조건자나 병렬화, 혼합정밀도 기법 등을 전혀 고려하지 않았다. 실제로 이런 고급 기술들을 적용하면, 반복 횟수를 크게 줄이거나, GPU 병렬연산을 활용해 대규모 문제를 더 빠르게 풀 수 있다. 반대로 밀집행렬 분해도 GPU 가속(CuBLAS/CuSolver) 하에서 훨씬 짧은 시간에 끝날 수 있어, 단지 $O(n^3)$이란 이유로 포기하긴 어렵다. 결국 “실행 환경(하드웨어)”, “문제 구조”, “사후조건”, “병렬 알고리즘” 등이 전체 성능을 종합적으로 결정한다.

실제 규모 확장과 추가 실험 방향

위 예시는 $n$이 수천수만 선에서 머무르지만, 산업 현장이나 대규모 연구 문제는 그보다 훨씬 큰 차원($n$ 수백만수천만)으로 확장된다. 그러면 다음과 같은 상황 변화를 체감하게 된다.

  • 직접해의 메모리 한계: 분해 과정에서 fill-in이 발생하면 수십~수백 GB의 메모리를 요구하게 되어, 일반 워크스테이션으로는 해결이 불가능해진다.

  • 분산 컴퓨팅 및 병렬화: 여러 노드를 동원해도, 직접해는 피벗팅·삼각화 과정에서 빈번한 통신이 요구되어 병렬 효율이 낮아질 수 있다. 반면에 반복해는 (사후조건자 구성에 따라 다르지만) 주로 행렬-벡터 곱셈과 벡터 연산 중심이라 병렬화가 비교적 용이하고, 도메인 분할을 통한 지역 연산이 가능하다.

  • 조건수와 사후조건자 문제: $n$이 클수록, 불균등하거나 복잡한 구조의 행렬이 나타나기 쉽고, 반복해 솔버에서는 사후조건자 없이 수렴이 매우 느려진다. 이에 따라 인컴플리트 LU(incomplete LU), AMG(Algebraic Multigrid) 등 강력한 사후조건자가 실질적 필수가 된다.

  • GPU 가속 스케일링: GPU 여러 장을 노드별로 연결하거나, GPU 클러스터를 구성해 크고 복잡한 문제를 해결하는 시나리오가 늘고 있다. 이때는 데이터 분할, 통신(예: NCCL, NVLink, InfiniBand) 최적화, 혼합정밀도 사용 여부 등이 핵심 변수가 된다.

추가로, 위 코드에서 test_solver_dense, test_solver_sparse에 GPU 가속을 도입하려면, CuPy나 PyTorch, TensorFlow와 같은 GPU 컴퓨팅 라이브러리를 적절히 연결해 사용해야 한다. Scipy 자체는 CPU 연산을 중심으로 하므로, GPU 기반 선형연산을 위해서는 별도의 설정과 라이브러리가 필요할 수 있다.

통합적 관점과 탐구

본 장의 사례와 예시 코드는 “단순 비교”를 통해 어떤 솔버가 어느 정도로 빠른지 감을 잡는 데 도움이 되지만, 실제 문제는 훨씬 더 복잡한 변수가 많다. 수치 안정성, 사후조건자 선택, 하드웨어 스펙, 문제 구조, 혼합정밀도 활용 등 다양한 요소를 고려해야 하며, 때로는 반복해가 절대적으로 우세하나, 어떤 경우에는 직접해가 오히려 낫기도 한다.

이후 남은 분량에서는, 더 큰 스케일의 예제나 도메인 특화된 문제 설정을 통해, 위 코드와 유사한 실험을 좀 더 확장된 범위에서 시도해 볼 것이다. 또한 Python 이외에 C++나 Octave 예제로도 비슷한 접근 방식을 시도해 볼 수 있으며, MPI/OpenMP 병렬화를 약간 도입한 간단한 스켈레톤 코드를 살펴볼 예정이다. 이를 통해 “이론적 복잡도”와 “현실 계산 시간” 사이의 간극이 구체적으로 어떻게 나타나는지, 그리고 어떤 선택이 왜 좋은지 체감하는 기회를 갖도록 하겠다.

간단한 병렬 프로그래밍 예시 (MPI, OpenMP)

선형시스템 풀이를 대규모로 확장하려면, 멀티코어 혹은 멀티노드 환경에서 분산 메모리 병렬화나 공유 메모리 병렬화를 적용해야 한다. 오픈소스 C/C++/Fortran 진영에서는 MPI(Message Passing Interface)와 OpenMP(Open Multi-Processing)가 사실상 표준처럼 자리 잡고 있다. 여기서는 간단히 “병렬로 행렬-벡터 곱셈을 하는 스켈레톤 코드”를 예시로 소개하며, 실제로 반복해 알고리즘에서 이러한 병렬화를 어떻게 적용하는지 감을 잡아본다.

OpenMP 예시

OpenMP는 하나의 공유 메모리 내에서 여러 스레드를 생성하여, 반복문 단위로 일을 나누는 데 유용하다. 다음 코드는 C++에서 CSR 포맷 행렬 $A$와 벡터 $x$를 곱해 결과 벡터 $y$를 구하는 간단한 예시다. 실제로는 블록 분할이나 SIMD 벡터화 등의 최적화가 더해지겠지만, 구조만 간단히 살펴보자.

위 코드에서는 #pragma omp parallel for를 통해 각 행에 대한 반복을 스레드들이 자동으로 나누어 처리한다. 동시에 캐시 사용을 극대화하기 위해, CSR 배열을 연속적으로 순회할 때 메모리 접근 패턴이 선형에 가깝도록 설계한다. 실제로는 OpenMP 스케줄링 전략(static, dynamic 등)을 조절하거나, NUMA 환경에서 힌트를 주는 등 세부 튜닝이 필요하다. 이를 반복해 솔버(예: CG)의 한 번 반복에 그대로 적용하면, 큰 크기의 희소행렬-벡터 곱셈을 멀티코어 CPU에서 병렬로 처리할 수 있다.

MPI 예시

MPI 환경은 멀티노드 혹은 멀티소켓에서 분산 메모리를 활용할 때 쓴다. 각 프로세스가 자기 메모리에 부분행렬(submatrix)과 부분벡터(subvector)를 가지고 있으며, 행렬-벡터 곱셈 시 경계에서 다른 프로세스가 갖는 $x$의 일부를 받아와야 할 수도 있다. 단순화를 위해, 전체 행렬 $A$를 “행 단위”로 나눠서 $p$개의 프로세스가 각자 $A_i$와 $x$, $y$의 해당 구간을 처리한다고 해 보자.

실제 코드 구현은 훨씬 더 복잡해진다. 각 프로세스가 참조해야 하는 열 인덱스가 “자기가 맡지 않은 부분”에 속한다면, 그에 해당하는 $x$를 주고받아야 한다. 이를 위해 MPI 전송(점대점 또는 집단 통신)을 적절히 배치하거나, 다양한 최적화 기법을 쓴다. 반복해 솔버에서는 매 반복마다 이 과정을 반복하므로, 통신 오버헤드를 최소화하는 것이 핵심 과제가 된다.

반복해 솔버 병렬 스켈레톤

CG 같은 반복해 알고리즘에서 핵심 단계는 행렬-벡터 곱셈, 벡터 내적, 벡터 차(합) 연산, 알파·베타 스칼라와 곱셈 등이 반복되는 형태를 띤다. 병렬화 관점에서 다음처럼 정리할 수 있다.

초기화: $r^{(0)} = b - Ax^{(0)}$, $p^{(0)} = r^{(0)}$, $k=0$

반복:

  1. $\alpha_k = \frac{\langle r^{(k)}, r^{(k)}\rangle}{\langle p^{(k)}, Ap^{(k)}\rangle}$

  2. $x^{(k+1)} = x^{(k)} + \alpha_k p^{(k)}$

  3. $r^{(k+1)} = r^{(k)} - \alpha_k A p^{(k)}$

  4. $\beta_k = \frac{\langle r^{(k+1)}, r^{(k+1)}\rangle}{\langle r^{(k)}, r^{(k)}\rangle}$

  5. $p^{(k+1)} = r^{(k+1)} + \beta_k p^{(k)}$

여기서 $\langle \cdot, \cdot \rangle$는 벡터 내적이며, 이 내적 계산도 각 프로세스가 자기 파티션 내 요소들을 곱한 뒤, 최종적으로 MPI를 통해 합산(reduction)해야 한다. 행렬-벡터 곱셈 $A p^{(k)}$ 역시 위에서 본 MPI 분산 코드를 사용한다. 이런 식으로 구성하면, 대규모 클러스터에서 CG가 병렬로 동작하게 된다.

병렬 효율의 난관

단순히 $n$을 $p$개의 노드에 나누면 “이론적”으로는 연산이 $\frac{n}{p}$로 줄 것 같지만, 실제로는 다음과 같은 이유로 성능 저하가 발생한다.

각 스텝마다 MPI 집단 통신(내적 계산, 벡터 갱신 등)이 동반되며, 노드 수가 많아질수록 통신·동기화 비용이 증가 행렬 구조가 불균등하면 어떤 노드는 많은 비영 요소를 처리해야 하고, 다른 노드는 적은 요소를 처리하는 부하 불균형이 발생 NUMA 효과나 GPU 간 통신 대역폭 차이 등 하드웨어적인 병목이 추가 이로 인해, 실험 결과 병렬 효율(Parallel Efficiency)이 몇 천~만 코어까지는 어느 정도 괜찮다가, 그 이상 확대하면 10% 이하로 떨어지기도 한다. HPC 애플리케이션을 개발하는 입장에서는 이 효율을 극대화하기 위해 도메인 분할 자동화, 비동기화 통신, 오버랩(overlap) 기법, 고급 사후조건자 설계 등을 종합적으로 고려해야 한다.

---적 위치

OpenMP나 MPI를 통한 병렬화 예시는 “행렬-벡터 곱셈”이라는 가장 기본이 되는 연산을 어떻게 분산 처리하는지 보여준다. 여기에 다중격자(Multigrid)나 영역분할 기법을 결합하면, 각 부분 격자 혹은 부분 도메인에서의 연산을 병렬 처리하고, 경계 정보만 교환함으로써 반복 횟수와 통신량을 균형 잡는 더욱 복잡한 아키텍처를 설계할 수 있다. 실제 산업용·연구용 소프트웨어(예: PETSc, Trilinos, hypre, MUMPS, SuperLU_dist 등)는 이러한 접근을 라이브러리 차원에서 잘 추상화해 두어, 사용자가 직접 MPI 코드를 작성하지 않고도 병렬 알고리즘을 활용할 수 있게 해 준다.

현실적으로, 병렬 코드 스켈레톤을 직접 다루며 모든 것을 최적화하는 일은 매우 어렵다. 그러나 내부 작동 방식을 이해해 두면, 대규모 문제나 특별한 구조가 있는 문제에서 “왜 병렬 성능이 기대보다 낮은가”, “어떻게 사후조건자나 도메인 분할을 조절해야 하는가” 등을 좀 더 정확히 진단할 수 있다.

사례 연구 4: 유체역학(나비에-스토크스) 문제

유체역학에서 가장 중요한 방정식 중 하나가 바로 나비에-스토크스(Navier-Stokes) 방정식이다. 이를 점성(粘性)이 있는 비압축성 유체 흐름에 대해 간단히 표현하면,

{ρut+ρ(u)u=p+μ2u+f,u=0,\begin{cases} \rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u}\cdot\nabla)\mathbf{u} = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f}, \\ \nabla \cdot \mathbf{u} = 0, \end{cases}

와 같은 형태가 된다. 여기서

  • $\mathbf{u}(x,t)$는 유체 속도 벡터,

  • $p(x,t)$는 압력 스칼라,

  • $\rho$는 밀도,

  • $\mu$는 동점성계수(또는 점성계수),

  • $\mathbf{f}$는 외력(중력, 원심력 등)을 나타낸다. 비압축성이므로 $\nabla\cdot \mathbf{u} = 0$라는 추가 제약이 붙는다.

시간적·공간적 이산화와 선형시스템

유체역학 문제를 수치해로 풀려면, 시간 이산화(예: 전진오일러, 후진오일러, Crank-Nicolson 등)와 공간 이산화(Finite Difference, Finite Volume, Finite Element) 방식을 결합하여, 특정 시점(또는 $\Delta t$마다)에 유체의 속도와 압력을 계산해야 한다. 공간을 $n$개의 격자나 요소로 나누고, 속도와 압력을 벡터 형태로 나열하면, 매 시점마다 “속도장(velocity field)과 압력장”을 풀기 위한 방정식 세트가 거대한 비선형 시스템을 구성하게 된다.

이 비선형 시스템을 (뉴턴-유형의) 반복법이나 분리 알고리즘으로 풀 때, 각 스텝마다 큰 규모의 “선형 서브시스템”을 풀어야 하는 과정이 자주 등장한다. 예를 들어, 분리 알고리즘 중 가장 많이 알려진 것은 “압력-속도 분리법(Projection method)” 계열이다. 간단히 요약하면, 임시속도 $\mathbf{u}^*$를 점성항(및 외력)을 고려하여 먼저 계산하고, 이후 $\nabla p$ 항을 통해 비압축성($\nabla\cdot \mathbf{u}=0$) 조건을 만족하도록 보정(projection)하는 방식이다.

이 과정에서 핵심은 압력에 관한 푸아송(Poisson) 방정식이 등장한다는 점이다. 즉,

2p=ρΔtu\nabla^2 p = \frac{\rho}{\Delta t}\,\nabla \cdot \mathbf{u}^*

와 같은 형태(또는 비슷한 변형)를 해석 영역 전체에 대해 풀어야 한다. 이는 앞서 소개했던 포아송 문제와 본질적으로 유사하며, 2D 또는 3D 공간에서 대규모 희소행렬 선형시스템으로 귀결된다. 또한 점성항이나 대류항을 다루는 부분에서도 유사하게, 일정 형태의 큰 선형시스템을 반복해 풀어야 할 수 있다(특히 implicit한 시간적 이산화를 사용할 때).

알고리즘 선택 이슈

유체 흐름 문제는 다음과 같은 특성이 있어서, 선형시스템 알고리즘을 선택할 때 주의해야 한다.

  1. 매 시점(time step)마다 대규모 선형시스템 반복 유체 시뮬레이션은 보통 여러 스텝에 걸쳐 진행된다. 고정 $\Delta t$를 사용해 총 $N_t$번 스텝을 수행하거나, 혹은 적응형(Adaptive) 시간 스키마를 써서 점차 바뀌는 $\Delta t$로 여러 번 계산한다. 각 스텝마다 적어도 한 번 이상의 대규모 선형시스템(예: 압력 푸아송 방정식)을 풀어야 하므로, 전체 계산 비용은 “한 번의 선형시스템 풀이 비용 $\times$ 스텝 수 $N_t$”로 커진다.

  2. 행렬 구조의 지속적 변화 점성항은 비교적 안정적이며 대체로 대칭 양의 정부호 구조(특히 Stokes 한계에서)와 비슷하지만, 대류항(advective term)이나 비선형 효과는 시점마다 달라진 속도장에 의해 업데이트된다. 결과적으로, 한 스텝에서 구한 해를 다음 스텝에서 바로 사용할 때, 행렬이 완전히 바뀌거나(비선형 뉴턴 반복에서 접선행렬이 달라지거나) 일부가 바뀔 수 있다. 이로 인해, “직접해(LU 등)를 한 번 분해해놓고, 여러 번 우변 $\mathbf{b}$에 대해 재활용”하기가 쉽지 않은 경우가 많다. 오히려 매 스텝마다 새로이 행렬에 대응되는 분해를 구해야 하거나, 반복해 기법이 자연스럽게 선호된다.

  3. 대규모 3D 격자와 병렬화 필요성 실제 공학적 유체 문제(예: 항공기 날개 주변 유동, 해양/대기 시뮬레이션, 터보기계 내부 유동 등)는 3D 격자를 수백만, 수천만 요소 이상으로 구성해야 정확도를 확보할 수 있다. 이런 경우 직접해의 $O(n^3)$ 분해(또는 희소 행렬에서의 대규모 fill-in)는 메모리나 연산 측면에서 사실상 불가능해진다. 따라서 Krylov 반복해(예: GMRES, BiCGSTAB 등)나 다중격자(Geometric/Algebraic Multigrid), 도메인 분할법(FETI, BDDC, etc.)이 핵심적인 역할을 한다.

  4. 압력 푸아송 방정식의 조건수 문제 비압축성 조건을 강제하기 위한 푸아송 방정식은, 대규모로 갈수록 상태공간이 커지고 경계 조건(벽면, 유입/유출, 자유 경계 등)이 복잡해진다. 이때 행렬의 조건수가 매우 악화될 수 있다. 따라서 사후조건자(preconditioner)가 없으면 GMRES나 CG 유사 알고리즘은 반복 횟수가 크게 늘어나, 실질적 계산량이 $O(n^2)$ 이상이 될 수도 있다. 다중격자 기법은 푸아송 방정식에 매우 잘 맞는 대표적 접근이지만, 격자가 비정형이거나 복잡 지형을 다루면 기하학적 다중격자(Geometric Multigrid)를 구성하기 어렵고, 대신에 AMG(Algebraic Multigrid)를 쓰게 된다.

IPCS, SIMPLE, PISO 등 분할 알고리즘에서의 선형시스템

CFD(Computational Fluid Dynamics)에서 비압축성 나비에-스토크스를 푸는 알고리즘으로 IPCS(Incremental Pressure Correction Scheme), SIMPLE(Semi-Implicit Method for Pressure-Linked Equations), PISO(Pressure Implicit with Splitting of Operators) 등이 널리 알려져 있다. 이들 알고리즘 모두 “압력 방정식”과 “속도 방정식”을 어느 정도 분리하여 번갈아가며 푼다는 구조를 가진다. 예를 들어 SIMPLE 과정을 간략히 쓰면,

  1. 임시 속도장 $\mathbf{u}^*$ 계산(압력 항 추정)

  2. 압력 보정 방정식 $A_p p' = b$ 풀기

  3. 속도장 보정 $\mathbf{u} = \mathbf{u}^* + \dots$

이때 2번 단계에서 등장하는 $A_p$가 (푸아송 방정식과 유사한 형태로) 대칭 양의 정부호 희소행렬이 되거나, 경우에 따라선 비대칭·비정형 행렬이 될 수도 있지만, 대체로 Krylov 방법과 사후조건자가 잘 들어맞는다. $A_p$가 매 스텝마다(혹은 몇 스텝마다) 달라질 수 있으므로, 직접해를 매번 새로 분해하는 것보다는 반복해가 일반적으로 유리하다. 대규모 병렬 계산 시에는 도메인 분할 + Krylov + 다중격자 + 사후조건자라는 전체적 틀로 구성된다.

수치 안정성과 CFL 조건

유체역학 시뮬레이션은 시간적 안정성 조건(예: CFL(Courant-Friedrichs-Lewy) 조건)을 주의 깊게 관리해야 한다. 너무 큰 $\Delta t$를 쓰면 수치 발산이 일어나고, 너무 작은 $\Delta t$는 선형시스템 풀이 횟수를 과도하게 늘려 전체 계산 시간이 증가한다. 이와 같은 매 시점의 스텝 크기와, 각 스텝에서 선형시스템에 적용되는 알고리즘의 수렴 속도는 상호 연관되어 있다. 뉴턴-유형 비선형 반복 횟수, Krylov 반복 횟수, 그리고 시간 스텝 수 $N_t$가 서로 균형을 이뤄야 효율적인 시뮬레이션이 가능하다.

예를 들어 “implicit” 스키마를 채택하면, $\Delta t$를 좀 더 크게 설정해도 수치 발산을 억제할 수 있으나, 그만큼 각 시점에서 풀어야 하는 선형시스템이 더 복잡해져서 반복 횟수가 늘어나거나, 다중격자·사후조건자의 부하가 커진다. “explicit” 스키마면, 각 스텝의 선형시스템이 비교적 단순할 수 있으나, $\Delta t$가 엄격한 CFL 조건에 묶여서 $N_t$가 매우 커진다. 어떤 방법이 전체 계산 시간 관점에서 최선일지는 문제 규모, 유동 특성(레이놀즈 수, 경계조건 등), 하드웨어 성능 등에 따라 달라진다.

예시 구현과 확장

나비에-스토크스 문제 전체를 Python이나 C++로 예제 코드로 담으려면 매우 방대하겠지만, FEniCS, OpenFOAM, deal.II 같은 오픈소스 CFD/FEM 프레임워크를 살펴보면, 내부적으로 다음 단계를 거쳐 선형시스템을 해결한다.

  1. 격자/메시 생성 (3D 도메인 분할)

  2. 유한요소나 유한체적 방법으로 나비에-스토크스 잔차 항들을 이산화

  3. 분할 알고리즘(IPCS, SIMPLE 등)에 따라 “임시 속도장” 계산 → “압력 방정식” 형성 → “속도장 보정” 등의 순환

  4. 압력 방정식이나 속도 방정식을 푸는 단계에서, 대규모 희소행렬 $\mathbf{A}$와 벡터 $\mathbf{b}$를 만들고, 라이브러리(PETSc, Trilinos, hypre 등)에 제공

  5. Krylov 솔버(GMRES, BiCGSTAB 등) + 사후조건자(AMG, ILU, etc.)로 연산

  6. 해가 수렴하면 다음 시점(또는 다음 뉴턴 반복)으로 진행

병렬 환경에서 동작할 경우, 격자나 행렬 데이터가 노드별로 분산되고, MPI를 통해 경계의 자유도 정보를 교환하며, 다중격자나 영역분할 알고리즘으로 스케일링을 확보한다. 실제로 현대 슈퍼컴퓨터에서 수천~수만 코어를 써서 수억 개 격자 이상을 시뮬레이션하는 사례도 존재한다. 이러한 초대형 시뮬레이션에서는, 앞서 다룬 “병렬화 기법”과 “사후조건자 최적화”가 필수이며, 혼합정밀도(mixed precision)나 적응형 재메시(adaptive remeshing) 등 온갖 고급 기법이 활용된다.

시사점

유체역학 문제는 단순히 포아송 식 하나만 푸는 것보다 훨씬 복잡하고, 비선형 대류항과 압력-속도 간耦合(coupling)을 관리해야 하므로, 한 번의 선형시스템 풀이가 아니라 “반복적, 다단계” 풀이 절차가 필연적으로 등장한다. 따라서 직접해의 반복 사용은 매우 비효율적일 수 있고, 대규모 문제에서는 “반복해 + 사후조건자 + 다중격자 + 병렬화” 구조가 주류를 이룬다. 특히 압력 방정식에서의 푸아송 타입 문제를 빠르게 해결하는 다중격자/AMG 방식은 유체 흐름 시뮬레이션에서 사실상 표준이다.

여기에 시간 스키마 선택, 안정성 조건, 뉴턴-유형 비선형 수렴 등 다양한 요소가 얽히므로, “수치해석 알고리즘”과 “PDE 알고리즘”이 긴밀히 교차한다. 결국 이 모든 과정을 구성하는 소프트웨어적인 방대한 인프라가 현대 CFD 및 HPC 응용의 기반이 된다.

사례 연구 5: 전자기(맥스웰) 해석 문제

유체역학 외에도 전자기장을 다루는 맥스웰(Maxwell) 방정식 해석에서도 큰 규모의 선형시스템이 반복해서 등장한다. 정전기, 정자기 문제는 비교적 단순한 포아송/라플라시안 방정식과 유사하지만, 전자파(고주파) 해석이나 안테나·회로 해석을 위해 맥스웰 방정식을 전체적으로 풀어야 할 때는 보다 복잡한 구조의 편미분방정식을 다룬다.

전자파 해석에서 시간영역(Time-Domain)이나 주파수영역(Frequency-Domain) 접근을 택하든, 공간 이산화(FEM, FDTD, FIT 등)를 택하든, 결국 대규모 연립방정식을 풀어야 한다. 대표적으로 주파수영역 FEM 방식에서, 편미분방정식을 이산화해 다음과 같은 선형시스템이 얻어지곤 한다.

KE=F\mathbf{K}\mathbf{E} = \mathbf{F}

여기서 $\mathbf{E}$는 전기장(또는 자기장) 분포를 의미하는 벡터, $\mathbf{K}$는 재료 물성(유전율, 투자율)과 주파수, 경계조건 등에 의해 결정되는 계수행렬, $\mathbf{F}$는 외부 소스(전류원, 외부 자극)에 해당한다. 문제 규모가 커지면, $n$이 수백만~수천만 차원에 달하기도 하며, $\mathbf{K}$가 복소수를 포함하거나 대칭이 아니거나(또는 준-에르미트 행렬 형태) 다양한 특성을 지니게 될 수 있다.

주파수 영역 vs. 시간 영역

주파수영역 해석에서는 Helmholtz 방정식 형태(파동 방정식을 변형한 것)를 풀어야 하며, 복소 계수행렬이 형성된다. 특히 유전체 경계나 금속 경계 처리, 경계조건(Perfect Electric Conductor, Perfectly Matched Layer 등)에 따라 행렬 구조가 달라진다. 이때 직접해($LU$ 등)는 $O(n^3)$ 연산과 큰 메모리를 요구하기에, 대규모 문제에서는 반복해(GMRES, BiCGSTAB 등)가 자주 사용된다. 그러나 이 행렬이 편의 양의정부호(symmetric positive definite)인 것은 아니므로, CG 대신 다양한 Krylov 방법이 고려되며, 사후조건자를 잘 설계해야 반복 횟수를 줄일 수 있다.

시간영역 해석(FDTD 등)에서는 매 시간 스텝마다 맥스웰 방정식을 부분적으로 선형화하거나, 완전히 직접 적분하는 기법을 사용할 수도 있다. 예컨대 Yee 격자를 사용하는 FDTD(Finite-Difference Time-Domain) 방식은 대체로 선형시스템 풀이 없이 “직접 갱신(Explicit Update)”을 반복하지만, 경계나 물성에 따라 세부적으로 암시적(Implicit) 계산이 필요한 부분이 있으면, 그 구역에서 대규모 연립방정식을 풀어야 한다. 이때에도 반복해 솔버가 활용되는 경우가 많다.

행렬 구조와 고주파 특성

고주파 전자기장 문제에서는 파장이 짧아지고, 도메인 내부에 수많은 자유도가 생기면서 행렬 차원이 매우 커진다. 동시에, 고유값 또는 스펙트럼 분포가 복잡해지므로, 단순한 사후조건자가 잘 작동하지 않을 수 있다. 도전성 경계(도체), 저손실 유전체 등 다양한 물리적 특성이 혼재하면, $\mathbf{K}$가 비대칭 또는 준-에르미트(semi-Hermitian) 형태가 되기도 한다.

이 경우 다중격자(Multigrid) 적용이 쉽지 않을 수 있다. 포아송형 문제처럼 단순 라플라시안 연산자가 아니라, 파동성분이 강하게 나타나는 운영자(operator)를 다루게 되기 때문이다. 그럼에도 불구하고 다양한 변형(예: 적응형 다중격자, 영역분할, 사후조건) 기법이 연구되며, 파동방정식 특성을 고려한 고급 사후조건자가 제안되기도 한다.

또한 주파수영역 해석에선, 행렬-벡터 곱셈의 복소 연산이 표준이므로, 실제 구현에서는 실수 배열 두 개(실수부, 허수부) 또는 복소수 타입을 다루어야 한다. 병렬화 역시 이 복소수 연산을 처리하도록 확장되어야 하며, MPI나 OpenMP에서 복소수 데이터 타입을 표준으로 지원하거나, 사용자가 구조를 설계해야 한다.

FFT 기반 솔버와 하이브리드 방법

맥스웰 방정식(특히 균일 매질이나 주기 구조)에서 빠른 푸리에 변환(FFT) 기반 기법을 사용하면, 복잡도와 메모리 요구량을 획기적으로 낮출 수 있는 경우도 있다. 예컨대 주기구조를 가진 광학 시뮬레이션에서는 광자결정(Photonic Crystal) 해석에 FFT를 활용해, 행렬 생성 자체를 회피하거나, 근사 연산으로 교환할 수 있다. 그러나 일반적인 비균질/복합 구조에서는 FFT만으로는 한계가 크고, 결국 FEM/FDTD 등을 조합한 하이브리드 방법을 모색한다.

이러한 하이브리드 방법은 시스템의 일부분(예: 균일 영역)에선 FFT 혹은 특수 풀이를 적용하고, 나머지 복잡한 부분(경계나 국부 구조)에선 FEM 또는 FD 기법으로 선형시스템을 형성해 풀도록 분할한다. 그러면 전체 연산량과 메모리 부담을 줄일 수 있으나, 각 부분을 어떻게 연결할지(계면 조건, 모드 매칭 등) 설계가 복잡해진다. 병렬화와 사후조건 기법도 이 분할 구조에 맞춰 최적화가 필요하다.

초대형 전자기 문제의 HPC 적용

실제 전자기 해석 소프트웨어(예: CST Studio, ANSYS HFSS, COMSOL RF 모듈 등) 내부에는 다양한 선형시스템 솔버가 내장되어 있다. 예를 들어 ANSYS HFSS는 FEM 기반 주파수영역 해석을 할 때, 유저가 직접해(LU) 혹은 반복해(Krylov, 다중격자, 영역분할) 등을 선택할 수 있도록 옵션을 제공한다. 유전체가 단순하고 메모리가 충분하다면 직접해가 유리할 수도 있지만, 큰 스케일에서 반복해가 사실상 표준이다.

병렬 HPC 환경에서 수백~수천 코어를 동원해 안테나나 레이더 단면(RCS) 해석, 고주파 필터나 광대역 해석 등을 진행하는 사례가 늘고 있다. 이런 초대형 시뮬레이션에서는:

  • 시스템 행렬이 복소수 대칭(또는 준-대칭)

  • 큰 주파수일수록 메시가 고밀도

  • 비균질 물성이 많으면 조건수 악화

  • 다중격자 적용이 어려우면, ILU(불완전 LU)나 블록 사후조건자, 혹은 영역분할 기법이 필수

따라서 연산 복잡도가 단순히 $O(n^3)$나 $O(Kn^2)$로만 나타나지 않고, 복소수 연산의 상수 계수나 통신 비용, 사후조건자의 성능 등이 복합적으로 작용한다. 더불어 파동문제 특유의 공진(resonance)이나 모드 구조, 경계 흡수 조건(PML) 등도 수렴 특성에 영향을 준다. 문제 세팅에 따라서는 반복 횟수가 기하급수적으로 늘어날 수도 있으므로, 물리 해석 단계에서부터 “어떤 주파수 구간을 어떻게 스위핑(sweeping)할 것인지” 전략적으로 접근해야 한다.

확장된 관점

전자기 문제는 물리적 양상이 다양하고 수치적으로도 난이도가 높기 때문에, 선형시스템 알고리즘 연구에서 끊임없이 새로운 방법들이 제안되는 영역 중 하나다. 다중스케일(multiscale) 구조, 준결정(quasicrystal) 패턴, 초재료(metamaterials) 같은 복잡한 물리계가 추가되면, 행렬 조건수가 더욱 나빠지고, 전산량도 거대해진다. 따라서 “반복해 + 병렬화 + 고급 사후조건”은 필수이고, 여기에 FFT, FMM(Fast Multipole Method), 랜덤화 기법 등을 혼합해 연산을 줄이는 아이디어가 계속 발전 중이다.

이렇듯 전자기 해석에서도 “선형시스템 풀이 알고리즘의 연산 복잡도”는 매우 중요한 연구 주제이며, 산업 현장에서 대규모 안테나·레이다·광학 시뮬레이션을 수행할 때 실무적 영향력이 크다. 앞서 살펴본 CFD(유체역학), 탄성해석, 포아송 방정식 등과 본질적으로 같은 “선형시스템에 대한 반복해와 직접해의 선택” 문제이지만, 복소계수, 비대칭성, 파동성(wave phenomenon) 등으로 인해 또 다른 관점의 난제가 추가되는 셈이다.

대규모 회귀·분류 문제와 선형시스템

기계학습(머신러닝) 분야에서도 선형시스템을 대규모로 풀어야 하는 상황이 자주 등장한다. 전통적인 통계학의 최소제곱회귀(OLS)나 로지스틱회귀, 서포트벡터머신(SVM) 등을 예로 들면, 문제 크기($n$ 혹은 특성 차원 $p$)가 커졌을 때 $O(n^3)$ 또는 $O(n^2)$ 연산이 쉽지 않아진다. 또한 빅데이터 시나리오에서 샘플 수가 수백만 이상이거나, 특징 차원 $p$가 수십만~수백만에 이를 수 있어, 표준적인 직접해(행렬 분해) 방식으로는 계산 불가능한 수준이 될 수 있다.

로지스틱회귀의 경우를 간단히 보면, $n$개의 샘플 $\mathbf{x}_i \in \mathbb{R}^p$와 레이블 $y_i \in {0,1}$가 주어졌다고 할 때, 음의 로그우도(negative log-likelihood)를 최소화하는 문제는

minβRp(i=1n[yilnσ(βxi)+(1yi)ln(1σ(βxi))]),\min_{\beta \in \mathbb{R}^p} \, \bigl( -\sum_{i=1}^n \bigl[y_i \ln \sigma(\beta^\top \mathbf{x}_i) + (1 - y_i)\ln \bigl(1-\sigma(\beta^\top \mathbf{x}_i)\bigr)\bigr] \bigr),

여기서 $\sigma(z) = \frac{1}{1 + e^{-z}}$는 시그모이드 함수이다. 이를 뉴턴(또는 준뉴턴) 방식으로 풀 때, 매 반복에서 근사적인 2차 서브문제를 풀어야 하며, 결국 선형시스템

(XWX)Δβ=X(yp)(\mathbf{X}^\top \mathbf{W} \mathbf{X})\, \Delta \beta = \mathbf{X}^\top \mathbf{(y - p)}

같은 형태(피셔 스코어 방정식 등)가 생긴다. 여기서 $\mathbf{W}$는 가중 행렬(대각 행렬), $\mathbf{p}$는 예측 확률 벡터, $\Delta \beta$는 업데이트량 등이며, $\mathbf{X}\in\mathbb{R}^{n\times p}$가 크면(또는 희소 행렬이 아니면) 이를 매 스텝 해석적으로 분해하기가 부담스러울 수 있다. 반복해 알고리즘으로 근사해나가거나, 확률적 기법(SGD, 미니배치 등)을 사용해 전체 선형시스템을 직접 풀지 않고 점진적으로 해를 구하는 접근이 확립되어 있다.

확률적 기법과 선형시스템

대규모 데이터를 처리할 때 확률적 경사하강법(SGD)이나 미니배치(mini-batch) 방법을 쓰는 이유는, 매 스텝에서 전체 $\mathbf{X}$를 다 보지 않아도 근사적인 방향으로 매개변수(예: $\beta$)를 업데이트할 수 있기 때문이다. 이는 “정확한 선형시스템을 풀어 매번 최적해로 가는” 방식이 아니라, “부분 표본이나 확률적 근사”로 수렴을 이끄는 방법이다.

반면, 어떠한 상황에서는 실제로 선형시스템을 풀어야 하는 정밀 해법이 요구된다. 예를 들어, 커널 서포트벡터머신에서 듀얼 문제를 풀 때, 듀얼 변수에 대해

Qα=1\mathbf{Q}\alpha = \mathbf{1}

꼴(또는 유사 형태)의 문제를 형성할 수 있는데, $\mathbf{Q} \in \mathbb{R}^{n\times n}$가 밀집행렬(커널행렬)이 되면 $O(n^3)$ 연산을 피하기 어렵다. 이 경우, 근사 커널 기법(Nystrom 근사, 랜덤 피처)이나 반복해 방법으로 부분적으로만 접근할 수 있다. 반복해를 쓰려면 한 번의 곱셈 $\mathbf{Q}\alpha$가 $O(n^2)$이고, 필요 반복 횟수가 $K$라면 $O(K n^2)$이 되므로, $K$가 작다면 직접해($O(n^3)$)보다 유리할 수도 있고, $K$가 크면 오히려 더 나빠질 수도 있다.

따라서 커널 기반 기법에서 “직접해 vs. 반복해 vs. 근사”는 매우 중요한 의사결정 포인트가 된다. $n$이 수천~수만 정도라면, GPU 가속으로 $O(n^3)$ 분해를 감당할 수 있는 경우도 있지만, 그 이상이면 거의 불가능하고, 반복해나 근사로 전환하는 경향이 크다.

희소성, 정규화, 그리고 사후조건자

머신러닝에서 라소(Lasso)나 엘라스틱넷(Elastic Net) 같은 정규화 모델을 쓰면, 해가 희소해질 수 있지만 계수행렬 $\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}$ 자체가 희소해지는 건 아니다(데이터 행렬 $\mathbf{X}$가 희소 구조를 갖는 경우는 예외). 다만, 이러한 정규화 항은 조건수를 개선하는 효과가 있어서, 반복해 알고리즘이 더 안정적으로 수렴할 수 있다. 예를 들어 릿지 회귀(ridge regression)에서는

(XX+λI)β=Xy(\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I})\,\beta = \mathbf{X}^\top \mathbf{y}

이 생성되는데, $\lambda > 0$가 크면 조건수가 낮아져 빠른 수렴이 가능해진다. 이때 “사후조건자 + 반복해” 조합이 강력한 방법이 될 수 있고, 실제 산업체나 대규모 데이터 분석에서 유용하다.

사후조건자 설계는 머신러닝 맥락에서도 여전히 통한다. 예를 들어, 데이터가 특정 클러스터 구조나 블록 구조를 갖는다면, $\mathbf{X}^\top \mathbf{X} + \lambda \mathbf{I}$를 블록 대각 형태로 근사하거나, Incomplete Cholesky 분해를 적용해 사후조건자로 활용할 수 있다. 문제 크기가 정말 크다면, 랜덤화 기법이나 스케치(sketch) 기법을 통해 더 작은 차원의 서브문제를 사후조건자로 삼는 연구도 진행 중이다.

분산 학습과 병렬 선형연산

빅데이터 시대에는 데이터가 여러 노드에 걸쳐 분산되어 저장되고, 학습 과정도 분산된 환경에서 이뤄진다. 이때 선형시스템을 직접 풀려면, 각 노드가 자기 데이터 부분을 이용해 부분 행렬을 형성하고, 그 결과를 MPI나 RPC(Remote Procedure Call) 기반으로 합치는 과정을 거쳐야 한다. 소위 “Parameter Server” 구조나 “Ring Allreduce” 방식 등을 통해, 그라디언트나 파라미터를 교환하기도 한다.

분산 환경에서 CG나 GMRES 같은 반복해 알고리즘을 사용하려면, 매 반복마다 네트워크 통신이 발생할 수 있다. 특히 전역 내적이나 벡터 합산 과정에서 모든 노드가 동기화해야 하므로, 노드가 수백~수천 개로 늘어나면 통신 병목이 발생하기 쉽다. 이 문제를 완화하기 위해, 비동기(Asynchronous) 분산 알고리즘이나 로컬 업데이트 후 주기적으로만 글로벌 동기화를 하는 알고리즘이 탐색된다.

머신러닝에서 자주 사용하는 SGD, 미니배치 SGD, Adam, AdaGrad 등은 사실상 선형시스템을 풀지 않고, 확률적 1차 방식으로 최적점을 근사한다. 그럼에도 불구하고, “최종적인 정밀 해”를 구하거나, “작은 크기의 서브문제”를 강건하게 풀어야 하는 상황에서는 여전히 선형시스템 해결이 필요하다. 예컨대 제한볼츠만머신(RBM), 가우시안 공정(Gaussian Process) 기반 모델, 커널 SVM, 혹은 Bayesian Inference에서의 후방분포 계산 등에선 정밀한 선형대수 해법이 활용된다.

예시 코드 (Python) : 랜덤화된 스케치

아래 예시는, 밀집행렬 $\mathbf{X}\in\mathbb{R}^{n\times p}$가 너무 커서 직접해가 부담스러운 상황을 가정하고, “랜덤 스케치” 기법으로 차원을 축소한 뒤에 선형시스템(릿지 회귀 형태)을 풀어보는 개념적 스케치 코드다.

이 코드는 $n$이 매우 크고 $p$도 큰 경우, 원래의 $O(\min(n^3, np^2))$ 직접해를 대신해, “$n$차원 데이터를 무작위 스케치로 sketch_size 차원으로 축소”한 다음, 축소된 문제에서 $O(p^3)$ 또는 $O(\texttt{sketch_size} \cdot p^2)$ 연산으로 해를 구한다. 정확도는 근사적으로 유지될 가능성이 높으나, 완전히 정확한 해는 아니다. 하지만 빅데이터 회귀나 정규화 모델에서는 이런 근사가 실무적으로 종종 사용된다.

만일 $p$가 수십만에 달해도 sketch_size가 충분히 작으면 계산이 가능하며, 이때 반복해 기법이나 사후조건자를 더 결합할 수도 있다. 결국 “대규모 선형시스템을 정확히 풀기 어려울 때, 문제 자체를 근사(차원 축소)하거나, 확률적 경사방법으로 접근”하는 기계학습 패러다임이 지배적으로 자리 잡고 있다.

시사점

머신러닝·딥러닝 관점에서 선형시스템 풀이에 대한 관심은, 빅데이터 시대 이후에도 여전히 중요하게 남아 있다. 확률적 1차 방법이 급부상한 지금도, 대규모 가우시안 프로세스나 커널 기법, 정교한 베이지안 방법을 다룰 때는 고전적 선형대수 해법이 결정적 역할을 하기 때문이다. 특히 실시간 예측보다는 오프라인 학습에서 높은 정밀도가 요구되는 상황, 또는 저차원의 부분문제를 정확히 풀어야 하는 상황 등에선, “어떤 알고리즘으로 어떻게 풀어야 효율적이고 안정적인지”가 핵심 질문이 된다.

계산환경도 GPU, 분산 클러스터, 혼합정밀도 등으로 빠르게 변화하고 있어, “직접해 vs. 반복해”의 구도가 단순히 빅오 표기로만 해결되지 않는다. 데이터 행렬의 구조(희소도, 랭크분포), 정규화나 사후조건자 사용 여부, 그리고 확률적 근사 기법의 적용 여부 등을 종합해, 궁극적으로 실제 계산 시간을 최소화하고 목적 함숫값을 만족하는 전략을 세우는 일이 중요해진다.

양자 컴퓨팅과 선형시스템

최근 고성능 계산(HPC)의 새로운 패러다임으로 떠오른 양자 컴퓨팅(Quantum Computing)에서도, 특정 유형의 선형시스템 풀이 문제를 가속화할 수 있는 가능성이 논의되고 있다. 대표적으로 Harrow-Hassidim-Lloyd(HHL) 알고리즘이 알려져 있는데, 이는 선형시스템 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 양자 상태로 인코딩한 뒤, 양자 회로 연산을 통해 $\mathbf{x}$를 근사적으로 구하는 방법이다. 이상적인 양자 게이트와 충분히 큰 양자 얽힘(entanglement)을 가정한다면, 고전적으로 $O(n^p)$의 시간이 걸리는 선형시스템 해결을 $\widetilde{O}(\log n)$ 정도로 단축할 수 있다는 이론적 결과가 소개되어 주목받았다.

하지만 실제로 이를 구현하는 것은 여러 가지 난제가 따른다. 양자 상태에 $\mathbf{A}$와 $\mathbf{b}$를 효율적으로 로드(loading)하고, 결과 벡터 $\mathbf{x}$를 다시 고전적 정보로 추출(read-out)하는 과정이 $O(n)$ 이상의 복잡도를 필요로 하면, 결국 전체 알고리즘 이점이 크게 줄어들 수 있다. 또한 HHL 알고리즘은 행렬이 좋은 조건수와 희소 구조(또는 적절한 블록 구조)를 가져야 양자 이점을 제대로 활용할 수 있다는 전제조건이 뒤따른다. 따라서 현재로서는 양자 컴퓨팅을 이용해 대규모 선형시스템을 실제로 고전적 알고리즘보다 빠르게 푸는 사례가 실용적으로 구현된 바는 매우 적다.

그럼에도 불구하고 양자 알고리즘 이론은 꾸준히 발전 중이며, 단순히 HHL에 국한되지 않고, Krylov 서브스페이스 방식과 양자 서킷을 접목하거나, 반복해 알고리즘의 특정 단계를 양자화(quantization)하여 병목 부분을 가속하려는 시도도 있다. 양자 하드웨어가 점점 발전하여 노이즈를 줄이고 큐비트 수를 늘려갈수록, 선형시스템 풀이에도 일부 영향이 있을 것이라는 기대가 존재한다.

자동 미분과 선형시스템

머신러닝, 최적화, PDE 역문제 등에서 자동 미분(Automatic Differentiation, AD) 기법이 널리 활용되면서, 역전파(backpropagation) 단계에서 암묵적으로 선형시스템을 푸는 상황이 간혹 발생한다. 예를 들어, 네트워크 구조가 복잡하거나, 물리기반 손실함수(Physics-Informed Neural Networks)에서 PDE 제약을 만족하는 항이 들어갈 경우, 부분적으로 선형시스템 풀이가 은닉된 형태로 등장할 수 있다.

특히 역문제(inverse problem)나 최적제어(optimal control)에서, “Adjoint State” 방정식을 풀어야 하는 상황에서는, 매 스텝마다 대규모 선형시스템이나 선형 연산을 역으로 추적해야 한다. 이 과정은 수치해석적 견지에서 보면 “사후조건 Krylov 솔버를 내부에 내장한 역전파”와 같은 형태를 띠게 된다. PyTorch나 TensorFlow 같은 딥러닝 프레임워크에서도 이러한 고급 요구사항을 만족하기 위해, 부분적으로는 SciPy/PETSc/Trilinos 등 외부 라이브러리와 연동하거나, 자체적인 선형알고리즘 루틴을 탑재하는 방향으로 확장되고 있다.

이러한 “AD + 선형시스템” 흐름은 최적화 문제를 다루는 자동미분-기반 라이브러리들(JAX, Julia의 Zygote, Ceres Solver, TAO 등)에서도 나타난다. 대규모 Hessian이나 Jacobian을 명시적으로 만들지 않고 곱셈 연산(또는 역곱)을 자동미분으로 계산하여, Krylov 반복에서 “행렬-벡터 곱”만 구현하는 식으로 접근하기도 한다. 이렇게 하면 행렬을 “표 explicitly” 만들 필요가 없어 메모리 사용량을 크게 줄이면서, 반복해 알고리즘으로 근사 해를 구할 수 있다.

멀티피직스(multi-physics) 문제와 결합 솔버

선형시스템 해석은 독립된 한 가지 물리 분야에만 머무르지 않는 경우가 많다. 예컨대 유체-구조 연성(fluid-structure interaction, FSI) 문제는 유체역학 문제와 탄성 문제를 동시 해석해야 하며, 열-유동-화학 반응이 결합된 다중물리(multi-physics) 현상도 흔하다. 이런 멀티피직스 문제는, 서로 다른 PDE 모형이 상호 연결되면서 훨씬 복잡한 연립방정식 시스템을 형성한다.

이들 문제를 전역 행렬 단위로 한꺼번에 풀면 초거대 차원의 선형시스템이 생기고, 각 하위 물리 모듈에서 해를 반복 교환하는 loose coupling 방식으로 접근할 수도 있다. 전자의 경우, 행렬 구조가 블록 대각블록(block-diagonal-block) 형태로 나타날 가능성이 커, “블록 사후조건자”나 “블록 Krylov 방법” 같은 기법이 연구된다. 후자의 경우, 각 물리 모듈이 독립적으로 자율적 반복해를 수행하고, 경계나 인터페이스에서만 서로 조건을 맞춰주므로, 도메인 분할과 유사한 병렬화 프레임워크를 적용하기도 한다.

예컨대, 열-유동 문제에서 “유동 방정식”과 “에너지 방정식”을 한꺼번에 풀면, 속도·압력·온도에 대한 거대한 벡터가 형성되고, 해당 연립방정식을 블록으로 나누어 보며 “온도 블록”과 “속도-압력 블록” 간 사후조건을 어떻게 설계할지 고민해야 한다. 이런 블록 사후조건자가 제대로 동작하면, 전체 반복 횟수가 크게 줄어든다. 다만 구현의 복잡도가 상승하기 때문에, 실제 산업용 소프트웨어에선 문제 규모와 성질에 따라 해법을 달리 선택한다.

동적 적응형 재메시(adaptive remeshing)와 선형시스템 갱신

PDE 해석에서 국소적으로 오차가 큰 부분만 격자를 세밀화(refinement)하거나, 변화가 작은 부분을 반대로 격자를 성글게(coarsening) 하는 적응형 재메시(adaptive mesh refinement) 기법을 사용하는 경우가 많다. 이럴 때, 시간 진행이나 뉴턴 반복 단계 중에 격자가 바뀌면서 계수행렬 $\mathbf{A}$도 함께 바뀌므로, 이미 계산해둔 분해나 사후조건자를 재활용하기 어렵게 된다. 매번 새로운 크기와 구조의 $\mathbf{A}$가 형성될 수 있기 때문이다.

반복해 알고리즘에서는, 사후조건자를 매번 새로 구성하지 않고 “이전 단계에서의 사후조건자를 변형해 재활용”하는 기법이 연구되고 있다. 이를테면, AMG 계층(코스 격자와 연산자)을 국소적으로만 업데이트하거나, incomplete LU에서 새로 생긴 fill-in에 대해서만 보강 정보를 추가하는 등의 방법을 시도한다. 이를 통해 재메시가 잦아도, 완전히 새로 사후조건자를 만드는 데 드는 오버헤드를 줄여, 전체 계산 시간을 단축한다.

이처럼 적응형 기법과 결합했을 때에도, “직접해를 매번 새로 분해하기”는 실질적으로 부담이 매우 큰 편이므로, 반복해 방식이 훨씬 더 보편화되어 있다. 다만, 적응형으로 인해 매 스텝에서 계수행렬이 크게 변형되는 상황이라면, 일반적인 사후조건자 설계가 어렵고, 수렴 보장이 흔들릴 수 있다는 점을 주의해야 한다. 결국, 문제 특성별로 기법을 세심하게 조합해야 최적의 성능을 얻을 수 있다.

후속 논의: 사례 종합과 추가 예제

이와 같이, 선형시스템 해결은 수많은 최첨단 분야(양자 계산, 자동 미분, 다중물리, 적응형 재메시 등)와 접목되며 계속 확장되고 있다. 본 장의 앞부분에서 다룬 전통적 직접해와 반복해의 연산 복잡도 비교, 희소성·사후조건자·병렬화·근사 기법 등은 이런 복합 시나리오에서도 핵심 개념으로 작용한다. 각각의 응용분야는 자체적으로 특화된 접근법과 소프트웨어 프레임워크를 발전시켜 왔지만, 그 근저에는 “선형시스템을 어떻게 효율적이고 안정적으로 풀 것인가”라는 공통 과제가 놓여 있다.

본 장의 남은 부분에서는 앞서 소개된 다양한 응용 사례를 종합하고, 코드 예시를 조금 더 구체화하여, 실제로 가우스 소거법·LU 분해·Krylov 방법·다중격자 등을 실행할 때 나오는 실제 계산 시간과 메모리 사용량, 그리고 오차 양상을 비교하는 표나 그래프를 예시로 보이려 한다. 또한 간단한 병렬 코드나 사후조건자 실험을 수행해, 독자들이 직접 시도해볼 수 있는 템플릿을 제공하고자 한다.

고급 사후조건자(Preconditioner) 설계

반복해 알고리즘의 성능을 획기적으로 높이는 핵심 도구 중 하나는 사후조건자(preconditioner)다. 사후조건자가 행렬의 조건수(condition number)를 획기적으로 낮춰주면, 반복 횟수($K$)가 현저히 줄어드는 효과가 있어, $O(K n^2)$ 복잡도를 큰 폭으로 감소시킨다. 하지만 사후조건자를 구성하는 과정 자체가 추가적인 비용(메모리, 연산)을 유발하므로, “얼마나 빨리, 작게 만들 수 있는가”가 결정적 이슈가 된다.

가장 간단한 예로는 대각 사후조건자(diagonal preconditioner)가 있다. $\mathbf{M}$을 $\mathbf{A}$의 대각성분만 모아놓은 행렬(대각행렬)로 잡으면, $\mathbf{M}^{-1}$ 계산은 쉽지만, 복잡한 문제에서 이 사후조건자는 효과가 미미하다. Incomplete LU(ILU)나 Incomplete Cholesky(IC) 같은 기법은 직접 $LU$, $Cholesky$ 분해를 단순화하여 일부 fill-in만 허용하고 나머지는 버리는 방식으로 $\mathbf{M}\approx\mathbf{A}$를 구성한다. fill-in의 정도를 제어하는 드롭(drop) 파라미터나 레벨(level)을 설정함으로써, 사후조건자의 정밀도와 계산 비용을 절충한다.

다중격자(Multigrid)나 영역분할(Domain Decomposition) 기법도 본질적으로는 사후조건자의 일종으로 해석할 수 있다. Krylov 반복($CG$, $GMRES$ 등)의 내적 연산마다 $\mathbf{M}^{-1}$을 곱해주면, 사실상 다중격자 순환(V-cycle, W-cycle 등)을 한 번 호출하는 구조와 유사하게 작동하기 때문이다. 이를 구현 차원에서 보면, Krylov 반복 루프 안에 “사후조건(Preconditioner Application)” 함수를 삽입해 두고, 여기서 AMG(Algebraic Multigrid) 한 번의 스무딩(smoothing)과 코스 격자 해결(coarse grid solve)을 수행하는 방식을 취한다.

사후조건자를 만드는 데 드는 시간 및 메모리를 최소화하기 위해, 최근에는 블록(작은 부분행렬) 단위로만 ILU를 수행하거나, 여러 단계에 걸쳐 부분분해를 수행하는 기법이 제안된다. 또, 행렬 구조가 주기적이거나, 랭크가 특정 영역에서 작게 유지되는 경우(예: H-행렬, 블록 저랭크 구조), 해당 블록들을 저랭크 근사해 사후조건자로 쓰는 시도가 이루어진다. 이런 경우, 큰 스케일에서 $O(n\log n)$ 또는 $O(n)$에 가까운 곱셈과 사후조건 연산이 가능할 수도 있으나, 구현 복잡도와 안정성 유지가 상당히 까다롭다.

고급 주제: Incomplete Factorization(ILU) 예시

대표적으로 Incomplete LU(ILU) 분해는 직교분해(QR)나 완전 $LU$ 분해보다 훨씬 적은 메모리로 대규모 희소행렬의 사후조건자를 만들 수 있다는 장점이 있다. 간단한 ILU(0)는 “대각선 아래·위의 원래 위치에 있던 비영 요소만 유지하고, 새로운 fill-in을 허용하지 않는다”는 방식이고, ILU(k)는 $k$ 레벨까지 fill-in을 허용한다. 드롭 전략이 추가되면 ILU($k$)-드롭 등의 변형이 나온다.

예를 들면 $n \times n$ 희소행렬 $\mathbf{A}$가 주어졌을 때, ILU(0)는 다음과 같은 단계로 구성할 수 있다.

  1. 가우스 소거 과정을 수행하되, 원래 $\mathbf{A}$에서 0이었던 위치는 분해 과정에서 새로 fill-in이 생겨도 그냥 버린다.

  2. 유지된 위치에 대해서만 $L$, $U$를 기록한다.

  3. $\mathbf{M} = \mathbf{L}\mathbf{U}$ (불완전 분해)라는 사후조건자를 얻는다.

이 $\mathbf{M}^{-1}$는 완전 분해 대비 정확도가 떨어져, $\mathbf{A}\mathbf{x}=\mathbf{b}$에서 반복 횟수를 그다지 많이 줄여주지 못할 수 있지만, 상당히 빠르고 메모리 점유가 작다. ILU(k)는 더 많은 fill-in을 허용해 정확도를 높이는 대신, 분해 생성 비용과 메모리 사용량이 증가한다. 대규모 문제에선 ILU가 사실상 표준으로 쓰이며, 다중격자와 결합하거나, 영역분할법의 로컬 사후조건자로 이용되기도 한다.

Python의 SciPy나 Trilinos, PETSc에서는 ILU(0)나 ILU(k) 루틴이 내장되어 있어, 반복해 솔버를 호출할 때 옵션으로 지정할 수 있다. 예컨대 PETSc에서는

와 같은 커맨드라인 옵션으로 ILU(2)를 사후조건자로 하는 GMRES를 간단히 설정할 수 있다. 실제 문제에서 PCG(Preconditioned CG)를 사용할지 GMRES를 사용할지, ILU(k) 레벨이나 드롭 허용치를 어떻게 정할지는, 실험을 통해 적정선을 찾는 식으로 결정한다.

Octave를 이용한 간단한 ILU 시연

다음 코드는 Octave에서 희소행렬을 만든 뒤, ILU 사후조건자와 함께 CG를 실행해보는 예시 스크립트다. Octave는 MATLAB과 문법이 거의 동일하므로, MATLAB 환경에서도 유사하게 동작한다. (실제 ILU는 Octave 코어에 내장되지 않았을 수 있으나, 패키지를 이용하거나 MATLAB에서 이용 가능하다.)

이 코드는 ILU(0) 대신 ilutpdroptol=1e-2를 주는 방식으로 불완전 분해를 구성한다. 실제론 ILU(0)나 ILU(k)를 직접 지정할 수도 있으며, 파라미터를 조정하면 fill-in 수준과 정확도를 제어할 수 있다. 이 사후조건자 적용(함수 M1)은 삼각행렬 $L$, $U$에 대해 전진·후진 대입을 한 번씩 수행한다. 반복해 알고리즘이 한 스텝 돌 때마다 사후조건자를 한 번씩 호출하므로, 너무 복잡한 분해를 써서 사후조건자가 매우 비싸면 총 연산량이 줄어들지 않을 수 있으니, 문제 규모와 구조에 맞춰 합리적으로 선택해야 한다.

고급 주제: 불균등 계수행렬과 블록 사후조건자

물리 해석 문제나 복합 물질 문제에서, 어떤 부분은 강성(K)이 매우 크고, 다른 부분은 작은 값을 갖는 식으로 계수가 불균등하게 분포하면, ILU 분해만으로 부족해질 수 있다. 이를 해결하기 위해 특정 큰 블록을 통째로 직접 분해하고, 그 외 영역은 ILU로 처리하는 “블록 사후조건자”가 연구된다. 예를 들어 대각블록(고강성 부분)을 $LU$로 정확히 풀고, 나머지 오프 블록에는 ILU(k)나 스무딩만 적용한다. 이를 “Schur 보완(Schur complement) 기반 사후조건자”라고 부르기도 한다.

다중격자와 영역분할에서도 유사한 사고방식이 쓰인다. 경계층(boundary layer)처럼 변화가 급격한 부분을 세밀한 영역으로 분할하여, 그 구역만큼은 직접해나 고정밀 사후조건자를 쓰고, 나머지는 다중격자나 ILU로 처리하는 식이다. 이런 기법은 매우 특정한 문제(예: 유체 유동의 박리영역, 고온영역, 접촉부위 등)에서 큰 성능 향상을 가져올 수 있다. 다만 사용자가 문제 물리를 잘 이해하고, 사후조건자 구조를 세심하게 설계해야 한다는 부담이 따른다.

대규모 병렬 환경에서의 ILU 한계와 대안

ILU 분해는 가우스 소거 방식을 기반으로 하기 때문에, 프로세스 간 순차적 의존도가 발생할 수 있다. 예컨대 ILU(0)를 행별로 생성할 때, 앞선 행에서 피벗팅된 결과가 다음 행으로 전파되어야 하는 형식이어서 병렬 효율이 낮아지는 문제가 있다. 그래서 멀티코어나 다중노드 환경에서 ILU를 쓰면 스레드가 대기하는 구간이 생길 수 있다.

이를 개선하기 위해, 블록 대각분할(block diagonal partition) 후 각 블록에서 ILU를 독립적으로 수행하고, 블록 간 경계는 다르게 처리하는 병렬 ILU가 제안된다. 또는 Additive Schwarz 방법(영역분할의 일종)을 사후조건으로 사용하는데, 각 프로세스가 자기 영역에 해당하는 부분행렬을 독립적으로 ILU 분해하고, 전체 반복에서 이를 합치는 병렬 스키마다. 다중격자도 자연스럽게 병렬화 가능하므로, 실제 초대형 문제에선 AMG나 FETI/BDDC 같은 영역분할 기법이 ILU보다 병렬 효율이 좋은 경우가 흔하다.

결국 병렬성이 강한 사후조건자(AMG, DD) vs. 구현이 간단하지만 병렬성이 약한 사후조건자(ILU)라는 구도가 종종 나타난다. 어떤 것이 전반적으로 빠른지 여부는 문제 구조, 코어(노드) 수, fill-in 정도, 네트워크 대역폭 등 여러 요소를 종합해야 판명된다.

사례 확장: 코드와 실제 측정

추가 실험으로, 위 Octave 예시에선 $n=50$ 대신 $n=200$ 이상으로 키워서 $n^2=40000$ 크기 정도의 2D 포아송 문제를 풀면, ILU-기반 PCG가 없는 경우(CG만)와 있는 경우(PCG)를 비교했을 때, 반복 횟수가 얼마나 차이 나는지를 확인할 수 있다. ILU 파라미터를 조정하면 반복 횟수가 줄어드나, ILU 분해 시간이 늘어나는 점도 동시에 관찰할 수 있다.

병렬 코드를 적용하려면, MPI나 OpenMP 기반의 구현이 필요하다. MATLAB은 병렬 연산 툴박스를, Octave는 병렬 패키지를 활용하거나, 아예 C++/Fortran으로 MEX 함수를 작성해 PETSc/Trilinos를 호출하는 형태로 병렬 선형연산을 붙일 수 있다. 작은 문제 규모에서는 큰 효과가 보이지 않을 수 있으나, 실제 대규모 해석 문제를 다룰 때는 사후조건자·병렬화·도메인분할 같은 고급 기법의 영향이 극적이다.

이처럼 사후조건자 하나만 파고들어도 수많은 기법(ILU 변형, AMG, 영역분할, 블록 분해 등)과 구현 이슈(병렬 스케줄링, fill-in 제어, 유연한 적용 등)가 얽혀 있다는 사실을 알 수 있다. 결국 수치해석적 시각에서 “반복해 알고리즘”의 복잡도는 $O(Kn^2)$ 형태라 한들, 그 $K$ 안에는 사후조건 적용 비용까지 포함되어야 하므로, 사후조건자가 $O(n^2)$ 혹은 $O(n\log n)$, 때론 $O(n)$에 가능하냐에 따라 전체 성능이 완전히 바뀐다.

이후 남은 분량에서는, 다중격자(AMG)와 영역분할(DD) 기법을 좀 더 깊이 살펴보고, 실제로 중규모 예제에서 실행한 뒤 성능 비교 결과를 제시할 것이다. 또한 C++/Fortran 기반의 PETSc, Trilinos, hypre 등을 어떻게 연결해 쓰는지, 사후조건자와 Krylov 솔버를 어떤 파라미터로 튜닝하는지가 전체 시뮬레이션 시간을 얼마나 좌우하는지 살펴볼 계획이다.

다중격자(AMG)와 영역분할(DD) 심화

앞서 사후조건자 기법의 핵심 예시로 다중격자(Multigrid)와 영역분할(Domain Decomposition, DD)을 간략히 짚어 보았지만, 실제로는 이 두 방법이 단순한 사후조건자 이상의 의미를 가진다. 다중격자는 별도의 전체 알고리즘(연산자 군)으로도 활용될 수 있고, 영역분할 역시 로컬 문제를 병렬로 해결하는 대규모 프레임워크이면서, 동시에 Krylov 알고리즘의 사후조건자로 동작할 수 있다. 특히 대규모 병렬 계산 환경에서 이 둘은 매우 중요한 역할을 담당한다.

다중격자의 종류와 적용 범위

기하학적 다중격자(Geometric Multigrid, GMG)는 물리공간에서 격자가 계층적으로 구성되는 상황을 가정한다. 예를 들어 2배씩 격자 해상도를 줄이면서 레벨을 내려가면, 가장 성긴 격자에서 저주파 오차를 잡고, 상위 격자에서 고주파 오차를 잡는 방식으로 스무딩을 반복한다. 정형 격자(Structured Grid)에선 구현이 비교적 간단하고, 이론적으로도 $O(n)$ 혹은 $O(n \log n)$ 복잡도 성능을 기대할 수 있다.

하지만 실제 응용에서는 복잡한 형상의 도메인이나 비정형 격자(Unstructured Mesh)를 쓰는 경우가 많다. 이때 기하학적 다중격자 구현이 어려울 수 있으므로, 자동으로 코스 격자를 구성하는 방법(Algebraic Multigrid, AMG)이 널리 사용된다. AMG는 $\mathbf{A}$의 행렬 계수만 보고, 비영 요소의 연결성을 분석한 뒤, “상위 격자(또는 F-노드)”와 “하위 격자(또는 C-노드)”를 구분하여, 연산자 제약·보간(restriction / prolongation)을 자동 생성한다. 물리공간에서의 기하학 정보가 부족해도 적용할 수 있다는 장점이 있으나, 그만큼 초기 설정과 매개변수 조정이 까다롭다.

AMG를 사후조건자로 쓰기 위해선, Krylov 반복 루프 내에서 한 번의 사후조건자 호출이 “AMG V-cycle” 혹은 “W-cycle” 한 번에 해당되도록 설정한다. 이를 통해 “조금 무거운 사후조건자”를 쓰면서 반복 횟수를 크게 줄이는 효과를 노린다. 실제로 대규모 포아송 방정식, 탄성해석, 전자기 푸아송 문제 등에서 AMG가 매우 성공적으로 쓰이고 있다.

AMG를 설계할 때, “강 연결(strong connection)” 기준, “coarsening 기법(CF 분할)”, “스무딩 연산자 선택(가우스-자이델, SOR, 레드-블랙 등)” 등을 세심히 조정해야 한다. 예를 들어 PETSc나 hypre, Trilinos의 ML/ MueLu 패키지 등은 각각 고유의 AMG 옵션을 제공하며, 문제에 따라 수십 가지 파라미터를 튜닝하여야 최적 효율이 나온다.

영역분할(Domain Decomposition) 기법

영역분할은 물리 공간을 여러 서브도메인(subdomain)으로 나눈 뒤, 각 도메인을 독립적으로 해석하고, 경계(interface)에서만 정보를 교환함으로써 전체 문제의 해에 수렴하는 방식이다. 대표적으로 FETI(Finite Element Tearing and Interconnecting), BDDC(Balancing Domain Decomposition by Constraints), Schur Complement DD 등이 널리 알려져 있다.

영역분할 알고리즘은 고도로 병렬화하기 쉽다는 장점이 있다. 각 서브도메인에 대한 부분문제를 독립적으로 풀 때, 직접해(예: LU)나 ILU, 혹은 다중격자 같은 방법을 혼합 적용할 수 있다. 이후 서브도메인 경계에서 “글로벌 접착 조건”을 만족하기 위한 스텝을 반복하므로, 통신(경계 정보 교환)은 서브도메인 간 경계에 한정되어, 전체 스케일이 커도 병렬 스케일링이 우수하다.

DD 기법을 Krylov 사후조건자로 쓰는 경우, 예를 들어 “Additive Schwarz 방법”이나 “BDDC”를 사후조건자로 설정할 수 있다. Additive Schwarz는 “각 서브도메인이 자기 부분행렬을 분해해 놓고, 전체 해 구간에서 이것을 합산해서 한 번의 사후조건자 적용을 수행”하는 구조를 갖는다. BDDC도 유사한 개념이지만 경계 제약 조건을 좀 더 정교하게 두어, 중복 구간(경계)에서 해석 정밀도를 높이고 반복 횟수를 줄이려 시도한다.

다중격자 vs. 영역분할

다중격자(MG)와 영역분할(DD)은 상호배타적인 기법이 아니라, 종종 결합하거나 섞어서 쓴다. 예를 들어 영역분할을 한 뒤, 각 서브도메인 내부 문제를 다중격자로 푸는 방식(DD+MG), 혹은 다중격자에서 코스 격자 레벨에서 또 영역분할을 쓰는 방식(MG+DD) 등이 가능하다. 실제 초대형 시뮬레이션에서는 이런 하이브리드 구성이 자주 등장한다.

양자택일로 비교해보면, 다중격자는 “문제 전역의 고주파/저주파 성분을 분리”해 효율적 해결을 지향하고, 영역분할은 “물리 공간을 여러 조각으로 나눠 병렬·지역화”에 초점을 둔다. 실제로 대규모 HPC 환경에서는 영역분할이 병렬 스케일링 측면에서 대단히 유리하며, 다중격자도 잘 병렬화할 수 있지만, 도메인 간 통신이 많아질 수 있다. 문제 형태(단순 포아송인지, 비선형 복합 PDE인지)에 따라 어느 쪽이 더 효과적인지는 달라진다.

실제 라이브러리 예시 (PETSc)

PETSc(Portable, Extensible Toolkit for Scientific Computation)은 고성능 수치해석 라이브러리로, 다양한 Krylov 솔버와 사후조건자를 제공한다. 사용자는 간단히 커맨드라인 옵션을 통해, 예:

이런 식으로 “MG 사후조건을 갖춘 CG 알고리즘”을 실행할 수 있다. 혹은

와 같이 “영역분할(ASM: Additive Schwarz) + 서브도메인 ILU” 조합을 사용하기도 한다. PETSc는 여러 노드(MPI) 환경에서 자동으로 영역분할 도구(ParMETIS, PT-Scotch 등)를 연동해, 서브도메인을 분할하고 각 부분행렬을 관리한다. 결과적으로 사용자는 최소한의 코드 변경으로 다양한 사후조건자·영역분할·다중격자를 시험해 볼 수 있다.

Trilinos, hypre, deal.II, FEniCS 등 다른 라이브러리도 유사한 개념 구조를 제공한다. 예를 들어 hypre는 BoomerAMG라는 AMG 루틴이 유명하고, Trilinos는 MueLu, ML 패키지를 통해 다양한 AMG 스키마를 제공한다. 영역분할로는 FETI, BDDC를 구현한 Ifpack, FROSch 등이 존재한다. 다중격자와 DD 기법을 혼합하는 Hybrid~Solver도 계속 개발되고 있다.

코드 예시: PETSc 예제 스켈레톤 (C++)

간단한 C++ 스켈레톤은 대략 다음과 같은 흐름을 갖는다. 실제 코드는 훨씬 복잡할 수 있으나, 개념만 요약해보면 다음 형태다.

코드를 컴파일 후 실행 시, 다음과 같은 명령줄 옵션을 줄 수 있다.

이 명령줄은 “4개 프로세스로 병렬 실행, CG + 다중격자 사후조건, 3 레벨, 각 레벨에선 Richardson 스무딩” 설정을 나타낸다. 영역분할을 쓰려면

처럼, “GMRES + Additive Schwarz(서브도메인 ILU)”로 바꿀 수 있다. PETSc 안에는 수많은 옵션이 있어서, 세부 파라미터까지 제어 가능하다. 유저는 문제 크기와 구조에 따라 여러 조합을 실험해 보고, 반복 횟수와 실행 시간을 측정해 최적 해법을 찾는다.

수치 실험과 파라미터 튜닝

영역분할과 다중격자 기법을 적용할 때, 수많은 파라미터(예: 영역분할 기법, AMG 스무딩 방법, 다중 레벨 수, 드롭 톨러런스, 오버랩 크기 등)를 조정해야 한다. 어떤 설정이 최선인지는 실제로 문제를 풀면서 반복 횟수, 총 계산 시간, 메모리 사용량, 그리고 수치 안정성을 관찰해야 알 수 있다. HPC 센터나 산업 현장에서는 수백~수천 코어를 동원해 이 파라미터들을 실험적으로 스캔(Scan)하면서 최적의 조합을 찾기도 한다.

또한 문제 특성(조건수, 희소 패턴, 비선형 coupling)에 따라 상이한 결과가 나와, 사용자 경험이나 축적된 노하우가 중요하다. 예컨대 특정 PDE 문제에서는 AMG가 매우 탁월한 성능을 보이는 반면, 어떤 복잡한 파동 문제에서는 AMG가 고전하고 영역분할이 더 낫기도 하다. 접촉 문제나 경계가 이동하는 문제, 멀티피직스 연성 문제 등에서는 각 서브문제가 서로 다른 솔버·사후조건자를 쓰면서 공동으로 수렴해야 해서, 훨씬 복합적인 시나리오가 전개된다.

결국 대규모 선형시스템 해석은 “알고리즘 + 구현 + 병렬화 + 문제 특성 + 사용자 튜닝”이 맞물린 종합 예술에 가깝다. 본 장에서 개념적 토대를 충분히 다지고, 실제 코드·라이브러리 예시와 더불어 성능 측정 결과를 해석하는 과정을 거치면, 자신에게 맞는 최적의 솔루션을 발견하는 데 큰 도움이 될 것이다.

실제 사례 종합과 추가 연구 방향

이제까지 다양한 기법(직접해, 반복해, 사후조건자, 다중격자, 영역분할 등)이 어떻게 적용되고 어떤 연산 복잡도를 갖는지, 그리고 실제 대규모 문제에서 성능이 어떻게 나타나는지 개략적으로 살펴보았다. 마지막으로, 본 장에서 다룬 내용 전반을 간단히 종합하고, 추가 연구와 실무에서 고려해야 할 몇 가지 방향을 제시한다.

사례 종합

앞서 제시한 사례들을 간략히 되짚어 보면, 공통적으로 “직접해 vs. 반복해” 구도가 핵심을 이룬다. 연산 복잡도만 놓고 보면, 직접해는 $O(n^3)$이고, 반복해는 $O(K n^2)$라고 쉽게 요약되지만, 실제로는 아래와 같은 변수가 뒤섞인다.

직접해의 경우:

  • $n$이 중간 수준(수천~수만)이라면 고성능 BLAS-3 병렬연산으로 빠르게 해결 가능.

  • 단 한 번의 풀이가 아니라, 여러 번 우변이 바뀌면서 해를 구해야 하는 경우도 유리. 분해 이후에는 $O(n^2)$로 해를 빠르게 구할 수 있기 때문.

  • $n$이 수십만~수백만 이상으로 넘어가면, $O(n^3)$이 현실적으로 불가능해지거나, 대규모 메모리 fill-in 문제로 인해 분해가 끝나기조차 힘듦.

반복해의 경우:

  • 한 번 반복에서 행렬-벡터 곱셈이 $O(n^2)$ (밀집행렬) 또는 $O(n)$ (희소행렬) 정도로 가능.

  • 반복 횟수 $K$가 커지면 결국 $O(Kn^2)$ 혹은 $O(Kn)$(희소인 경우)로 계산 비용이 달라져, $K$가 조건수나 사후조건자 품질에 의해 크게 좌우됨.

  • 조건수가 매우 좋지 않으면 $K$가 $O(n)$ 가까이 증가할 위험이 있으며, 이럴 때는 사실상 $O(n^3)$에 이르는 연산이나 더 큰 비용을 치를 수도 있음.

  • 대규모 희소행렬 문제에서 사후조건자(ILU, AMG, DD 등)를 잘 쓰면 $K$를 거의 상수 수준으로 묶을 수도 있음. 이 경우 $O(n)$ 혹은 $O(n\log n)$에 준하는 놀라운 성능이 가능.

추가로 병렬성, 메모리 구조, 문제 특성(물리적 PDE냐, 통계 회귀냐, 커널행렬이냐 등), 그리고 정확도/안정성 요구사항에 따라 어느 알고리즘이 최적이 될지는 크게 달라진다. 실제 산업·연구 환경에서는, 다음과 같은 판단 흐름이 일반적이다.

주어진 문제의 스케일($n$), 희소성 여부, 대칭성·양의정부호 여부, 그리고 해석 시간 한계·메모리 한계 등을 먼저 파악하고, 최적·근사 기법 중에서 하나를 고른다. 문제 규모가 “중간”이라면(예: $n$이 수천~수만), 직접해나 반복해 둘 다 시도해보며, 특정 하드웨어에서 최적화를 수행한다. 문제 규모가 “초대형”이라면, 사실상 반복해 + 사후조건 + 병렬화가 필수적이어서, 더 이상 직접해가 후보가 되기 어렵다. 이때 사후조건자를 어떻게 만들고, 어떤 병렬 분산 구조를 쓸지가 핵심 설계 포인트가 된다.

후속 연구·실무 방향

  1. 병렬·분산 튜닝 자동화 HPC 라이브러리(PETSc, Trilinos, hypre 등)는 각종 사후조건자, Krylov 솔버, 도메인 분할, AMG 파라미터를 매우 세밀하게 제어할 수 있도록 옵션을 제공한다. 이 옵션들이 많아질수록, 사용자의 작업량도 커지는데, 이를 자동화(automatic tuning)하는 연구가 활발하다. 예컨대 소프트웨어가 문제 구조를 자동 분석해 “현재 행렬에서 ILU(2) + GMRES가 좋겠다”, “다중격자 레벨 수는 4가 적당한다” 같은 결정을 내려주도록 하는 것이다. 이런 전문가 시스템은 아직 완성도가 높지 않지만, 미래에는 “자동 솔버 선택 프레임워크”가 한층 진화해, 일반 사용자도 대규모 선형시스템을 효과적으로 풀게 될 가능성이 있다.

  2. 혼합정밀도·비정밀 아키텍처 활용 GPU나 특수 AI 칩은 배정밀도(double precision)보다 단정밀도(single), 반정밀도(half)의 연산 효율이 매우 높다. 따라서 반복해 과정 초기나, 오차가 여전히 큰 단계에서는 저정밀로 빠르게 연산하고, 후반부에만 고정밀로 전환하는 혼합정밀도 기법이 고려된다. 또, 새로운 아키텍처(정수 기반 가속기, 유사 아날로그 컴퓨팅, 양자컴퓨팅 등)에서도, 선형시스템을 근사로 해결하는 시도가 이어지고 있다. 이들은 아직 초기 단계지만, 오래된 $O(n^3)$ 혹은 $O(Kn^2)$ 벽을 넘을 새로운 길을 열어줄 수도 있다.

  3. 멀티피직스·적응형 문제에서의 동적 행렬 실제 문제는 대규모 비선형, 시변(時變), 연성(連成) 형태를 띠는 경우가 많다. 하나의 선형계수행렬을 계속 사용하는 것이 아니라, 매 시간 스텝 혹은 매 뉴턴 반복마다 행렬이 갱신되어 버린다. 이렇게 빈번히 갱신되는 환경에서는 사후조건자나 분해를 매번 다시 만들기 어렵다. 이에 대한 해법으로 ‘partial factorization update’나 ‘AMG 계층 부분 재구축’, ‘증분(incremental) ILU’ 등의 아이디어가 시도되고 있다. 이 분야는 연구 여지가 매우 크며, 성공적 접근을 통해 동적 문제에서도 효율적인 선형시스템 풀이가 가능해진다.

  4. 이산화(direct discretization) vs. 근사화(Approx. Methods) 조합 PDE 기반 문제에서, 격자를 더 정교하게 만들수록 선형시스템 규모가 폭발적으로 늘어난다. 반면, 물리적으로 매우 세밀한 결과가 필요하지 않은 구간에서는 격자를 줄이거나, 저랭크 근사 혹은 모델 축소(reduced-order modeling)를 사용해 전체 문제 차원을 낮출 수 있다. 머신러닝에서도, 초대형 커널행렬을 직접 쓰는 대신, 랜덤 스케치나 Nystrom 근사로 행렬을 압축하여 $O(n^2)$ 이상의 연산을 회피한다. 이런 접근이 확대되면, 전통적 $O(n^3)$·$O(Kn^2)$ 알고리즘을 완전히 대체할 수도 있다.

  5. 하이브리드 HPC 환경 최근 슈퍼컴퓨터는 CPU+GPU만이 아니라, 여러 종류의 가속기나 프로세서를 혼합한 이기종(heterogeneous) 아키텍처로 구성되는 추세다. 이를 효율적으로 사용하려면, 특정 연산은 GPU가 맡고, 다른 부분은 CPU가 맡고, 데이터 이동은 최소화하는 식의 정교한 스케줄링이 필수다. 선형시스템 풀이는 워낙 핵심 연산들이 많으므로, 여러 장의 GPU, 다양한 프로세서 코어, 그리고 빠른 노드 간 네트워크를 어떻게 최적으로 묶어 쓰느냐가 고성능 달성의 관건이 된다.

간단 실험 예시 확장

마지막으로, 실제로 학습자가 이번 장에서 소개된 알고리즘을 시험해 보려면 다음과 같은 확장 방향을 권장할 수 있다.

Python 또는 Octave/Matlab 환경에서, 작은·중간 규모(수천~수만 차원)의 문제를 대상으로

  • 가우스 소거법(직접해) vs. CG, GMRES 등 반복해를 각각 실행

  • 사후조건자(ILU, 간단 다중격자(Geometric) 등) 적용 전후로 반복 횟수와 시간 측정

  • 문제 구조(대각우세, 심한 비대칭, 불균등 계수 등)를 바꿔가며 결과 비교

이 과정을 병렬 라이브러리(예: Python의 multiprocessing 혹은 C++/Fortran에서 MPI 등)와 연계하면, 코어 수 늘림에 따른 계산 시간 변화를 확인할 수 있다. 또한 GPU가 이용 가능한 환경(CuPy, PyTorch, CUDA C/C++)이라면, 일부 연산(행렬 분해, 행렬-벡터 곱)을 GPU에서 수행했을 때의 성능 이점을 체감해볼 수 있다.

좀 더 고급 수준에 도달하면, PETSc, Trilinos, hypre 같은 전문 HPC 라이브러리를 설치해, 명령줄 옵션만으로도 수십 종의 솔버·사후조건자를 조합할 수 있다. 이때 스케일 테스트(예: $n$을 10배씩 키워가며 실행)와 병렬 스케일 테스트(프로세스 수를 1배, 2배, 4배 등으로 늘리며 실행)를 해 보면, 이론적 복잡도와 실제 실행 시간의 괴리가 어디에서 발생하는지 명확히 알 수 있다.

중간 정리

현대 수치해석에서 선형시스템 해결은 수많은 응용 분야(물리 시뮬레이션, 머신러닝, 최적화, 역문제 등)에서 변함없이 중요한 과제로 남아 있다. 본 장에서 다룬 연산 복잡도 이론, 직접해·반복해 알고리즘, 사후조건자·병렬화 개념, 그리고 다양한 응용 사례들은 향후 독자가 실제 프로젝트나 연구를 진행할 때 밑바탕이 될 것이다.

예제 코드를 통한 학습 방향

본 장에서 다루었던 내용을 요약한 뒤, 실제로 코드를 실행해 보며 이해를 심화할 수 있는 몇 가지 예시를 제안한다. 간단한 스크립트 예제부터 본격적인 HPC 라이브러리 사용까지, 점진적으로 접근하는 순서를 권장한다.

1. Python/Octave/Matlab의 기본 함수 활용

  • 직접해 vs. 반복해

    • 임의 크기의 밀집행렬(또는 간단한 희소행렬)을 만들고, np.linalg.solve(Python)나 linsolve, spsolve(Octave/Matlab 등)를 써서 직접해를 구한다.

    • 같은 문제에 대해 CG, GMRES 등을 반복해로 풀어보고, 연산 시간과 반복 횟수를 측정한다(사후조건자 없이).

    • 문제 크기 $n$을 서서히 늘려가면서 실행 시간을 로그-로그 그래프로 그려, 이론적 $O(n^3)$, $O(Kn^2)$와 대조해본다.

  • 사후조건자 도입

    • 희소행렬 포아송 예제(2D/3D)를 구성한 뒤, ILU(또는 Incomplete Cholesky) 사후조건자를 적용한 PCG(Preconditioned CG)나 GMRES를 실행한다.

    • 드롭 톨러런스, ILU 레벨, 스무딩 횟수 등을 달리해서 반복 횟수와 실제 계산 시간을 비교한다.

    • (Octave/Matlab) ilu, ichol 함수, (Python SciPy) spilu 등을 활용할 수 있다.

2. 간단 병렬화 테스트

  • OpenMP 병렬화

    • CSR 형식 행렬과 벡터 곱셈 루틴을 C++로 작성하고, #pragma omp parallel for를 붙여서 멀티코어 속도 향상을 체험한다.

    • 문제 크기 $n$을 늘리거나 스레드 수를 바꿔가며, 성능 스케일링을 관찰한다.

  • MPI 병렬화

    • 여러 노드(혹은 여러 프로세스)에서 문제 일부를 나누어 저장하고, 행렬-벡터 곱셈 과정에 필요한 경계 데이터를 MPI로 교환한다.

    • 간단한 CG 알고리즘을 전체 MPI 환경에서 수행해 보고, 코어 수나 노드 수를 늘릴 때 시간(또는 반복 횟수)이 어떻게 변하는지 측정한다.

이 과정에서 “통신 vs. 계산” 비중이 어떻게 달라지는지 감을 잡을 수 있고, 실제로 높은 병렬 효율을 내려면 어떤 최적화가 필요한지 느낄 수 있다.

3. 전문 HPC 라이브러리(PETSc, Trilinos, hypre 등) 사용

  • 기본 예제 컴파일

    • PETSc 예제를 내려받아 컴파일하고, -ksp_type cg, -pc_type ilu, -pc_factor_levels 1, -pc_factor_levels 2 등의 옵션을 바꿔가며 실행한다.

    • -ksp_type gmres, -pc_type mg (AMG) 조합도 시험해 본다.

    • 문제 크기를 단계적으로 늘리며, 반복 횟수와 전체 시간, 메모리 사용량을 로그에 기록한다.

  • 영역분할(Additive Schwarz, FETI, BDDC) 기법

    • -pc_type asm, -sub_pc_type ilu처럼 영역분할 + ILU 조합을 사용해 본다.

    • 코어 수(MPI 프로세스 수)를 늘리면서 실행 시간을 비교하고, 병렬 스케일링 그래프를 그려본다.

이런 실습을 통해 “이론적 복잡도”와 “실제 계산”의 차이를 체감할 수 있으며, 어떤 옵션이 수렴 횟수를 획기적으로 줄여주는지, 혹은 어느 순간 병렬 통신이 병목으로 작용하여 스케일링이 제한되는지도 확인 가능하다.

4. 문제 특성 바꾸기

  • 조건수 변경

    • 대각우세(diagonally dominant)를 강화하거나, 서로 다른 물성(강성) 구간을 주어 행렬의 조건수를 악화시키는 식으로 실험한다.

    • 반복해 솔버의 수렴이 얼마나 민감한지 확인한다. 사후조건자가 없으면 사실상 수렴하기 어려운 상황도 마련해 볼 수 있다.

  • 비대칭·복소행렬

    • 유체역학(나비에-스토크스)나 전자기(Helmholtz) 문제에서 비대칭·복소 행렬을 만들어, CG 대신 GMRES나 BiCGSTAB 등으로 풀어 본다.

    • AMG가 비대칭 행렬에서 잘 작동할지, 영역분할이 더 좋은 결과를 낼지 직접 비교한다.

이렇게 문제 구조를 바꿀 때, 라이브러리 옵션이나 사후조건자 종류도 맞춤형으로 세팅해야 하며, 그렇지 않으면 수렴이 매우 느려지거나 아예 실패하기도 한다.

학습 마무리

위와 같은 실습 과정을 거치면, 본 장에서 배운 “직접해 vs. 반복해”, “사후조건자와 병렬화”의 중요성이 구체적으로 체감된다. 실제 응용 현장에선 훨씬 더 큰 문제와 다양한 물리·기계학습 모델이 얽히지만, 기본 원리는 동일하다. 문제 규모가 커질수록, 올바른 알고리즘 선택과 튜닝이 수십~수백 배의 성능 차이를 만들 수 있다는 사실을 잊지 말아야 한다.

Last updated