# 희소 행렬(Sparse Matrix)과 효율적 연산 #### 희소 행렬의 개념과 의의 행렬의 크기가 매우 큰 문제를 다룰 때 대부분의 계수는 실제로 0이 아닌 항이 매우 적고, 나머지는 0인 경우가 많다. 이러한 행렬을 희소 행렬이라고 부른다. 즉, 전체 원소의 수에 비해 0이 아닌 원소의 수가 매우 적을 때 희소 행렬의 특성이 두드러진다. 이는 물리학, 컴퓨터 그래픽스, 대규모 공정 시뮬레이션, 전산 유체 역학(CFD) 등에서 공통적으로 나타나는 현상이다. 희소 행렬에 대해서는 0이 아닌 원소들이 어떻게 분포되어 있는지가 연산 및 메모리 사용량을 크게 좌우한다. 또한 희소 행렬을 위한 연산 기법을 설계하면, 실질적으로 필요한 연산의 양이 크게 줄어들고 메모리 사용량도 절약할 수 있다. 따라서 희소 행렬에 대한 효율적 표현과 연산 기법은 대규모 선형 방정식을 풀거나 고차원 문제를 다룰 때 필수적으로 고려되어야 한다. #### 희소 행렬의 주요 저장 방식 희소 행렬을 실제로 다룰 때에는 모든 원소를 $n\times n$ 형태로 그대로 저장하지 않는다. 필요 이상의 메모리를 차지할 뿐만 아니라, 0이 아닌 원소에 대한 정보만 추출하기가 어렵고 연산 속도 또한 비효율적이다. 따라서 주로 사용하는 희소 행렬의 대표적 저장 방식으로는 Coordinate(CoO) 형식, Compressed Row Storage(CRS 또는 CSR) 형식, Compressed Column Storage(CCS) 형식 등이 있다. 모든 형식의 공통된 목표는 0이 아닌 원소의 위치와 값을 효율적으로 추적하여 곧바로 연산에 이용하는 것이다. #### Coordinate(CoO) 형식 Coordinate 형식은 0이 아닌 항이 위치한 (행 인덱스, 열 인덱스)와 그 값을 한 쌍으로 묶어서 별도의 배열에 보관한다. 일반적으로 행 인덱스 정보를 저장하는 배열, 열 인덱스 정보를 저장하는 배열, 그리고 0이 아닌 원소의 실제 값을 저장하는 배열 총 세 개를 사용한다. $nnz$를 전체 0이 아닌 항의 개수라고 하면, 이 형식에서는 다음과 같은 배열이 필요하다. $$ \begin{align} \text{RowIndices} \in \mathbb{R}^{nnz} \\ \text{ColIndices} \in \mathbb{R}^{nnz} \\ \text{Values} \in \mathbb{R}^{nnz} \end{align} $$ 행렬 $\mathbf{A}$에서 0이 아닌 항을 $(i\_k, j\_k)$라고 하고 그 값을 $a\_{i\_k j\_k}$라고 하자. 그러면 위 배열들은 다음과 같은 정보를 담는다. $$ \text{RowIndices}\[k] = i\_k, \quad \text{ColIndices}\[k] = j\_k, \quad \text{Values}\[k] = a\_{i\_k j\_k} $$ 이러한 형태의 장점은 단순하고 직관적이라는 점이다. 특히 초기 값 설정이나 행렬의 요소가 임의로 추가 혹은 삭제되는 경우에 적합하다. 다만 행렬-벡터 곱 연산과 같은 연산 과정을 수행할 때, 데이터가 특정 행 또는 열을 기준으로 정렬되어 있지 않을 수 있으므로 추가적인 정렬이나 인덱스 처리가 필요하다. #### Compressed Row Storage(CRS 또는 CSR) 형식 CRS 형식은 대부분의 희소 연산 라이브러리에서 널리 채택되는 형식이다. 같은 행에 속하는 0이 아닌 항의 열 인덱스 및 값을 연속된 형태로 저장하고, 해당 행이 어디에서 시작되는지를 별도의 포인터 배열로 관리한다. 일반적으로 다음과 같은 세 개의 배열을 쓴다. $$ \begin{align} \text{RowPtr} \in \mathbb{R}^{n+1} \\ \text{ColIndices} \in \mathbb{R}^{nnz} \\ \text{Values} \in \mathbb{R}^{nnz} \end{align} $$ 여기서 $n$은 행렬의 차원(즉 행의 개수), $nnz$는 0이 아닌 항의 총 개수다. $\text{RowPtr}\[i]$에는 $i$번째 행의 0이 아닌 항이 $\text{ColIndices}$와 $\text{Values}$ 배열 상에서 어디서부터 시작되는지에 대한 인덱스 정보가 들어 있다. 예를 들어 $i$번째 행의 0이 아닌 항들은 $\text{ColIndices}\[\text{RowPtr}\[i]]$부터 $\text{ColIndices}\[\text{RowPtr}\[i+1]-1]$까지 연속 저장되며, 값들은 $\text{Values}\[\text{RowPtr}\[i]]$부터 $\text{Values}\[\text{RowPtr}\[i+1]-1]$에 해당한다. 행렬-벡터 곱 $\mathbf{y} = \mathbf{A}\mathbf{x}$을 CRS 형식으로 수행할 때 각 행별로 열 인덱스에 맞춰 적절히 곱셈 및 덧셈을 수행하면 되므로 매우 효율적이다. CRS 형식으로 행렬 $\mathbf{A}$를 표현한다면, $i$번째 행에 대한 결과 $\mathbf{y}\[i]$를 구하는 과정은 대체로 다음과 같이 전개된다. $$ \mathbf{y}\[i] = \sum\_{k=\text{RowPtr}\[i]}^{\text{RowPtr}\[i+1]-1} \text{Values}\[k],\mathbf{x}\[\text{ColIndices}\[k]] $$ 이는 $i$번째 행에 있는 0이 아닌 항들만 모아서 곱을 수행하는 방식이므로, 모든 행을 순차적으로 처리하면 된다. 이처럼 CRS 형식에서는 같은 행에 속한 원소들이 연속되어 있기 때문에 행을 기준으로 한 스칼라 곱이나 벡터 곱 등에 아주 적합하다. #### Compressed Column Storage(CCS) 형식 CRS 형식이 행을 기준으로 압축 저장했다면, CCS 형식은 열을 기준으로 압축 저장한다는 점만 다르다. 따라서 $\text{RowIndices}$, $\text{ColPtr}$, $\text{Values}$ 배열을 사용하며, $\text{ColPtr}\[j]$가 $j$번째 열에 속하는 0이 아닌 항이 어디서부터 시작되는지 가리킨다. 연산에서 열 위주의 접근이 필요한 경우, 예컨대 열을 기준으로 데이터 처리가 잦은 알고리즘에서는 CCS 방식을 채택한다. #### 희소 행렬의 시각적 이해 아래에는 $5\times 5$ 크기의 예시 행렬에서 0이 아닌 항들이 어떻게 배치될 수 있는지 mermaid를 이용해 개략적으로 표현해 보자. {% @mermaid/diagram content="graph TB A((A11)) --> B((A12)) A((A11)) --> D((A14)) B((A12)) --> C((A13)) C((A13)) --> E((A15)) D((A14)) --> E((A15)) " %} 행렬의 구체적인 계수 값을 나타낸 것은 아니지만, 어떤 위치에 0이 아닌 항이 몰려 있는지를 그림으로 이해할 수 있다. 실제로는 문제마다 특정 구조가 뚜렷하게 나타나기도 하며, 대칭 행렬인지 여부, 밴드 구조인지 여부, 블록 대각 형태인지 여부 등을 통찰할 수 있다. #### 희소 행렬 연산의 난이도와 특성 희소 행렬 연산에서 가장 큰 어려움 중 하나는 일반적인 Dense 형태와 달리, 연산 과정에서 새로운 0이 아닌 항이 발생하는 “fill-in” 현상을 제어하는 것이다. 예컨대 직교화 과정이나, LU분해, Cholesky 분해, QR분해 등의 직·간접 선형 연산을 할 때, 원래 0이었던 위치가 분해 과정에서 0이 아닌 항이 되어버리는 상황이 자주 발생한다. 이로 인해 예기치 못한 메모리 사용량 증가가 나타날 수 있으며, 연산 시간 또한 크게 늘어날 수 있다. 따라서 대규모 시스템을 다룰 때는 행렬의 특성을 파악하고, 행렬의 대칭성이나 대각 블록 구조, 대각 우세 여부 등을 최대한 활용해야 한다. 행렬의 모양을 재배열(즉 행과 열의 순서를 재정렬)하여 fill-in을 줄이는 기법이 크게 발전되어 왔다. 대표적으로 Cuthill-McKee 알고리즘이나 Minimum Degree 알고리즘 등이 존재하며, 이를 사용하여 희소 분해를 수행하는 과정에서 새롭게 생기는 0이 아닌 항의 수를 가능한 한 억제한다. 또한 주어진 희소 구조에 특화된 사전조건자(Preconditioner)를 구성하면, 반복법에서의 수렴을 빠르게 만들어준다. 따라서 희소 행렬의 구조적 특성을 최대한 이용하는 것이, 단순히 하드웨어 성능을 높이는 것 못지않게 중요하다. #### 대규모 시스템에서의 메모리 고려 사항 희소 행렬에 대한 계산을 수행할 때, $n$의 값이 매우 클 수 있으므로(수만, 수십만, 혹은 그 이상) 전체 메모리 사용량이 중요한 이슈가 된다. 일반적인 Dense 행렬은 $n^2$개의 원소를 그대로 저장해야 하므로 대규모 문제에서는 사실상 메모리에 적재조차 하기 어렵다. 그러나 희소 행렬 구조를 잘 파악하고 압축 저장을 사용하면, 실제 0이 아닌 항들만 관리하므로 $nnz$가 $n^2$에 비해 훨씬 작을 때에는 막대한 메모리 절감을 얻을 수 있다. 또한 계산 과정에서 발생하는 임시 버퍼나 중간 연산 결과가 차지하는 메모리도 잘 고려해야 한다. 예컨대 반복법을 활용하는 경우에는 매 반복 단계마다 $\mathbf{A}\mathbf{x}$ 형태의 연산이 일어나는데, 이를 효율적으로 수행하려면 희소 형식이 반드시 잘 정비되어 있어야 한다. 분해법을 사용하는 경우라면 중간 분해 과정에서 fill-in에 따른 추가 저장 공간도 합산해 보아야 한다. 그래서 재정렬 알고리즘으로 행렬 구조를 변경하거나, 적절한 분할 기법으로 문제를 여러 개의 블록으로 쪼개어 병렬화 전략을 쓰기도 한다. #### 희소 구조 활용의 실제 예 희소 행렬 구조를 적절히 활용하면, 계수가 매우 큰 방정식을 효율적으로 푸는 것이 가능하다. 앞으로 이어질 내용에서 구체적인 예시와 계산 알고리즘, 추가적인 최적화 기법 등을 더욱 자세히 살펴보겠다. 예를 들어 Octave나 Python(CSR 라이브러리), C++(Eigen, Trilinos 등) 같은 수치해석 라이브러리를 통해 코드 차원에서 구체적으로 어떻게 구현되는지 살펴볼 예정이다. 이때 각 라이브러리가 채택한 내부 저장 형식, 그리고 행렬-벡터 곱 수행 방식, 전처리기법, 재정렬 등에 대해서도 다룰 것이다. #### fill-in 현상과 재정렬 희소 행렬의 구조가 유지되려면 분해나 곱셈 같은 연산 과정에서 발생하는 fill-in 현상을 억제해야 한다. 예를 들어 LU 분해나 Cholesky 분해 과정에서 원래 0이었던 위치가 비(非)영이 항으로 채워지는 일이 생길 수 있는데, 이를 효과적으로 줄이지 않으면 연산 및 메모리 사용량이 기하급수적으로 늘어난다. 다음의 단순한 예시로 fill-in 개념을 살펴보자. 행렬 $\mathbf{A}$가 희소 구조를 가졌다고 하자. 어떤 LU 분해를 통해 $\mathbf{A} = \mathbf{L} \mathbf{U}$를 구한다고 할 때, $\mathbf{L}$과 $\mathbf{U}$의 특정 위치에 새로운 항이 생기는지 여부가 fill-in을 결정한다. 분해 과정의 어느 단계에서 $$ \begin{align} a\_{ij} \leftarrow a\_{ij} - \sum\_{k=1}^{j-1} l\_{ik} u\_{kj},\\ \text{(또는 이에 준하는 연산)} \end{align} $$ 형태의 갱신이 이뤄진다고 하자. 만약 이 갱신 결과가 0이 아니게 되면, 그 위치에 원래는 없던 항이 새로 생긴다. 이 항이 반복적으로 누적되면 희소 구조가 훼손되고, $nnz$의 증가는 곧바로 메모리와 연산량 증가로 이어진다. fill-in을 줄이기 위해서는 분해 시작 전에 행과 열의 순서를 재배열하여, 0이 아닌 원소들이 모여 있는 영역을 최대한 대각 근처로 모으거나, 기존의 0 위치가 계속해서 0으로 유지되도록 유도한다. 이런 절차를 “재정렬(Reordering)”이라고 부른다. 대표적으로 Cuthill-McKee(CM), Reverse Cuthill-McKee(RCM), AMD(Approximate Minimum Degree), Nested Dissection(ND) 등이 널리 사용되며, 특정 기법을 선택하는 것은 행렬의 구조나 대칭성 여부, 그리고 알고리즘 목적(직접 분해법 vs 반복법)에 따라 달라진다. #### Sparse Direct Solver와 Symbolic Factorization 희소 행렬에 대한 직접 분해법(Direct Solver)을 구현할 때는 대규모 행렬을 실제로 숫자 연산하기 전에 분해 과정에서 발생할 fill-in의 양을 미리 대략적으로 예측하고, 필요한 공간을 미리 할당하기도 한다. 이를 “Symbolic Factorization”이라고 하며, 실질적인 수치 연산(Numeric Factorization) 전에 행렬의 그래프 구조를 바탕으로 fill-in 위치를 탐색한다. Symbolic Factorization을 통해 분해 시 발생할 여유 공간을 사전에 확보해 두면, 분해 과정 중에 동적 메모리 할당을 최소화할 수 있고, 대규모 계산에서 반복적인 메모리 할당/해제를 피하게 되므로 성능에 이점이 있다. Nested Dissection이나 Minimum Degree 계열 재정렬 알고리즘은 이 Symbolic Factorization 과정에서 fill-in을 유의미하게 줄여 주므로, 대칭 스파스 행렬(특히 양의 정부호 행렬)에서는 일종의 표준 기법으로 자리잡았다. #### 희소 행렬과 반복법(Iterative Solver) 직접 분해법이 매우 큰 연산량이나 메모리 사용을 동반한다면, 반복법을 고려해야 한다. 반복법(Iterative Method)은 $\mathbf{A}\mathbf{x} = \mathbf{b}$를 풀기 위해 초기 추정치 $\mathbf{x}\_0$로부터 시작하여 점진적으로 해에 근사해 나가는 방식이다. CG(Conjugate Gradient), GMRES(Generalized Minimal Residual), BiCGSTAB(BiConjugate Gradient Stabilized) 등이 대표적인 예다. 이때 매 반복 단계마다 필요한 가장 중요한 연산은 행렬-벡터 곱 $\mathbf{A}\mathbf{x}$이다. 따라서 희소 구조가 제대로 압축돼 있지 않으면, 매번의 곱셈이 큰 비용을 유발하게 된다. 반면 CRS나 CCS 같은 형식을 사용하면 0이 아닌 항에 대해서만 곱셈을 수행하므로 반복법 전체의 계산 효율이 크게 높아진다. 희소 구조에 특화된 반복법에서는 사전조건자(Preconditioner) 또한 중요한 요소다. 대표적으로 ILU(Incomplete LU) 분해나 ILUT(ILU with threshold), IC(Incomplete Cholesky) 등이 쓰이는데, 이는 실제로 분해를 완전히 수행하지 않고 일부 fill-in을 허용하거나, 일정 수준 이하의 항은 버리도록 설계된다. 이로써 분해된 결과가 충분히 희소 구조를 유지하며, 반복법의 수렴 속도를 개선한다. ILU(0)는 0이 아닌 원소에 대해서만 분해를 수행하고 새로운 fill-in은 전혀 허용하지 않는 방식이며, ILU(k)는 대각선으로부터 k-밴드 이내의 fill-in은 허용하는 식이다. #### 병렬화와 분산 환경 희소 행렬 연산은 병렬화와 분산 환경에서 더욱 중요하다. 대규모 해석 문제나 초고차원 계 시스템을 다룰 때, 단일 노드의 메모리와 연산 성능만으로는 한계가 있기 때문이다. 희소 행렬을 행 단위 또는 블록 단위로 잘 분할(Partition)하여 여러 프로세서(또는 스레드)에 분산시키고, 각 단위끼리 중복 계산 및 통신을 최소화해야 한다. 행렬 $\mathbf{A}$가 대칭 스파스 형태라면, 분해 과정에서도 병렬화를 지원하는 여러 라이브러리(PARDISO, MUMPS, SuperLU, PETSc, Trilinos 등)를 사용할 수 있다. 병렬화 관점에서의 fill-in 제어는 더욱 중요한 의미를 갖는다. 불필요한 fill-in이 증가할수록 통신과 동기화 부담이 커지므로, 분산 메모리 환경에서 노드 간 메시지 전달량이 커진다. 따라서 재정렬 알고리즘은 단순히 fill-in을 줄이는 것뿐 아니라, 클러스터/그래프 파티셔닝을 통한 균등 분배와 통신 최소화가 함께 고려된다. Metis, Scotch, ParMetis, PT-Scotch 등은 이런 분산 환경의 파티셔닝 및 재정렬 라이브러리로 널리 알려져 있다. #### 블록 구조와 하이브리드 표현 일부 응용에서는 희소 행렬 중에서도 블록 단위로 희소한 형태를 갖는 경우가 있다. 예컨대 유한 요소법(FEM) 문제에서, 하나의 요소가 여러 자유도를 갖고 이들이 대응되는 부분을 하나의 블록으로 묶어 나타내면, 블록 대각 행렬(Block Diagonal Form)이나 블록 밴드(Block Banded Form)이 만들어진다. 이때 각 블록을 별도의 Dense 행렬처럼 다루고, 블록 간은 희소하게 연결하는 방식을 사용하면, 한편으로는 블록 내부에서의 효율적 연산(SIMD, 벡터화)을 활용하고, 다른 한편으로는 전체 구조적 희소성을 유지할 수 있다. 이런 하이브리드 표현을 사용하면, 특히 고차원 FEM 문제나 다중 물리(Multi-Physics) 연계 시뮬레이션에서 계산 효율과 메모리 효율을 동시에 얻을 수 있다. 예를 들어 CRS 형식을 확장하여 “Block Compressed Row Storage(BCRS)”를 구현하거나, 별도 라이브러리(예: PETSc의 BAIJ 형식)로 블록 스파스 구조를 관리할 수 있다. 각 블록은 작은 Dense 행렬이므로, BLAS 레벨-3 루틴을 적극 활용하여 연산에 대한 효율을 높일 수 있다. #### 숫자 안정성과 병합 희소 행렬 연산에서 숫자 안정성(Numerical Stability) 문제도 종종 부상한다. 예를 들어 ILU 분해나 Incomplete Cholesky 분해 시, fill-in을 인위적으로 제한하면 행렬의 스펙트럼 특성(예: 어떤 고윳값의 분포)이 극단적으로 변형될 수 있다. 이런 경우 반복법에서 수렴이 급격히 느려질 수 있고, 부동소수점 오차 누적으로 인해 불안정 현상이 발생할 수 있다. 그래서 threshold 기반의 ILU(ILUT) 기법이 개발되었고, 일정 임계값 이하의 fill-in만 버리도록 하여 숫자 안정성을 어느 정도 보장한다. 또한 CG나 GMRES 같은 방법을 수행할 때, $\mathbf{A}\mathbf{x}$의 희소 구조를 유지하려면 중간 결과와의 병합 과정에서 0이 아닌 항이 새로 생기지 않도록 관리해야 할 수도 있다. 그러나 실제 알고리즘은 중간 단계에서 희소 벡터끼리의 덧셈, 내적, 스칼라 배 등 다양한 연산이 수행되므로, 연산 결과가 완전히 동일한 희소 구조를 유지한다는 보장은 없다. 그래도 벡터 차원에서는 행렬에 비해 희소 구조가 단순하므로, 특별히 벡터 연산을 위한 별도의 압축 방식을 쓰거나, 동적으로 0을 제거하는(즉 tolerance 이하의 항을 0으로 취급하는) 방식을 쓸 때도 있다. #### 실제 코드 구현 예시(C++ / Eigen 라이브러리) 아래 예시는 C++과 Eigen 라이브러리를 사용하여 희소 행렬을 정의하고, 특정 연산을 수행하는 과정을 간단히 보여준다. Eigen은 내부적으로 CRS 방식에 가까운 Compressed Storage 방식을 사용한다. ```c++ #include #include int main() { // 행렬 크기 int n = 5; // 희소 행렬 정의 (Eigen에서 row-major, col-major 선택 가능) Eigen::SparseMatrix A(n, n); // 삽입 방식 A.insert(0,0) = 10.0; A.insert(0,1) = 2.0; A.insert(1,0) = 3.0; A.insert(2,2) = 5.0; A.insert(3,4) = 1.0; A.insert(4,3) = -1.0; // 압축(실제 내부 구조 압축) A.makeCompressed(); // 임의의 벡터 정의 Eigen::VectorXd x(n), b(n); x.setOnes(); // 행렬-벡터 곱 b = A * x; std::cout << "Result: " << b.transpose() << std::endl; return 0; } ``` 이 코드는 크기가 5×5인 희소 행렬 $\mathbf{A}$에 0이 아닌 항을 몇 개 삽입하고, $\mathbf{A}\mathbf{x}$ 연산을 수행한다. $A.makeCompressed()$ 호출을 통해 내부적으로 $\text{RowPtr}$, $\text{ColIndices}$, $\text{Values}$에 해당하는 압축 구조가 확정된다. 입력 단계에서 원소를 삽입할 때 중복 항이 발생할 수도 있는데, Eigen은 이들을 합산 처리하거나 예외를 일으킬 수도 있으므로 API를 사용할 때 주의한다. 위 예시에서처럼 실제 희소 행렬을 다룰 때는, 0이 아닌 항을 어느 시점에 입력하고, 언제 압축(혹은 스토리지 확정)을 할지가 중요하다. 입력이 잦은 상태에서 잦은 압축을 수행하면 오히려 오버헤드가 커질 수 있다. 반면 모든 항을 한 번에 삽입한 다음, 최종적으로 한 번만 압축하는 방식을 쓰면 효율적이다. #### GPU 가속과 SpMV 대규모 선형계 문제를 GPU(Graphic Processing Unit) 환경으로 옮기면 행렬-벡터 곱(SpMV: Sparse Matrix-Vector multiplication)을 병렬로 가속할 수 있다. GPU는 일반적으로 메모리 대역폭이 CPU에 비해 훨씬 크고, 대규모 스레드를 병렬로 실행할 수 있으므로, 희소 행렬 연산에서 상당한 이점을 얻을 수 있다. SpMV를 GPU에서 최적화하려면, 우선 행렬을 CRS(또는 CSR) 형식으로 옮긴 뒤, $i$번째 행의 시작·끝 인덱스를 빠르게 접근해 $\mathbf{x}$의 요소와 곱을 수행해야 한다. GPU에서는 스레드 하나가 행 하나의 연산을 담당하게 구성하는 방식을 자주 쓴다. 즉, 각 스레드는 해당 행에 있는 0이 아닌 항들을 순회하여 곱셈·덧셈을 하고, 그 결과를 $\mathbf{y}\[i]$에 저장한다. 이렇게 하면 큰 차원의 경우 상당히 높은 처리량을 얻을 수 있다. 한편 희소 행렬이 불균일하게 분포되어 있으면, 어떤 행에는 0이 아닌 항이 많고 어떤 행에는 매우 적을 수 있다. 이런 경우 GPU 스레드 간의 부하가 크게 편차가 날 수 있으므로, Warp divergence나 Load imbalance가 발생할 수 있다. 이런 현상을 완화하기 위해서는 행을 재정렬하거나, ELLPACK(ELL) 방식 혹은 Hybrid(ELL+COO) 방식을 선택하는 등 다양한 기법이 연구되어 왔다. CUDA, HIP, OpenCL 등 다양한 플랫폼에서 Sparse BLAS 루틴을 제공하고 있으므로, 이를 이용하면 SpMV 구현 시 편의를 누릴 수 있다. #### ELLPACK(ELL) 형식과 Hybrid 형식 ELL 형식은 모든 행이 갖는 0이 아닌 항의 최대 개수를 $k$라 하고, 이를 기준으로 $n\times k$ 형태의 배열을 만들어 각 행마다 실제로 존재하는 0이 아닌 항과 열 인덱스를 저장한다. 나머지 부분은 단순히 0으로 패딩한다. 이렇게 고정된 열 방향 배열로 만들어놓으면, GPU에서 행 단위의 연산이 매우 간단해진다. 그러나 행 별 0이 아닌 항 수의 편차가 클수록 패딩으로 인해 메모리 낭비가 심해질 수 있다. 이 때문에, ELL 형식을 부분적으로만 사용하고, 나머지는 COO 형식 등으로 보관하여 Hybrid 형식을 구성하기도 한다. 예를 들어 수많은 행에서 0이 아닌 항이 일정 수준 이하로 안정적으로 분포해 있으면, 그부분은 ELL로 관리하고, 예외적으로 더 많은 항을 갖는 행은 별도의 COO 영역으로 관리한다. 이를 통해 ELLPACK의 단순한 구조를 유지하면서도, 메모리 낭비와 부하 불균형 문제를 어느 정도 완화할 수 있다. #### GPU 병렬화 시 고려 사항 GPU에서 희소 연산을 수행할 때 고려해야 할 사항이 다양하다. 먼저, $\mathbf{A}\mathbf{x}$에 사용되는 벡터 $\mathbf{x}$와 결과 $\mathbf{y}$는 글로벌 메모리에 있어야 하며, SpMV 연산 동안 여러 스레드가 동시에 접근하게 된다. 이때 $\mathbf{y}$를 작성(Write)하기 위한 경쟁(Write conflict)은 발생하지 않는 편이지만, $\mathbf{x}$를 읽어오는 과정에서 캐시 동작 패턴을 최적화하려면 열 인덱스(ColIndices)가 정렬되어 있는 편이 유리할 수 있다. CRS 형식에서도 GPU 최적화를 위해 행별로 열 인덱스를 오름차순으로 정렬해 두는 경우가 많다. 또한, 분산 메모리 구조에서 GPU 여러 대가 연계되어 스파스 연산을 수행하면, 노드 간 통신(예: MPI)과 GPU 내의 병렬화 방식을 결합하는 복잡한 환경이 된다. 이때 희소 행렬을 어떻게 분할하고, 각 노드에서 관리하는 부분행렬(또는 부분열 행렬)을 어떻게 효율적으로 전송·결합할지가 중요한 이슈다. SpMV 연산 후에는 파편화된 결과를 다시 모아서 전체 $\mathbf{y}$ 벡터를 복원해야 하므로, 노드 간 통신량을 최소화하기 위해 파티셔닝 알고리즘을 신중히 적용한다. #### 희소 BLAS와 라이브러리 활용 Dense BLAS 표준처럼, 희소 선형대수 연산에 대해 표준화된 Sparse BLAS 라이브러리도 활발히 사용된다. 예를 들어 Intel MKL, cuSPARSE(NVIDIA), oneAPI, AMD HIPSPARSE 등의 라이브러리는 CSR, CSC, COO, ELL, HYB와 같은 다양한 형식을 지원한다. 행렬-벡터 곱은 물론이고, 특정 분해(ILU, IC)나 트라이앵글 풀이(Triangular Solve) 등을 빠르게 수행하도록 최적화되어 있다. 이런 라이브러리를 활용하면, 사용자는 단순히 CSR 배열을 넘기고 함수 하나만 호출해도 SpMV가 실행된다. 내부적으로는 GPU가 사용하는 최적화된 커널이 동작하거나, CPU 벡터화(AVX, AVX-512 등)를 활용하는 경로가 선택된다. 대규모 코드에서 희소 연산을 반복적으로 수행할 때, 반복법의 모든 단계는 결국 SpMV나 Sparse Triangular Solve 같은 연산을 여러 번 호출하게 되므로, 라이브러리 최적화 수준이 문제 해결 시간에 결정적 영향을 미친다. #### 도메인 분할과 멀티 피지컬 시뮬레이션 수치해석에서 희소 행렬이 큰 비중을 차지하는 분야 중 하나가 유한 요소법(Finite Element Method)이다. 물리적 영역(domain)을 여러 개의 소영역(Element)으로 분할해, 각 요소에서 형성되는 국소 행렬(또는 스텐실)을 조합하면 전역 행렬 $\mathbf{A}$가 완성된다. 이 과정에서 대부분의 요소가 서로 상호작용하는 범위는 국소적이기 때문에, 전역적으로 0이 아닌 항이 매우 제한된 곳에만 형성되는 희소 행렬이 만들어진다. 또한 열-유동, 열-구조, 전자기-열, 음향-구조 등 복합 물리(멀티 피지컬) 문제를 풀 때, 블록 구조의 희소 행렬이 탄생한다. 여러 물리 방정식을 한 번에 풀어야 하는 시스템은 블록 행렬 형태로 표현되며, 각 블록 내부는 서로 다른 변수 집합을 의미한다. 이런 상황에서는 블록 계수 간 연결을 잘 파악해 재정렬하거나 블록 전처리(블록 ILU, 블록 대각 사전조건자 등)를 적용하는 전략이 쓰인다. #### 희소 행렬과 그래프 이론 희소 행렬은 행과 열을 그래프의 정점(Vertex)이라 보고, 0이 아닌 항을 에지(Edge)로 해석하면 자연스럽게 그래프 구조가 된다. 예를 들어 대칭 행렬이라면 무방향 그래프로 해석할 수 있으며, 비대칭 행렬이라면 유향 그래프를 고려하게 된다. 이렇게 변환해 놓으면, fill-in 현상은 그래프에서 새로운 에지가 형성되는 일로 볼 수 있고, 재정렬은 그래프의 정점 순서를 다시 정하는 문제와 동일하다. Nested Dissection 기법은 그래프에서 분할 정복(Divide and Conquer) 개념을 적용하여, 작은 소구역으로 나눈 뒤 경계 정점에 대한 순서를 재구성하는 방식으로 fill-in을 줄이는 기법이다. Metis, Scotch, ParMetis, PT-Scotch 등은 그래프 파티셔닝 및 재정렬에 정평이 나 있으며, 이들은 희소 행렬의 대규모 연산에서 사실상 표준 툴로 자리잡았다. #### 응용 예시(Octave) 아래는 Octave에서 희소 행렬을 정의하고 반복법(CGM)으로 $\mathbf{A}\mathbf{x}=\mathbf{b}$를 푸는 간단한 예시다. Octave 역시 내부적으로 희소 행렬 구조를 효율적으로 관리하며, 기본 반복법 함수를 제공한다. ```octave % Octave code snippet % 희소 행렬 생성 n = 5; i = [1; 1; 2; 3; 4; 5]; j = [1; 2; 1; 3; 5; 4]; v = [10; 2; 3; 5; 1; -1]; A = sparse(i,j,v,n,n); % b 벡터 정의 b = [1; 2; 3; 4; 5]; % Conjugate Gradient로 Ax=b 풀이 % (단, A는 대칭 양의 정부호 행렬이 아니므로 CG가 불안정할 수 있음) x0 = zeros(n,1); [x, flag, relres, iter] = cgs(A, b, 1e-6, 100, [], [], x0); disp("Solution x = "), disp(x'); ``` 위 예시는 $5\times 5$ 희소 행렬을 Octave의 sparse 구조로 정의하고, $\mathbf{b}$를 만들어서 `cgs` 함수를 호출한다. 실제로는 CG 알고리즘은 대칭 양의 정부호 행렬에 적용하는 것이 보통이므로, 이 예시에선 꼭 CG가 안정적으로 수렴한다는 보장은 없다. 하지만 희소 연산 흐름을 이해하는 데는 도움이 된다. `A = sparse(i,j,v,n,n)` 호출 시, 내부적으로 COO에 가까운 초기 입력을 받아 메모리 내에서 압축 구조를 만든다. 이후 `cgs` 수행 과정에서 행렬-벡터 곱이 반복적으로 일어날 때, Octave는 희소 형식의 장점을 활용한다. #### 고차원 문제로의 확장 차원이 수만\~수십만을 넘어가는 고차원 문제에서도, 희소 행렬 기법은 사실상 유일한 해결책에 가깝다. Dense 행렬을 전부 저장하려면 $n^2$ 규모의 메모리가 필요하지만, 희소 구조를 이용하면 $nnz$가 $n$의 선형적 스케일로 늘어날 뿐이고, 효율적인 연산이 가능해진다. 물론 분해법을 쓸 경우 $nnz$가 분해 과정에서 급격히 늘어날 가능성이 있으므로, 주로 반복법을 먼저 고려하게 된다. 반복법에서도, 필요한 연산의 대부분은 $\mathbf{A}\mathbf{x}$ 형태의 곱이므로, 이 곱을 병렬화(또는 GPU화)하여 빠르게 돌리고, 사전조건자를 적절히 구성하면 상당히 효율적으로 대규모 시스템을 풀 수 있다. 대표적인 대규모 사례로는 해석학습(Recommender System), 그래프 랭킹, 기계학습의 대규모 데이터 행렬 분해, 물리 시뮬레이션(CFD, 구조해석, 전자기 해석 등), 금융(파생상품 모델링) 등이 있다. #### 중간 요약 및 향후 흐름 지금까지 희소 행렬의 대표적 저장 방식(CoO, CRS, CCS 등), fill-in 문제와 재정렬, GPU 병렬화, 반복법과 사전조건자, 라이브러리 활용, 그리고 구체적 코드 예시 등을 폭넓게 살펴보았다. 대규모 선형계 문제는 실제 산업과 학계에서 매우 빈번히 등장하므로, 희소 구조를 제대로 이해하고 구현할 줄 아는 능력은 곧 수치해석 전반에서 핵심 역량으로 직결된다. 뒤이어 이어질 내용에서는 희소 시스템에 특화된 선형대수 알고리즘(직접법, 반복법, 혼합형 전략 등)을 좀 더 구체적으로 다루고, 다양한 응용 사례를 살펴보겠다. #### 다중 수준(multilevel) 방법과 AMG(Algebraic Multigrid) 희소 행렬이 대표적으로 나타나는 상황 중 하나는 편미분방정식(PDE)을 유한 요소나 유한 차분 기법 등으로 이산화할 때다. 문제 크기가 커질수록 유한 요소 격자의 크기(노드 수)가 기하급수적으로 증가하고, 이에 따라 선형계의 차원도 매우 커진다. 이 경우 단순 반복법이나 고전적 사전조건자만으로는 수렴 속도가 크게 저하되기 쉽다. 이때 다중 수준(multilevel) 방법을 적용하면, 고차원 격자를 단계적으로 얇은(또는 조밀하지 않은) 격자로 투영(projection)하여 더 작은 차원의 문제로 접근한 뒤, 다시 상위 격자로 해를 보정(prolongation)하는 과정을 반복할 수 있다. 이렇게 하면 각 레벨(level)에서 해를 보완함으로써, 단일 레벨에서 고차원의 잔차를 해소하는 것보다 훨씬 빠른 수렴을 기대할 수 있다. 다중 수준 방법의 대표적인 구현 중 하나가 ‘Algebraic Multigrid(AMG)’다. AMG는 기하학적(Geometric) 정보를 명시적으로 사용하지 않고, 행렬 $\mathbf{A}$ 자체의 계수 행렬로부터 강약 연결(strong/weak connection) 정보를 자동으로 추출하여 ‘보조 격자(coarse grid)’ 구조를 구성한다. 따라서 임의 형태의 격자나 일반 그래프 형태 문제에도 적용 가능하다는 장점이 있다. AMG는 대규모 희소 선형계의 해법으로 널리 알려져 있으며, 다양한 라이브러리(Hypre, ML/Trilinos, BoomerAMG 등)에서 구현되어 있다. #### Additive Schwarz, Domain Decomposition, 그리고 병렬화 또 다른 유력한 접근법으로 Domain Decomposition(DD) 기법이 있다. 해석 영역을 여러 개의 부분영역(subdomain)으로 분할하고, 각 부분에서 부분계(subsystem)를 해결한 뒤, 경계 접합부에서 조건을 조정하면서 전체 문제를 풀어가는 방식이다. 이때 각 부분영역별 선형계가 상당히 독립적인 형태로 분리되므로, 병렬화에 유리하며, 메모리 및 연산을 분산할 수 있다. DD 기법은 필연적으로 희소 행렬 구조를 유지하며, 부분영역의 크기나 인접 관계, 그리고 경계 조건을 어떻게 처리하는지가 핵심 이슈다. Additive Schwarz 방법은 여러 부분영역을 동시에(병렬로) 풀어 놓고, 얻은 해를 합쳐나가는 전략이다. 보통 DD 기법과 함께 사전조건자로 활용하여, 반복법에서의 수렴을 가속하는 구성이 많다. 예를 들어, 전역 문제에서 얻은 잔차를 각 부분영역에서 국소적으로 해석하고, 그 결과를 합산(additive)해서 해를 갱신한다. 이는 전역 문제를 직접 푸는 것보다 훨씬 가벼운 연산으로 나누어 처리할 수 있고, 부분영역 내에서의 희소 구조(또는 블록 구조)를 활용해 효율적으로 접근할 수 있다. #### Krylov 하위공간 반복법과 고급 사전조건자 희소 행렬을 다루는 반복법은 기본적으로 Krylov 하위공간(Krylov Subspace)을 단계적으로 확장하며 해를 추적하는 방식이다. CG, BiCG, GMRES, BiCGSTAB, TFQMR, IDR 등 여러 알고리즘이 존재하며, 각각의 장단점이 다르다. 공통적으로 스펙트럼 분포가 광범위하거나 조건수가 큰 문제에서는 사전조건자(Preconditioner) 없이는 수렴이 매우 느려진다. 이를 해결하기 위해 ILU(불완전 LU), IC(불완전 Cholesky) 외에도, 다중 수준 사전조건자(AMG), 블록 대각 사전조건자, Domain Decomposition 계열 사전조건자 등을 적극 활용한다. 고급 사전조건자를 구성할 때도, 희소 행렬의 구조적 특성이 중요하며, 분해 과정에서의 fill-in을 어느 정도 허용할 것인지, 혹은 임계값(threshold)를 어디까지 허용할 것인지 등을 조정해야 한다. 지나치게 fill-in을 허용하면 사전조건자가 “과도하게 Dense”해져서 메모리/연산비가 커지고, fill-in을 너무 제한하면 사전조건 효과가 낮아져서 수렴이 늦어진다. 이 균형점을 찾아내는 것이 숙련된 실무자의 역할이 되기도 한다. #### 수치 안정성(Numerical Stability)과 복합 정밀도(Precision) 희소 연산이 대규모로 진행될 때는 부동소수점 정밀도(단정도, 배정도, 혹은 혼합 정밀도) 설정이 중요하다. 완전 배정도로 모든 연산을 수행하면 수치 오차는 상대적으로 줄어들지만, 메모리 사용량과 연산 비용이 증가한다. 반대로 단정도(half, single)로 모든 연산을 하면 빠르지만, 수렴 불안정이나 정밀도 손실 위험이 커진다. 최근에는 배정도와 단정도를 혼합하는 Mixed Precision 기법이 다수 연구되었다. 예컨대 전처리 과정(예: ILU, AMG)을 단정도로 수행하여 빠르게 근사한 사전조건자를 구한 뒤, 주 반복 단계(Krylov 반복)에서는 배정도를 사용해 최종 해를 refine하는 식이다. 이 기법은 GPU 연산에서 단정도 연산이 훨씬 빠른 점을 활용하며, 동시에 수치 안정성을 어느 정도 확보하는 절충안으로 주목을 받고 있다. #### 희소 행렬 문제에서의 에러 분석과 수렴 판정 희소 행렬에 대해 반복법을 사용할 때, 잔차 $r = \mathbf{b} - \mathbf{A}\mathbf{x}$의 크기와 함께 미지수 벡터 $\mathbf{x}$ 자체의 오차를 동시에 모니터링하는 것이 좋다. 어떤 문제에서는 $|r|$이 줄어도 실제 해 $\mathbf{x}$와의 상대 오차가 커질 수 있으며(특히 $\mathbf{A}$가 ill-conditioned인 경우), 반대로 $|r|$만으로 충분히 정확하게 수렴을 판정할 수도 있다. 실제 계산에서는 $|r|$가 지정된 공차(tolerance) 이하가 되면 수렴했다고 간주하는 방식이 일반적이다. 이때 흔히 사용하는 공차는 $10^{-6}$, $10^{-8}$, $10^{-12}$ 등 문제에 따라 달라진다. CFD나 구조해석과 같은 물리 시뮬레이션의 경우, 원하는 결과 물리량(변형량, 압력, 속도 등)이 특정 민감도 이하의 차이만 있어도 충분하다면 상대 오차 기준을 조금 완화할 수도 있다. 반면 고정밀 금융 계산이나 예민한 최적화 문제는 더 엄격한 공차를 요구한다. #### 동적(Adaptive) 전략 희소 행렬 문제를 반복법으로 풀다가, 특정 구간에서 수렴이 더뎌지면 사전조건자를 동적으로 바꾸거나(Adaptive preconditioning), 재정렬을 시도하거나, 더 높은 차수의 투영기법이나 다중 수준 기법으로 전환하는 전략도 가능하다. 이처럼 계산 도중에 선택 알고리즘을 전환하는 방식을 ‘동적(Adaptive)’ 접근이라 부른다. 복잡한 물리 시뮬레이션에서는, 비선형 해석 과정에서 야코비(Jacobian)이 바뀔 때마다 희소 구조가 달라질 수 있으므로, 분해나 재정렬을 몇 번이고 반복하는 경우도 많다. 계산 시간이 길어질 수 있으나, 대규모 문제에서 잘못된 사전조건자나 재정렬을 유지하는 것보다는, 중간에 업데이트하는 편이 전체 계산 시간을 줄이는 결과가 될 수 있다. #### 실무 최적화 팁 실제 프로젝트에서 희소 행렬을 다룰 때는 다음과 같은 점을 신중히 고려한다. * **배열 및 메모리 할당** 0이 아닌 항($nnz$)을 정확히 혹은 대략적으로 예측해서, 충분한 배열 크기를 확보한 뒤 연산을 시작한다. 동적 할당/해제가 빈번하면 성능이 급격히 떨어진다. * **재정렬 / 파티셔닝** 분해법을 사용하는 경우 특히 fill-in 억제를 위한 재정렬이 필수적이다. 또한 대규모 병렬/분산 환경에서는 영역 파티셔닝(도메인 분할)과 그래프 파티셔닝을 적절히 조합한다. * **반복법과 사전조건자 선택** 대상 행렬의 성질(대칭 양의 정부호, 비대칭, 구조적 대칭 등)에 따라 CG, GMRES, BiCG 계열 등을 선택하고, 필요한 사전조건자(ILU, IC, AMG, DD 등)를 맞춤 설정한다. * **수치적 유의사항** ill-conditioned 문제가 예상되면 더 높은 정밀도나, 반대로 혼합 정밀도 기법을 고려한다. 반복 중 잔차가 급격히 튀거나(NaN 발생 등) 발산하면 사전조건자 또는 재정렬 전략을 재검토해야 한다. * **라이브러리와 하드웨어 활용** 최신 CPU/GPU의 벡터화 기능(AVX-512, CUDA, HIP 등), 또는 Sparse BLAS 라이브러리를 적극적으로 이용한다. 분산 메모리 환경이라면 MPI, OpenMP, CUDA/OpenACC 등 멀티 모델을 결합할 수도 있다. #### 간단한 C++ 예: ILU 전처리와 GMRES 직접 예시로, 희소 행렬에 대해 ILU 전처리를 적용한 뒤 GMRES를 수행하는 구조를 보여주겠다. 다음은 C++의 Trilinos 라이브러리를 예시로 든 간소화된 개념 코드다. ```c++ #include "Epetra_ConfigDefs.h" #include "Epetra_SerialComm.h" #include "Epetra_CrsMatrix.h" #include "Epetra_Vector.h" #include "Teuchos_ParameterList.hpp" #include "Ifpack.h" #include "AztecOO.h" int main() { Epetra_SerialComm Comm; // 행렬 차원 int n = 5; Epetra_Map Map(n, 0, Comm); Epetra_CrsMatrix A(Copy, Map, 0); // 희소 행렬 삽입 // 예: A(0,0) = 10, A(0,1)=2, ... double val; int col; val = 10.0; col = 0; A.InsertGlobalValues(0, 1, &val, &col); val = 2.0; col = 1; A.InsertGlobalValues(0, 1, &val, &col); val = 3.0; col = 0; A.InsertGlobalValues(1, 1, &val, &col); val = 5.0; col = 2; A.InsertGlobalValues(2, 1, &val, &col); val = 1.0; col = 4; A.InsertGlobalValues(3, 1, &val, &col); val = -1.0; col = 3; A.InsertGlobalValues(4, 1, &val, &col); A.FillComplete(); // b, x 벡터 Epetra_Vector b(Map), x(Map); for (int i = 0; i < n; i++) { b[i] = (double) (i+1); x[i] = 0.0; } // ILU 전처리(Ifpack) Teuchos::ParameterList ifpackParams; ifpackParams.set("fact: level-of-fill", 0); Ifpack Factory; Ifpack_Preconditioner* Prec = Factory.Create("ILU", &A, 0); Prec->SetParameters(ifpackParams); Prec->Initialize(); Prec->Compute(); // GMRES(AztecOO) AztecOO solver; solver.SetUserMatrix(&A); solver.SetLHS(&x); solver.SetRHS(&b); solver.SetPrecOperator(Prec); solver.SetAztecOption(AZ_solver, AZ_gmres); solver.Iterate(100, 1.0e-8); // 결과 출력 x.Print(std::cout); delete Prec; return 0; } ``` 이 코드는 개념적 예시이며, 실제 빌드 환경이나 링크 설정, 라이브러리 버전에 따라 달라질 수 있다. 핵심은 다음과 같다. * `Epetra_CrsMatrix`로 희소 행렬을 생성하고, 0이 아닌 항을 여러 차례 `InsertGlobalValues`로 넣은 뒤 `FillComplete()`로 압축 구조를 확정한다. * `Ifpack` 팩토리를 통해 ILU 전처리기를 생성하고(“fact: level-of-fill”=0은 ILU(0)임), `Initialize()`와 `Compute()`를 호출하여 분해를 수행한다. * `AztecOO` 솔버를 이용해 GMRES 반복법을 수행한다. ILU 전처리기는 `solver.SetPrecOperator(Prec)`로 등록된다. * 반복이 끝나면 해 벡터 `x`를 출력한다. 실무에서는 여기서도 재정렬(인덱스 순서 변경), fill-in 레벨, 임계값 등 다양한 파라미터를 조절할 수 있으며, 행렬 구조나 문제 규모에 따라 Sparse BLAS를 GPU 버전으로 바꿔 쓸 수도 있다.