대규모 행렬에 대한 수치적 솔루션 요약
직접 해법과 대규모 문제에서의 어려움
밀집 행렬에 대한 직접 해법은 $\mathbf{A}\mathbf{x} = \mathbf{b}$ 형태의 선형 방정식을 풀 때, 일반적으로 LU 분해 혹은 Cholesky 분해와 같이 행렬 분해 기법을 적용하는 과정을 거친다. 그러나 문제 차원이 매우 커지면 메모리 사용량과 연산 횟수가 급격히 증가하므로, 대규모 행렬을 다룰 때에는 행렬의 구조적 특성(스파스 행렬, 대각우선 행렬 등)을 최대한 활용할 필요가 있다. 직접 해법은 이론상 유한 단계 안에 해를 구하지만, 스파스 구조가 충분히 희박하지 않다면 분해 과정에서 발생하는 fill-in이 심해져서 실제 연산 복잡도가 매우 커질 수 있다.
직접 해법에서는 $\mathbf{A}$를 적절히 분해하는 과정을 거치는데, 예를 들어 비특이이고 일반적인 정방 행렬에 대해 피벗팅을 고려한 LU 분해를 수행하면
형태로 표현하고, 순차적으로 전방 대치와 후방 대치를 통해 해 $\mathbf{x}$를 구한다. 행렬 $\mathbf{A}$가 대칭 양의정부 행렬(symmetric positive definite)이면 Cholesky 분해를 사용하여
형태로 나타내어 계산량을 절감할 수도 있다. 하지만 이때에도 대규모 문제에 대해서는 fill-in을 줄이기 위한 리오더링 기법이나 멀티프로세서 환경에서의 병렬 연산 기법이 필수적으로 고려되어야 한다.
행렬의 크기가 방대해지면 메모리 계층 구조와 병렬 연산 모델을 효과적으로 설계해야 한다. 캐시 친화적(block-based) 알고리즘이나 압축 포맷(CSR, CSC 등)의 활용은 커다란 연산 비용을 획기적으로 줄이는 데 기여한다. 실제 HPC 환경에서는 행렬을 여러 노드에 분산하여 저장하고, 통신량 최소화를 목표로 스케줄링을 진행한다. 직접 해법은 안정성과 정확도가 높지만, 용량과 계산 속도 측면에서 큰 제약이 따른다.
반복 해법과 Krylov 하부공간 기법
대규모 문제에서는 주로 반복(iterative) 해법이 널리 활용된다. 반복 해법은 초기값 $\mathbf{x}_0$에서 시작하여, 잔차(residual) 벡터 $\mathbf{r}_k = \mathbf{b} - \mathbf{A}\mathbf{x}_k$를 개선하는 방식으로 해에 수렴시키는 절차를 거친다. 이때 Jacobi나 Gauss-Seidel 등 고전적 방식은 간단하지만 큰 문제에서 빠르게 수렴하지 않거나 수렴 보장이 느슨할 수 있다. 이를 보완하기 위해 Krylov 하부공간(Krylov subspace) 기법을 사용한 GMRES, CG, BiCGSTAB 등이 중점적으로 활용된다.
Krylov 하부공간 기법에서 중요한 아이디어는 잔차 벡터 $\mathbf{r}_k$와 $\mathbf{A}\mathbf{r}_k$ 등을 순차적으로 생성하여, 미리 정의된 차원을 갖는 하부공간에서 최적 근사해를 찾는 방식으로 수렴성을 높인다는 점이다. 예를 들어 CG 방법의 경우 대상 행렬이 대칭 양의정부라는 가정 하에서, 잔차 벡터와 검색 방향(search direction)을 반복적으로 조정하여 빠른 수렴을 유도한다. GMRES는 일반적 행렬에 대해서도 적용 가능하지만 잔차의 직교화 과정에 의한 메모리 사용량 증가 문제가 있으므로, 대규모 문제에서는 Restarted GMRES(m) 방식이나 다른 변형이 자주 쓰인다.
Krylov 기법에서는 적절한 사전조건(preconditioner)을 선정하는 과정이 필수적이다. 사전조건화는 원래 문제 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 $\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}$ 형태로 바꾸어, 스펙트럼 분포를 개선하고 반복법의 수렴을 가속화한다. 사전조건 행렬 $\mathbf{M}$으로는 간단한 대각 행렬부터 ILU(불완전 LU) 분해처럼 보다 정교한 행렬 분해 기반 기법까지 다양한 방법이 시도된다. ILU(k)는 희소구조 유지 정도를 조절하기 위해 fill-in 수준인 k를 설정하여, 적당한 균형점을 찾는다.
고차원 병렬 환경과 분산 메모리
반복 해법을 적용할 때에도 결국 매 반복 단계에서 $\mathbf{A}\mathbf{v}$ 형태의 행렬-벡터 곱셈이 자주 호출된다. 분산 메모리 병렬 환경에서는 행렬과 벡터가 노드별로 나누어 저장되므로, $\mathbf{A}\mathbf{v}$ 연산 시 통신 비용이 성능 병목을 일으키게 된다. 적절한 도메인 디컴포지션(domain decomposition)과 데이터 재배치, 비정형 메시에 대한 그래프 파티셔닝 기법 등으로 통신량을 최소화해야 한다. 대규모 HPC 환경에서는 프로세서 개수에 따른 강한 스케일링(strong scaling)과 약한 스케일링(weak scaling)을 함께 평가하여 알고리즘 효율성을 검증한다. 실제 구현에서는 OpenMP, MPI, CUDA 등의 다양한 병렬 프로그래밍 모델이 고려된다.
Octave 환경에서 Krylov 기반의 반복 해법 실험을 해보고 싶다면 아래와 같은 예시로 시작할 수 있다.
위 코드는 무작위 스파스 행렬을 생성한 뒤, 초기값을 $\mathbf{x}_0 = \mathbf{0}$로 두고 GMRES로 풀어 보는 간단한 스켈레톤이다. 실질적 문제에 적용하려면 사전조건 행렬을 추가로 지정해야 하며, 보다 큰 규모의 예제에서는 분산 메모리 환경에서의 병렬 연산(예: MPI)을 고려해야 한다.
도메인 분할 기법
물리적 해석 영역을 여러 하위 도메인으로 분할하고, 각 서브도메인에 할당된 소규모 문제를 개별적으로(또는 병렬적으로) 풀어 전체 해를 조합하는 접근을 도메인 분할(domain decomposition) 기법이라 한다. 크게 Dirichlet-Neumann 방식, Robin-Robin 방식, Schur 보충(Schur complement) 방식 등이 제안되어 왔으며, 공학 문제에서 주로 보이는 불균일한 물리적 특성이나 복잡한 경계 조건에 대해 잘 대응할 수 있다는 장점이 있다.
도메인 분할법의 핵심은 각 서브도메인 내에서는 국소적인 해를 빠르게 구하고, 서브도메인 간 경계 조건을 만족시키기 위한 글로벌 적합성(global compatibility) 조건을 반복적으로 조정해 가면서 전체 문제의 해를 얻는 것이다. 예를 들어 2차원 영역을 다룰 때, 영역을 여러 블록으로 나누고 각 블록 내의 미분방정식을 선형화한 뒤 스파스 선형계 $\mathbf{A}_i \mathbf{x}_i = \mathbf{b}_i$ 형태로 구성한다. 여기서 $\mathbf{A}_i$는 각 서브도메인에 해당하는 내부 자유도와 경계 자유도를 함께 포함한다. 전역 문제와 서브도메인 문제 사이를 연결하기 위해 Schur 보충 기법을 적용하면
형태에서 경계 블록만을 다루는 축소 문제(보완 문제)를 유도하고, 이를 반복적으로 업데이트하여 전체 해를 구하게 된다.
도메인 분할법은 병렬성이 높고 국소 계산이 집중되는 방식이어서, 대규모 HPC 환경에서 스케일링을 달성하기 수월한 편이다. 하지만 경계 혹은 서브도메인 간 접합 지점에서 발생하는 수렴 지연이나 인접 도메인들 간의 조화를 맞추기 위한 통신 비용이 문제의 스케일에 따라 커질 수 있으므로, 적절한 그래프 파티셔닝 알고리즘과 서브도메인별 사전조건의 효율적 설계를 함께 고려해야 한다.
멀티그리드 기법
멀티그리드(multigrid)는 PDE 유도 문제에서 가장 널리 활용되는 고속 전역 반복 해법 중 하나이다. 일반적으로 격자해석(Finite Element/Finite Difference 등)에서 만들어진 스파스 행렬에 대해, 미세 격자(fine grid)와 더 거친 격자(coarse grid)를 오가면서 잔차를 효율적으로 소거하는 전략을 쓴다.
멀티그리드 절차에서는 보통 V-cycle이나 W-cycle을 적용하여, 미세 격자에서 얼마간 스무딩(smoothing) 반복을 수행한 뒤 거친 격자로 제약(prolongation/restriction)을 통해 문제 크기를 축소시킨다. 거친 격자에서 계산된 근사해 정보를 다시 미세 격자로 확대(prolongation)하여 해를 개선하는 식이다. 이때 거친 격자 문제 자체도 또 다른 멀티그리드를 적용하여 더 작은 격자에서 풀어낼 수 있다. 멀티그리드는 저주파(장파장) 오류 성분을 제거하는 데 탁월하므로, 전역적인 수렴을 빠르게 달성한다.
멀티그리드 접근은 미분 방정식에 최적화된 형태로 개발되어 왔지만, 선형계 일반에도 대체로 적용 가능하다. 기하학적(gMG) 접근은 격자 구조가 명확한 상황에서 자연스럽게 여러 레벨의 그리드를 구성한다. 반면 대규모 비정형 스파스 행렬에 대해서는 기하학적 관점이 어려우므로, 강약 연결(strong/weak connection)을 기반으로 AMG(Algebraic Multigrid) 기법을 적용한다. AMG는 행렬 계수 자체만을 기반으로 상위 레벨에서 하위 레벨로 갈수록 집합을 묶어(coarsening) 더 작은 문제를 형성하며, 적절히 반복 및 보간(interpolation) 연산을 정의해 준다.
멀티그리드를 사전조건으로 도입하여 Krylov 하부공간 기법과 결합하면, 매우 뛰어난 수렴성을 기대할 수 있다. 예를 들어 AMG를 CG나 GMRES의 사전조건으로 적용하면, 고차원 문제에서도 반복 횟수가 크게 줄어드는 효과가 있다. 그러나 멀티그리드 레벨 구성이나 스무딩 기법, 보간 연산자 설계 등 구현 세부가 복잡하고, 문제 특성에 따라 적절한 파라미터 튜닝이 필요하다.
사전조건 병합 전략
도메인 분할 사전조건법(DD preconditioner)이나 멀티그리드 기반 사전조건법, ILU 등 서로 다른 사전조건을 함께 결합하여 복합적인 효과를 노리는 전략도 대두되고 있다. 예를 들어 부분영역 사전조건과 AMG 사전조건을 계층적으로 병합하면, 지역적 해부와 전역적 잔차 제거 효과를 동시에 달성할 수 있다. 구체적인 구현은 문제 크기, 스파스 패턴, 병렬 하드웨어 구성, 계산 시간 대비 정확도 요구사항 등을 모두 고려해야 한다.
아래에 Python으로 기초적인 AMG 사전조건을 활용한 CG 예시 스켈레톤을 제시한다. 실무에서는 pyAMG, hypre, PETSc 등의 라이브러리가 종종 쓰인다.
위 예제에서는 간단한 1차원 2차 미분 방정식 계수를 대각선이 2, 양 옆 대각선이 -1인 삼중대각 행렬로 구성하여 $\mathbf{A}$를 만든 뒤, pyAMG를 이용해 RUUGE-STüBEN 방식의 AMG 사전조건을 생성하고 CG로 문제를 푸는 과정을 보여 준다. 대규모 문제일수록 사전조건의 설계와 병렬 구현이 중요해지며, 실전에서는 수십~수백 개 이상의 프로세서 노드를 활용해 전체 부하를 분산하는 식으로 확장된다.
GPU 가속 및 혼합 정밀도 기법
행렬 연산은 병렬성이 높아 GPU(Graphics Processing Unit) 가속의 이점을 크게 누릴 수 있다. 대규모 문제에서 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 풀 때, GPU에 대량의 연산을 오프로드하여 CPU-GPU 간 협업을 진행하면 전반적인 시간 단축 효과가 매우 커진다. 특히 반복적 방법에서는 매 반복마다 행렬-벡터 곱셈이 반복 수행되므로, GPU의 대규모 병렬 처리 능력이 극대화된다. 그러나 GPU 메모리 대역폭은 여전히 한계가 있고, CPU와의 통신 경로 역시 병목이 될 수 있으므로, 효율적인 행렬 저장 형식(CSR, ELLPACK 등)과 최적화된 커널 설계가 필수적이다.
GPU 환경에서 혼합 정밀도(mixed-precision) 방식을 적용할 경우, 계산의 주요 단계는 단정밀도(single precision)나 반정밀도(half precision)로 빠르게 처리하고, 최종 보정 단계에서 배정밀도(double precision)로 해를 정밀하게 보강하는 전략을 취할 수 있다. 이는 연산 성능이 높은 저정밀도 코어를 적극 활용하되, 결과 해의 정확도를 떨어뜨리지 않으려는 시도로, 최근 고성능 컴퓨팅에서 각광받는 기법 중 하나다. 반복적 해법이나 다단계(multi-stage) LU 분해 등에 혼합 정밀도가 적용될 수 있는데, 예를 들어 아래와 같은 형태의 반복적 보정(Iterative Refinement)과 결합하면 실질 오차를 안정적으로 줄여나갈 수 있다.
이 과정을 여러 차례 반복하여 $\mathbf{x}_\text{refined}$를 갱신하면, 궁극적으로 배정밀도 수준의 정확도를 달성하면서도 전체 계산량은 상당 부분 저정밀도 연산을 활용하여 시간 비용을 크게 절감할 수 있다.
GPU상에서 혼합 정밀도 연산을 수행하기 위해서는 벤더별 라이브러리(CUDA, ROCm 등)나 고수준 프레임워크(CuPy, PyTorch, TensorFlow 등)의 저정밀도 BLAS/LAPACK 지원 기능을 적절히 호출해야 한다. 다만 직접 LU 분해를 GPU 상에서 혼합 정밀도로 구현할 경우 피벗팅 처리, fill-in 관리, 컬럼 병합 전략 등에 대한 심도 있는 연구가 필요하다.
고급 반복적 보정과 적응형 정밀도 전략
표준적인 반복적 보정(Iterative Refinement) 기법은 잔차 벡터를 배정밀도에서 다시 풀어 수정항 $\mathbf{\delta}$를 계산하는 구조를 가진다. 여기에 적응형 정밀도(adaptive precision) 전략을 결합하면, 문제 상황에 따라 저정밀도와 고정밀도 연산 간 전환을 동적으로 수행할 수 있다. 예를 들어 잔차가 충분히 작을 때는 단정밀도 곱셈으로도 오차 감소가 빠르지만, 잔차가 일정 임계값 이하로 떨어지면 배정밀도 곱셈으로 세부 보정을 수행하는 식이다. 이를 통해 전체 알고리즘의 병렬 효율성을 유지하면서 메모리 전송량과 연산 시간을 균형 있게 줄일 수 있다.
전통적인 부분 피벗팅(partial pivoting)이나 스케일링(scaling) 기법은 배정밀도에서 수치적으로 안정적인 해를 구하는 데 핵심적인 장치이지만, 혼합 정밀도 환경에서는 일부 피벗팅 동작이 저정밀도 연산에서 발생할 수 있으므로, 결과적으로 불안정이 커질 수도 있다. 따라서 대규모 문제에서는 사전에 행렬의 행 또는 열을 재배열(reordering)하거나, 일괄적인 pivot-free LU 분해를 시도하는 등의 사전조건 기법을 병행하기도 한다.
대규모 HPC 소프트웨어 프레임워크
현실적으로 대규모 스파스 행렬 문제를 풀기 위해서는 PETSc, Trilinos, hypre 등 고성능 병렬 연산 지원 라이브러리를 활용하는 것이 일반적이다. 이러한 프레임워크들은 분산 메모리(MPI) 환경에서의 데이터 분할과 통신, 다양한 사전조건(AMG, Block Jacobi, ILU 등) 및 반복 해법(CG, GMRES, BiCGSTAB 등)을 손쉽게 구성할 수 있도록 고수준 API를 제공한다. GPU 가속 기능과 혼합 정밀도 연산 모듈도 지속적으로 확장되고 있으며, 3차원 구조나 적응형 격자(Adaptive Mesh Refinement) 문제에 대한 실무 지원도 활발하다.
Trilinos의 AztecOO, Ifpack, MueLu 모듈, PETSc의 PCMG, PCHYPRE, PCLU 등은 도메인 분할 및 멀티그리드 기반 사전조건을 고도로 추상화하여, 사용자 입장에서 반복 해법을 다양하게 시험해 볼 수 있게 해 준다. 예컨대 PETSc를 이용해 선형계 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 풀이하려면,
와 같이 실행 인자를 통해 쉽게 알고리즘을 교체할 수 있다. 실제 코드는 C, C++ 또는 Python 바인딩을 통해 작성하고, 내부적으로는 PETSc 라이브러리가 MPI를 동원하여 병렬 연산을 수행한다. 행렬 구조를 효율적으로 관리하기 위해 AIJ(압축 행렬), BAIJ(블록 압축 행렬) 등 다양한 스파스 포맷을 지원한다.
실험 및 디버깅에서의 고려사항
대규모 문제를 다루다 보면 수렴이 지연되거나, 어떤 특정 노드에서 과도한 부하가 발생하는 로드 언밸런스(load imbalance) 문제, 혹은 수치적 불안정성으로 인한 해의 발산(divergence)이 관찰될 수 있다. 따라서 병렬 성능 분석 툴(예: Tau, Scalasca, Intel VTune 등)을 사용하여 병목이 되는 부분을 찾고, 도메인 분할 과정이나 사전조건 설정을 반복 조정하면서 최적의 성능을 모색해야 한다. 실제로 PDE 유도 문제에서 격자 세분화가 극단적으로 진행될 때는 경계 조건이나 계수 분포가 숫자적으로 문제를 더 어렵게 만들 수 있으므로, $\mathbf{A}$의 상태수(condition number)에 대한 사전 조사, 스케일링, 임계 영역에 대한 거친 격자 보강 등 다양한 기법을 시도해야 한다.
아래에는 C++에서 PETSc를 이용해 분산 메모리 환경에서 간단한 선형계 GMRES 해법을 설정하는 예시 스켈레톤을 보여 준다. 실제 대규모 문제에서는 사용자 정의로 행렬 구조를 적절히 분산할 필요가 있으며, 사전조건이나 반복자 설정은 명령줄 옵션 또는 API를 통해 유연하게 변경할 수 있다.
이런 식으로 PETSc 코드를 작성해 두면, 실제 실행 시에
등의 옵션을 가감하여 반복자와 사전조건을 실험적으로 조정할 수 있다. 대규모 HPC 시스템에서 노드 수를 늘림에 따라, 스케일링 테스트를 수행해 성능 향상도를 관찰하는 방식으로 튜닝을 진행한다.
하이브리드 접근과 로우랭크 근사
대규모 스파스 행렬을 효율적으로 풀기 위해 직접 해법과 반복 해법을 혼합한 하이브리드 방식이 모색되기도 한다. 예를 들어 큰 행렬을 여러 부분블록으로 분할한 뒤, 각 블록 내부는 스파스 직접 해법(LU, Cholesky 등)으로 풀고 블록 간 상호 작용은 반복 해법으로 처리함으로써, 국소 영역의 성능과 전역적 확장성을 동시에 추구한다. 이런 분할-반복 접근법은 유한요소법이나 볼륨요소법을 적용한 PDE 문제에서 특히 자주 언급된다.
로우랭크 근사(low-rank approximation)는 엄밀한 스파스 분해 대신, 근사 행렬 $\mathbf{A}_\text{approx}$를 사용하여 연산량과 메모리를 크게 절감하는 기법이다. 예를 들어 H-행렬(Hierarchical matrix)이나 $H^2$-행렬, FMM(Fast Multipole Method) 등에 기반을 두어, 원래 행렬을 계층적으로 블록 분할하고 각 블록을 저랭크 행렬 형태로 가까이 근사한다.
이런 구조를 갖는 블록들이 계층적으로 중첩되어 전체 행렬을 구성할 수 있다고 가정하면, 곱셈 및 분해 연산의 지배적인 비용이 기존 스파스 알고리즘보다 훨씬 감소한다. 대규모 3차원 PDE 솔버나 고차원 적분방정식에서 FMM 기반의 H-행렬 기법이 도입되는 이유가 여기에 있다.
H-행렬은 블록 계층을 구성할 수 있는 기하학적(또는 그래프) 정보가 명확해야 하며, 각 블록 내에서 행렬을 어떠한 랭크로 근사할지(트렁케이션 임계값) 결정하는 것이 관건이다. 또한 해가 반복적으로 개선되는 과정에서 행렬 근사 오차가 수렴성을 저해하지 않도록, 사전조건처럼 적절히 혼용하여 사용하는 전략이 필요하다. 이를 통해 스파스 직접 해법 대비 fill-in을 억제하고, 반복 해법 대비 수렴 가속을 동시에 기대할 수 있다.
통신 회피(communication-avoiding) 알고리즘
반복 해법에서 가장 빈번하게 발생하는 연산은 $\mathbf{A}\mathbf{v}$ 형태의 행렬-벡터 곱셈이지만, 실제 분산 메모리 환경에서는 노드 간 통신이 연산 수행보다 더 큰 병목이 될 수 있다. 따라서 통신량 최소화를 목표로 하는 통신 회피 알고리즘(communication-avoiding CG, CA-GMRES 등)이 개발되고 있다.
이들 기법의 주요 아이디어는 선형 독립 벡터들을 한꺼번에 로컬로 구축한 뒤, 단일 통신 교환으로 향후 여러 번 쓸 정보를 미리 집약해 두는 것이다. 예를 들어 Krylov 방법에서 직교화(Arnoldi, Gram-Schmidt 등)를 수행할 때, 기존에는 단계마다 부분벡터를 브로드캐스트하거나 리듀스(reduce)하여 직교 조건을 갱신했다면, 통신 회피 알고리즘에서는 여러 스텝을 묶어서 한 번에 대규모 통신을 처리한다. 이렇게 통신 횟수를 획기적으로 줄이면, 고속 네트워크 환경에서도 확장성이 개선된다. 다만 한 번에 묶는 스텝 수가 늘어날수록 국소 메모리 사용량이 커질 수 있고, 직교화 오차 누적이 발생하기 때문에 적절한 절충이 필요하다.
파이프라인화(pipelining) 접근
반복법의 잔차 갱신과 내적, 직교화 등의 단계가 순차적으로 일어나는 전통적 구조에서는, 다음 연산을 시작하기 전까지 통신이 완료되어야 하므로 CPU/GPU 코어가 빈 대기 시간을 갖는 경우가 많다. 파이프라인화 기법은 이러한 단계를 가능한 한 겹쳐서 진행함으로써, 계산과 통신을 동시에 수행하도록 스케줄링한다. 예를 들어 스레드 일부는 이전 단계에서 계산된 벡터를 GPU로 전송하거나, 다른 노드에 통신을 수행하고, 나머지 스레드는 현재 단계의 계산을 지속한다.
이런 파이프라인 구조는 GPU·CPU 혼합 환경에서 특히 중요하다. 각 반복 단계에서 데이터를 GPU로 보내 곱셈을 수행하고, 그 결과를 CPU로 돌려받아 집계하는 과정을 순차적으로 처리하면 오버헤드가 크지만, 파이프라인을 잘 구성하면 통신과 계산을 상당 부분 겹칠 수 있다. CUDA 스트림(stream) API나 ROCm HIP 스크림 등 저수준 기능을 적극 활용해야 하며, 통신 라이브러리(MPI, UCX, libfabric 등)와의 상호작용도 파이프라인에 맞추어 튜닝해야 한다.
mermaid를 이용한 병렬 솔버 워크플로우 예시
아래는 Krylov 반복 해법을 대규모 병렬 시스템에서 수행하는 과정을 개략적으로 나타낸 다이어그램이다. 스파스 행렬 $\mathbf{A}$와 벡터 $\mathbf{b}$가 노드별로 분산 저장되어 있다고 가정한다.
이 다이어그램은 Krylov 반복에서 반복적으로 호출되는 SpMV(행렬-벡터 곱), 국소 계산, 전역 통신을 순환 구조로 묘사한다. 파이프라인 기법이나 통신 회피 알고리즘을 적용하면, B와 C, D 사이의 통신/동기화 단계를 최소화 혹은 중첩하여 전체 시간을 단축할 수 있다.
체크포인팅(checkpointing)과 고장 회복(fault tolerance)
규모가 큰 HPC 환경에서는 수많은 노드 중 일부가 계산 도중 장애를 일으킬 위험이 상존한다. 대규모 반복 해법이 진행되는 과정에서 장애가 발생하면, 초기에 수집한 해를 전부 잃어버릴 수도 있으므로, 주기적 체크포인팅과 결합된 복원(fault recovery) 매커니즘이 중요하다. 반복해법은 순차적 계산을 거듭하여 잔차를 줄여 가는 구조이기 때문에, 특정 지점에서의 중간 해($\mathbf{x}_k$)와 잔차 정보, 사전조건 변수 등을 안전한 저장소에 저장해 두는 방식으로 대비한다.
체크포인트 파일을 자주 생성하면 I/O 부하가 커지고 계산 흐름이 자주 중단되어 성능 저하가 발생할 수 있으므로, 장애 확률과 성능 감소 간 균형을 맞출 주기가 필요하다. 자가 복구(self-stabilizing) 기법이나 ABFT(Algorithm-Based Fault Tolerance) 기법을 통해, 핵심 반복 단계에서 단순한 검증(예: Checksums)을 수행하고 오류 검출 시 국소적으로 보정하는 접근도 연구되고 있다.
전처리(reordering)와 후처리
스파스 행렬에서 fill-in을 감소시키거나 대각 우세 구조를 강조하려면 행·열 순서를 적절히 바꾸는 것이 효과적이다. 대칭 양의정부 행렬의 경우 RCM(Reverse Cuthill-McKee)나 METIS/Scotch 같은 그래프 파티셔닝 알고리즘을 통해 밴드폭을 줄이거나, 도메인 분할에 적합한 구조로 변환한다. 이 작업을 전처리 단계에서 수행하면, 직접 해법이나 ILU 계열 사전조건에서 발생하는 불필요한 fill-in이 감소하거나, 반복 수렴이 빨라진다.
문제에 따라 스케일링(scaling) 작업도 선행될 수 있다. 행렬 $\mathbf{A}$의 행이나 열별로 크기가 매우 다르다면, 일련의 대각 행렬 $\mathbf{D}_r, \mathbf{D}_c$ 등을 사용해
형태로 재구성해 두면, $\mathbf{A}$의 비편향적(better-conditioned) 형태가 확보되어 반복법 수렴에 긍정적인 영향을 미친다. 단, 해 $\mathbf{x}$를 최종적으로 복원하기 위해 $\tilde{\mathbf{x}}$에서 역스케일링 과정을 거쳐야 한다.
이처럼 대규모 선형계 문제를 풀 때에는 전·후처리, 하이브리드 사전조건, 병렬·분산 알고리즘, 로우랭크 근사, 통신 최적화, 체크포인팅 등을 종합적으로 고려하여 효율과 안정성을 함께 확보해야 한다.
Newton-Krylov 접근과 비선형 문제의 대규모 해법
물리·공학 현상에서는 비선형 편미분 방정식이나 복잡한 경계 조건으로 인해, 단순한 선형계가 아닌 대규모 비선형 문제를 다루어야 하는 경우가 많다. 이를 풀기 위해 Newton-Krylov 계열 기법이 널리 활용된다. 즉, 고전적 뉴턴 방법에서 야코비 행렬(또는 그 근사)을 구한 뒤 선형 보조문제
을 매 뉴턴 스텝마다 풀어 $\delta \mathbf{u}_k$를 구하고, 이를 통해 해를 갱신한다. 여기서 $\mathbf{J}(\mathbf{u}_k)$는 $\mathbf{F}$의 야코비 행렬이며, $\mathbf{F}(\mathbf{u}) = \mathbf{0}$ 형태의 비선형 문제를 대상으로 한다.
대규모 문제에서 $\mathbf{J}(\mathbf{u}_k)$가 매우 커질 수 있고, 해당 행렬이 스파스 구조를 지니더라도 직접 분해로는 감당하기 힘든 연산량이 발생하기 쉽다. 따라서 뉴턴 반복 각 단계에서 보조 선형계
를 완벽히(직접) 풀지 않고, Krylov 하부공간 반복법(CG, GMRES 등)을 사용하여 근사해를 구하는 방식을 적용한다. 이를 Inexact Newton 방식이라고 부른다. 선형계 풀이 과정에서 사전조건을 잘 설계하면, 빠른 수렴과 함께 뉴턴 단계의 전체 횟수도 크게 줄어든다.
야코비-벡터 곱 $\mathbf{J}(\mathbf{u}_k),\mathbf{v}$는 형식미분(analytic differentiation) 대신 수치 도함수(finitedifference)나 색칠 알고리즘(coloring technique)을 이용해 효율적으로 계산하기도 한다. 특히 자동 미분(AD, Automatic Differentiation) 도구를 활용하면, 복잡한 다항식 혹은 프로그래밍 로직에 대해 기호적 혹은 소스 변환 방식으로 야코비을 도출하여 계산 과정을 자동화할 수 있다. 대규모 분산 환경에서는 야코비-벡터 곱을 노드별로 분산 처리하며, GPU 병렬화를 병행해 성능을 극대화한다.
Newton-Krylov 알고리즘에서 수렴 기준은 보통 비선형 잔차 $|\mathbf{F}(\mathbf{u}_k)|$가 일정 임계값 이하가 되는지를 확인하며, 선형 보조문제 해법의 정확도(내부 반복 수)와 외부 뉴턴 단계 수 사이에 균형이 중요하다. 즉, 선형계 해를 지나치게 정밀하게 구하려고 하면 불필요한 계산 부담이 늘고, 반대로 너무 느슨하게 풀면 뉴턴 단계가 많이 소요되어 결과적으로 총 계산 시간이 길어질 수 있다. 이를 위해 Inexact Newton에서는
같은 동적인 공차(dynamic tolerance) 설정을 두어, 초기에는 선형계 해법을 비교적 부정확하게(저정밀도) 풀다가 비선형 잔차가 줄어들수록 해법 정밀도를 높이는 방식을 취할 수 있다.
견고한 사전조건과 잔차 기반 적응
대규모 비선형 시스템을 풀 때는, 야코비이 바뀔 때마다 사전조건도 함께 갱신하거나 적어도 부분적 업데이트를 진행해야 한다. 완전한 ILU 분해나 멀티그리드 구조를 매 뉴턴 단계마다 재구성하는 것은 비용이 막대하므로, Schur 보충이나 도메인 분할 사전조건을 점진적으로 갱신하거나, AMG 프로롤롱/리스트릭션 연산자를 이전 단계 결과를 활용해 재활용하는 접근이 시도된다.
비선형 반복에서 야코비이 변동되면, Krylov 해법의 수렴 특성도 달라질 수 있다. 따라서 사전조건 역시 잔차나 단계별 조건수 변화를 모니터링하여, 필요할 때만 다시 구성하는 잔차 기반 적응(residual-based adaptation) 방식이 연구되고 있다. 이로써 매 단계에서 과도한 사전조건 구성 시간을 쓰지 않고도, 적정 수준의 수렴 가속 효과를 유지할 수 있다.
Matrix-Free 방법
Newton-Krylov 접근에서 야코비 행렬을 명시적으로 구성하지 않고, 벡터 곱 $\mathbf{J}(\mathbf{u}_k),\mathbf{v}$만을 직접 계산하는 Matrix-Free 방법이 자주 사용된다. 예를 들어 유한요소 해석에서 $\mathbf{F}(\mathbf{u})$를 구성하는 잔차 연산을 코드로 구현해 두었다면, $\mathbf{J}(\mathbf{u}_k),\mathbf{v}$는
꼴의 수치 근사로 얻거나, 또는 자동 미분을 이용해 탄젠트 연산 형태로 직접 산출할 수 있다. 이 경우 야코비 행렬이 매우 크고 스파스 패턴이 복잡한 문제에서도, 실제로 행렬 요소를 전부 메모리에 보관하지 않아도 되므로 계산과 저장 자원을 크게 절약할 수 있다.
그러나 Matrix-Free에서는 사전조건 역시 야코비이 없어도 구성 가능한 형태여야 하므로, 단순 대각 스케일링 사전조건부터 제한된 형태의 ILU나 멀티그리드 기법을 결합하는 등 제약이 따른다. 대규모 3D 유한요소 해석에서 이 방법을 채택하면, 물리 모델의 local assembly를 통해 필요한 국소 행렬-벡터 연산만 처리하고, 전역 행렬 형태는 직접 저장하지 않는 식으로 병렬 메모리 사용량을 최소화할 수 있다.
Newton-Krylov-Schwarz(NKS) 기법
도메인 분할법과 Newton-Krylov를 결합한 NKS 방법은, 비선형 PDE 문제를 분산 메모리 환경에서 대규모로 풀기 위한 대표적 전략이다. 전역 시스템에 대한 뉴턴 반복을 수행하면서, 각 뉴턴 단계에서 발생하는 보조 선형계는 Krylov 반복을 적용하고, 그 사전조건으로는 도메인 분할에 기반한 Schwarz 유형 사전조건을 쓰는 식이다. 도메인 내의 부분영역별로 국소 뉴턴해를 구하거나 국소 선형계 해법을 직접 수행하고, 전역적으로는 Krylov 반복을 통해 경계 조건 정합을 맞춰가는 구조가 전형적이다.
NKS 접근의 효율성은 서브도메인 크기(국소 문제 스케일), 중첩(overlapping) 정도, 경계 연결 전략 등에 따라 달라진다. 물리적으로 큰 기여를 하는 부분영역에 대해 좀 더 정교한 해법(정확한 직접 해법 등)을 적용하고, 비교적 영향이 작거나 균일한 구간에는 간단한 반복 해법을 쓰는 식으로 리소스를 차등 배분하기도 한다. 이를 통해 전체 계산 과정에서 통신량을 억제하고, 국소 해법의 병렬 효율을 높여 전반적인 스케일링을 달성한다.
Block Nonlinear 방법과 PFASST 등
비선형 해석에서 시공간을 블록으로 확장하여, 시간 차원을 포함한 병렬화를 강화하는 연구들도 있다. 예를 들어 다단계 병렬 시간 적분 알고리즘인 PFASST(Parallel Full Approximation Scheme in Space and Time)는 멀티그리드 아이디어를 시간 차원에도 적용하여, 한 번의 반복 과정에서 여러 시점(time-step)을 병렬로 계산하며 전역 잔차를 동시에 줄이는 전략을 취한다. 대규모 병렬 시스템에서 공간뿐 아니라 시간 차원까지 분할하면, 고전적인 순차적 시간 적분 방식의 병목을 완화할 수 있다.
이런 Block Nonlinear 접근은 PDE 문제의 시공간 해석, 비선형 다물리耦合(multi-physics coupling), 과도 문제(transient problems) 등에서 응용 여지가 크지만, 구현 복잡도가 높고 스테빌리티 분석이 까다롭다. 그러나 초고성능 컴퓨팅 시대에 성능 한계를 극복하기 위한 유망한 시도로 계속 진화 중이다.
알고리즘·컴파일러·하드웨어 공동 최적화
최근에는 CPU 아키텍처가 멀티·매니코어화, GPU 병렬화, HBM(High-Bandwidth Memory) 채택, 캐시 계층 변화 등으로 급속히 변화하고 있어, 대규모 행렬 연산 알고리즘 역시 하드웨어 특성에 맞추어 미세 조정이 필수적이다. 예를 들어 멀티그리드나 Krylov 반복에서 벡터화(vectorization)에 특화된 BLAS 루틴을 얼마나 활용할 것인지, GPU 텐서 코어(Tensor Core)의 매트릭스 연산 능력을 혼합 정밀도로 어떻게 극대화할 것인지, 통신 병목을 줄이기 위해 RDMA(Remote Direct Memory Access) 기능을 어떤 방식으로 사용할 것인지 등 미세 조율이 이뤄진다.
또한 컴파일러 수준에서 LTO(Link Time Optimization), IPO(Interprocedural Optimization), 자동 벡터화, GPU 커널 병합 등 다양한 최적화 기법이 적용될 수 있다. 특정 플랫폼이나 베어 메탈 환경에서 최적화된 라이브러리를 활용하는 것이, 범용 라이브러리를 쓰는 것보다 훨씬 나은 성능을 보장할 때가 많다. 따라서 대규모 수치 해석 프로젝트에서는 알고리즘 설계, 병렬 프로그래밍, 컴파일·라이브러리 튜닝, 클러스터 설정 등 각 단계를 유기적으로 결합하여 최적 성능을 이끌어내야 한다.
자동 미분(AD)과 모델링 프레임워크
대규모 선형계나 비선형 계를 풀 때, 적절한 야코비(혹은 그 근사)의 정확성과 계산 효율을 동시에 확보하는 것이 매우 중요하다. 주어진 지배방정식과 경계/초기 조건으로부터 잔차 벡터 $\mathbf{F}(\mathbf{u})$를 형성하는 과정은, 일반적으로 수치 해석 코드 상에서 다양한 보조 변수를 다루게 된다. 물리학, 기계공학, 전산유체역학(CFD), 전자기 해석 등 여러 응용 분야에서는 보통 수십~수백만 개 이상의 자유도를 가지는 시스템을 구성하며, 그 과정에서 방정식의 복잡성으로 인해 기호적 미분을 직접 수행하기 어려울 때가 많다.
이럴 때 자동 미분(Automatic Differentiation, AD) 기법이 유용하다. 소스 변환(Source Transformation) 방식의 AD는 주어진 코드(예: C, C++, Fortran 등)를 파싱하여 미분 연산을 삽입한 새 코드를 자동 생성하고, 역전파(Reverse Mode)나 순전파(Forward Mode)를 통해 야코비 항을 효율적으로 추적한다. 연산 그래프(Computation Graph) 수준에서의 오버로드(Operator Overloading) 기반 AD 역시 C++ 템플릿 메타프로그래밍 등을 사용하여 구현할 수 있다. 이를 통해 $\mathbf{J}(\mathbf{u}_k)$가 명시적 스파스 구조로 만들어지거나, 적어도 $\mathbf{J}(\mathbf{u}_k),\mathbf{v}$ 형태의 곱셈이 자동으로 지원된다. Matrix-Free 해법을 구현할 때, 특히 역전파형 AD가 빠른 값을 보이는 경우가 많다.
복잡한 다물리 시뮬레이션 환경에서는 전산유체(CFD), 열전달, 구조해석(고체역학), 반응물질 이동 등 서로 다른 물리 모델들이 결합하여 커다란 비선형 시스템을 형성한다. 이런 멀티피직스(multi-physics) 시나리오에서 매번 손수 미분 방정식을 코드화하고 야코비을 구성하는 것은 개발·유지보수 비용이 매우 높다. 따라서 고수준 모델링 프레임워크(FEniCS, deal.II, FreeFEM, Sundance 등)나 상용 PDE 솔버(ANSYS, COMSOL 등)의 내부 루틴을 통해, 멀티피직스 연산자와 야코비 조립을 자동화하는 방법을 적극 검토한다. 이때에도 내부적으로는 간접적 자동 미분 기술이 활용되어, 반복해법에 필요한 열·행 계산을 투명하게 제공한다.
p-적응도와 분할 정렬(Adaptive Refinement)
유한요소법(FEM) 기반의 대규모 문제에서, 단순히 격자(grid)의 간격을 세분화(h-적응)하는 것만이 아니라, 요소 차수를 높여(p-적응) 정밀도를 높이는 전략도 수행된다. 실제로 요소 차수가 달라지면 지역적으로 스텐실(stencil) 구조가 변동되고, 행렬 내의 연결 패턴과 계수 값이 크게 달라질 수 있다. 따라서 초기 문제 설정에서 하위 차수(예: 1차 또는 2차)로 시작했다가, 특정 지역에서의 해 오차가 크게 추정되면 그 영역만 고차 다항(예: 3차 이상)으로 올려 메쉬를 재구성하기도 한다.
이런 적응 과정(hp-adaptive 등)에서 행렬 구조가 반복적으로 변경될 때, 사전조건이나 도메인 분할도 함께 재설정해야 하므로 계산이 복잡해진다. 일괄적 리팩터링으로 대규모 LU 분해를 매번 다시 하는 것은 비현실적이어서, 이전 단계에서 얻은 행렬 구조와 해 정보를 기반으로 차근차근 업데이트하는 점진적(adaptive) 알고리즘 설계가 필요하다. 블록 구조나 멀티그리드 레벨 구성을 기반으로 하면, 각 레벨에서 적응된 그리드 정보를 재활용하는 스킴을 고려할 수 있다. 예컨대 코어 영역은 이미 고차 요소로 충분히 해상도가 커서 더 이상 세분화할 필요가 없고, 경계 근방만 세분화 혹은 고차화를 수행하여 전체 행렬의 크기와 랭크를 제어한다.
특이행렬과 준특이 문제
물리적으로 시스템이 강체 운동(rigid body motion)을 허용한다거나, 경계 조건이 불충분하여 행렬 $\mathbf{A}$가 정칙이 아닌 경우(즉, 특이(singular)하거나 준특이(near-singular)한 경우)가 발생하기도 한다. 이때는 $\mathbf{A}\mathbf{x} = \mathbf{b}$ 해가 유일하지 않거나 아예 존재하지 않을 수 있다. 실제로 큰 스파스 행렬 문제에서 행 랭크(행 공간)와 열 랭크(열 공간)가 일치하지 않는 상황은 자주 나타난다.
이런 특이·준특이 문제에 대해서는 축소 문제(reduced system)나 구속 조건을 추가로 도입하여, 계수 행렬을 유의미하게 invertible한 형태로 변환해야 한다. 예를 들어 물리적으로 강체 운동 방향에 대응하는 몇 개의 자유도를 고정(Dirichlet 구속)해 주거나, 라그랑주 승수(Lagrange multiplier)를 삽입하여 경계나 점 구속식을 추가하는 식이다. 또는 적절한 정규화(regularization) 항을 부여해 행렬이 완만하게 양의정부 구조를 띠도록 유도하기도 한다.
반복 해법 입장에서는 준특이도가 클수록 조건수가 커져 수렴이 느려지므로, 가중 least-squares 해법이나 투영 기법을 조합하여 적절한 해를 찾을 수 있도록 설계해야 한다. 직교 투영이나 미분 기법을 사용해 영공간(null space)을 분리해 두고, 나머지 부분에서만 해법을 진행하는 방식도 연구되고 있다.
희소 행렬의 구조적 데이터베이스
기초 연구나 알고리즘 검증 단계에서, 다양한 형태의 희소 행렬에 대해 알고리즘 성능을 시험해 보고자 할 때는, well-known 스파스 행렬 저장소(예: SuiteSparse Matrix Collection, formerly UF Sparse Matrix Collection)를 참조한다. 이 데이터베이스에는 그래프 문제, 유한요소 문제, 최적화 문제 등 여러 응용 분야에서 추출된 실제 희소 행렬들이 방대하게 저장되어 있어, 해당 행렬에 대한 차원, 비율(sparsity ratio), 대칭성, 조건수 추정치, 스트럭처 종류 등이 함께 제공된다.
이런 벤치마크 행렬을 사용하면, 특정 반복 해법이 어떤 유형의 문제에서 수렴이 빠른지 혹은 느린지, ILU나 AMG 사전조건이 어느 스파스 구조에 유리한지 등을 쉽게 비교 실험할 수 있다. 또한 동일 행렬에 대해 GPU 가속을 적용했을 때와 CPU만 썼을 때의 성능 차이를 정량적으로 측정해 볼 수도 있다. 대규모 시뮬레이션 프로젝트에서는 이처럼 다양한 테스팅 과정을 거쳐, 최적의 해법 조합(사전조건 + 반복자 + 병렬화 기법)을 찾아내는 것이 일반적이다.
정량 평가와 메트릭
병렬 선형 솔버의 성능을 평가하기 위한 대표적 지표로는 벽시계 시간(wall-clock time), 플롭스(FLOPS), 메모리 사용량, 부하 불균등도(load imbalance), 통신량, 병렬 효율성 등이 있다. 단순히 중앙처리장치(CPU) 스레드와 노드를 늘렸을 때 전체 수행 시간이 어떻게 변하는지를 보는 강한 스케일링(strong scaling) 테스트와, 문제 크기도 함께 늘려 가며 처리 시간을 확인하는 약한 스케일링(weak scaling) 테스트가 모두 중요한 척도가 된다.
대규모 HPC 시스템에서, 이론상 최대 부동소수점 연산 속도의 일정 비율 이상을 지속적으로 달성하는 것은 쉽지 않다. 실제로는 메모리 접근 패턴이 병목을 일으키거나, 노드 간 통신이 잦아지면 프로세서가 대부분의 시간을 대기 상태로 보내게 된다. 따라서 다양한 프로파일링 툴(Scalasca, Intel Trace Analyzer, NVProf, Nsight Compute 등)을 병행하여, 어떤 부분(행렬-벡터 곱, 내부 직교화, 통신, I/O 등)이 병목이 되는지를 파악하고 각 부분의 최적화를 시도한다.
GPU 환경에서는 코어별로 스레드를 얼마나 효율적으로 배치했는지, 워프(divergence) 문제는 없는지, 공동 메모리(shared memory)나 텍스처 메모리를 잘 활용하는지 등이 관건이 된다. 혼합 정밀도 연산을 쓸 경우, 텐서 코어나 SIMD 레지스터 구조를 잘 활용했는지도 점검한다. 이런 일련의 튜닝 과정을 지속해야, 대규모 선형계 풀이에서 안정적이고 빠른 솔루션을 얻을 수 있다.
추세와 미래 방향
대규모 행렬 연산 기법은 계속 진화하고 있으며, Exascale 시대의 HPC 환경에서는 수많은 프로세서와 가속기(GPU, TPU, FPGA 등)의 협업을 고려해야 한다. GPU나 신경망 NPU(neural processing unit) 같은 특수 코어에서 행렬 연산을 초병렬로 수행하는 환경이 일반화될수록, 혼합 정밀도 기법이나 통신 회피 알고리즘 같은 고차원적 최적화가 더욱 필수적인 요소로 자리 잡는다. 또한 양자 컴퓨팅(quantum computing)의 가능성이 제기되는 현 시점에서, 전통적 디지털 HPC와 양자 알고리즘을 어떻게 조합해 더욱 효율적인 대규모 행렬 솔루션을 구성할지도 흥미로운 연구 주제가 될 것이다.
PDE 제약 최적화와 대규모 KKT 시스템
전산과학 분야에서 대규모 행렬을 다루게 되는 또 다른 대표적 응용은 PDE에 의해 제약되는 최적화 문제다. 예컨대 다음과 같은 형태의 문제를 생각할 수 있다.
여기서 $\mathbf{u}$는 상태(state) 변수를, $\mathbf{p}$는 설계(design) 혹은 제어(control) 변수를 나타내며, $\mathbf{F}$는 PDE 및 경계 조건을 반영한 잔차 연산자다. 이 문제를 해석적 혹은 수치적으로 풀어내려면, 최적화 알고리즘(예: SQP, Sequential Quadratic Programming)이나 증분적 방법(예: Lagrange-Newton) 등을 통해 거대한 KKT 시스템을 풀어야 한다. KKT 시스템이란 다음과 같이 상태·접근방정식(Adjoint Equation), 그리고 최적성 조건을 한데 묶은 형태의 선형계 혹은 준선형계를 말한다.
$\boldsymbol{\lambda}$는 라그랑주 승수(접근 변수)이며, $\mathcal{L}(\mathbf{u},\mathbf{p},\boldsymbol{\lambda})$는 라그랑주 함수다. 문제 규모가 커지면 이 계수 행렬 또한 초대형 스파스 구조를 지니므로, 효율적인 선형 솔버가 핵심이다.
이런 PDE 제약 최적화에서 흔히 쓰이는 방법 중 하나가 상태·접근 기반의 축소 접근(state/adjoint approach)으로, $\mathbf{u}$와 $\boldsymbol{\lambda}$만으로 전방 PDE와 접근 PDE를 풀어 그 결과를 토대로 $\mathbf{p}$ 업데이트 방향을 결정한다. 이 경우에도 매 단계에서 거대한 선형계(혹은 비선형계)를 푸는 작업이 필요하기 때문에, Krylov 반복법, 도메인 분할, 멀티그리드, 혹은 사전조건 기법을 적절히 결합해 계산량을 줄여야 한다. 한편, full KKT 시스템을 직접 풀어버리는 approach도 가능하나, 여기서는 블록 분할·사전조건(예: Schur 보충 방식을 이용한 분해) 등 고급 설계가 필수적이다.
Adjoint 접근과 민감도 해석
설계 변수 $\mathbf{p}$의 차원이 크면, PDE 제약 조건을 만족시키며 목적함수 $\Phi(\mathbf{u}, \mathbf{p})$의 기울기 $\nabla_{\mathbf{p}}\Phi$를 구하는 과정이 계산적으로 큰 부담이 된다. 이때 접근(Adjoint) 방정식을 활용하면, 다수의 설계 변수가 있을 때에도 상태·접근 PDE 각 한 번씩만 풀어 전체 기울기를 효율적으로 얻게 된다.
즉, 상태 방정식
를 풀어 $\mathbf{u}^*$를 구하고, 이어 접근 방정식
를 풀어 $\boldsymbol{\lambda}$를 얻는다. 그 뒤
공식으로 목적함수 기울기를 한 번에 구할 수 있다. p\mathbf{p}가 수천~수백만 차원에 달하더라도, 접근 방정식을 한 번만 해결하면 되므로 효율성이 매우 높다. 문제는 상태·접근 방정식이 모두 대규모 스파스 선형계로 귀결된다는 점이며, 이를 해결하는 데에 앞서 언급한 고성능 반복 해법과 분산 메모리 계산 기법이 다시 등장한다.
역문제와 정규화 기법
PDE 제약 최적화는 역문제(Inverse Problem)의 범주로도 볼 수 있다. 실험이나 관측을 통해 얻은 데이터를 토대로, PDE 계수나 경계 조건, 초기 조건 등을 추정하려는 문제들은 보통 관측 데이터 오차, 측정 지점 불충분, 물리적 노이즈 등의 영향으로 인해 매우 ill-posed하다. 따라서 해가 유일하게 정해지지 않거나 수치적 불안정이 발생한다.
이를 극복하기 위해 Tikhonov 정규화, TV(Total Variation) 정규화, 또는 다른 형태의 제약 항(regularization term)을 목적함수에 추가하는 전략을 쓴다. 예를 들어
형태로 정규화 항 $|\mathbf{p}|^2$를 더해, 해의 불안정성을 완화시키고 일정한 매끄러움이나 소규모 성분을 유지하도록 유도한다. 그 결과 KKT 시스템이 양의정부(또는 준양의정부) 구조로 가까워지면서 선형 솔버의 조건수가 개선되고 수렴 속도가 빨라지기도 한다. 단, $\alpha$ 등 정규화 계수를 선택하는 과정에서는 수많은 모의실험 혹은 적응형(A-optimal, L-curve 등) 기법이 쓰여야 한다.
PDE 제약 최적화 소프트웨어 프레임워크
대규모 PDE 제약 최적화 문제를 풀기 위해 PETSc, Trilinos, FEniCS, deal.II, Firedrake 등과 같은 오픈소스 프레임워크들이 발전해 왔다. 이들에는 선형·비선형 솔버, 접근 알고리즘, 자동 미분, 멀티피직스 적응 등이 통합적으로 구현되어 있다. 예컨대 PETSc에는 TS(시간 적분) 모듈과 SNES(비선형 방정식 솔버) 모듈, 그리고 TAO(최적화) 모듈이 존재하여, 한 코드 베이스에서 PDE 해석, 비선형 해석, 최적화 과정을 모두 수행할 수 있다.
Trilinos에서도 NOX(비선형 해석), LOCA(고유치·분기해석), ROL(최적화 라이브러리)가 통합적으로 연동되며, Epetra/Tpetra로 표현되는 대규모 스파스 행렬·벡터 연산을 저수준에서 지원한다. FEniCS나 deal.II 같은 FEM 라이브러리는 문제 정의와 해석 과정을 고수준 언어(Python, C++ 등)로 간결하게 표현할 수 있게 해 주고, 내부적으로 자동 미분과 효율적 스파스 연산, 병렬화(MPI, OpenMP, CUDA 등)를 제어한다. 이런 툴들은 멀티피직스 coupled PDE 해석이나 다단계 비선형 반복, 그리고 대규모 KKT 시스템 풀이를 상위 레벨에서 편리하게 다루는 데 최적화되어 있다.
커스텀 병렬 구조의 접목
이렇게 잘 정돈된 프레임워크를 쓰지 않고, 프로젝트의 독자적 요구사항에 맞춰 직접 코드를 짜는 경우에도, 분산 메모리(MPI) + 멀티스레드(OpenMP) + GPU(CUDA/HIP) 혼합 구조에서 PDE 제약 최적화 알고리즘을 구현할 수 있다. 상태 방정식 풀이와 접근 방정식 풀이가 동일한 물리 도메인 위에서 이뤄진다면, 서로 같은 분산 데이터 구조(도메인 분할, 메쉬 분할)를 재사용할 수 있으며, 로컬 어셈블리(assembly) 과정과 전역 통신 패턴도 거의 동일하다. 이런 점을 활용하면 상태·접근 해석을 병렬적으로 동시에 수행하거나, 파이프라인 기법을 적용하는 등 추가적인 속도 향상이 가능하다.
다만 상태와 접근 해의 정확도, 그리고 $\nabla_{\mathbf{u}}\mathbf{F}$ 계산 과정이 꼬여 있어서 디버깅과 유지보수가 어렵고, 문제 규모가 커지면 미묘한 수치적 불안정성이 발생할 수 있으므로, 철저한 검증과 로그 추적이 뒤따라야 한다. 실제 운영 단계에서는 체크포인트/리스타트 기능을 갖춰, 중간 단계에서의 해·접근 벡터를 저장해 둔 뒤 장애가 발생해도 다시 이어 풀 수 있도록 대비한다.
다목적(Multi-Objective) 또는 확률적(Stochastic) 최적화
주어진 PDE 시스템에 대해 단 하나의 스칼라 목적함수만 다루는 것이 아니라, 여러 개의 서로 다른 목적함수(예: 비용·효율·안정성)를 동시에 고려해야 할 수도 있다. 이때는 Pareto 최적성 개념을 적용한 다목적(Multi-Objective) 최적화가 성립하며, 각 목적함수를 가중합하거나 ε-제약 방식을 쓰는 등 다양한 접근이 가능하다. 그 과정에서도 여전히 주요 연산은 대규모 스파스 선형계의 반복적 풀이이므로, 앞서 언급한 알고리즘적·병렬적 기법이 동일하게 핵심이 된다.
확률적(Stochastic) 혹은 불확실성 양을 반영하는 최적화 문제에서는, 불확실성 변수를 샘플링하여 몬테카를로(Monte Carlo) 방식으로 여러 번 PDE를 풀어 평균·분산·VaR(Value at Risk) 등 통계적 지표를 최소화하는 접근이 일반적이다. 그러나 샘플 개수가 매우 많아지면, 총 PDE 풀이 횟수가 방대해지므로 배치(batch) 방식 혹은 멀티레벨 몬테카를로, 확률적 뉴턴 방법, 멀티그리드 가속 기법 등을 복합적으로 적용해 계산 비용을 억제해야 한다.
이렇듯 대규모 선형계 솔버는 단순한 PDE 해석을 넘어, PDE 제약 최적화나 역문제, 확률적 모델링 등 고차원 응용의 기반 모듈로서 활용된다. 그리고 이 모든 확장이 “얼마나 빠르고, 안정적이고, 확장성 있게 $\mathbf{A}\mathbf{x}=\mathbf{b}$를 풀어낼 수 있는가”라는 문제로 직결된다는 사실이, 결론적으로 대규모 행렬 풀이 알고리즘의 중요성을 더욱 부각시킨다.
정리 및 추가 언급 사항
지금까지 살펴본 대규모 행렬 솔루션 기법들은 단순한 선형계 풀이에서부터, 비선형 PDE, 도메인 분할, 멀티그리드, 뉴턴-Krylov, 도메인 분할 사전조건, 멀티피직스 및 PDE 제약 최적화 등의 범주를 망라한다. 다양한 문제 설정과 하드웨어 환경에서 어떤 알고리즘이 최선인지 판단하려면, 다음과 같은 요소들을 종합적으로 고려해야 한다:
문제 구조: 스파스 패턴, 대칭성 여부, 양의정부 여부, 다물리耦合, 시간종속성 등
행렬 크기 및 분할된 도메인 수, 희소도와 fill-in 발생 양상
병렬 환경: CPU 멀티코어, GPU, 클러스터 노드 수 및 통신 성능, 메모리 대역폭
사전조건: 도메인 분할, ILU, AMG, 하이브리드 로우랭크 근사, Matrix-Free, DD(도메인 분할) + MG 등
알고리즘적 세부: Krylov 반복( CG, GMRES, BiCGSTAB, CA-GMRES ), 혼합 정밀도, 파이프라인, 통신 회피
소프트웨어: PETSc, Trilinos, hypre, FEniCS, deal.II 등 제공 프레임워크와의 호환성
튜닝: 캐시 블록 최적화, GPU 커널 설계, MPI 통신 레이아웃, 체크포인팅, 로깅 및 오류 복구
유지보수·개발 비용: 자동 미분(AD) 활용, 모델링 프레임워크, 사용자 커스텀 코드 혼합 여부
아래는 간단히 Octave를 이용해 비교적 큰 희소행렬을 반복 해법으로 푸는 스크립트를 예시로 든 것이다. 실제로는 수십만~수백만 차원의 스파스 행렬까지 가능하며, 적절한 사전조건 모듈 또는 병렬 버전 연동이 뒤따라야 한다.
이 예제는 단순한 스파스 대각선 행렬을 생성하여 GMRES 반복풀이를 시도한다.
flag = 0이면 수렴 성공
flag = 1이면 반복 횟수(maxIter)를 초과, etc.
현실적 대규모 문제에서는 도메인 분할(병렬 계산)이나 사전조건 등을 추가해야 하며, 메모리 요구사항도 훨씬 커지므로 적절한 HPC 환경에서 MPI 또는 CUDA 등이 연동된 라이브러리를 사용한다.
Last updated