이전 글에서는 벡터와 행렬의 기본 구조 및 연산(스칼라 곱, 벡터 덧셈, 행렬-벡터 곱, 행렬-행렬 곱)을 구현하며 러스트 기반 수학/과학 계산 라이브러리의 토대를 마련했습니다. 이번 글에서는 고급 선형대수 연산으로 한 걸음 더 나아가 LU 분해, 역행렬(Inverse), 그리고 선형 방정식 시스템(Ax = b) 해결 등을 구현해보겠습니다.
이러한 연산은 머신러닝, 최적화, 시뮬레이션, 통계 분석 등 다양한 분야에서 핵심적인 역할을 하며, BLAS/LAPACK 같은 고성능 라이브러리가 제공하는 기능을 러스트에서 직접 구현하고 최적화하는 과정은 러스트 기반 수학 라이브러리를 한 단계 성장시킵니다.
목표 연산 정리
- LU 분해(LU Decomposition): 정방 행렬 A를 하삼각 행렬 L과 상삼각 행렬 U로 분해하는 과정. 이를 통해 역행렬 계산, 선형 방정식 풀이 등 다른 문제를 효율적으로 해결할 수 있습니다.
- 역행렬(Inverse) 계산: LU 분해를 이용하여 A의 역행렬 A⁻¹을 구해봅니다. 단, A가 비가역일 경우(치환 행렬 사용 시 pivot이 0)에 대한 예외 처리도 중요합니다.
- Ax = b 선형 방정식 풀이: LU 분해를 통해 Ax = b 형태의 연립 방정식을 효율적으로 풉니다.
우리는 주로 f64 타입 행렬에 초점을 맞추겠지만, 향후 제네릭으로 확장할 수 있는 구조를 염두에 둘 수 있습니다.
LU 분해 구현
LU 분해에서는 부분 피봇팅(Partial Pivoting) 기법을 사용해 수치안정성을 확보하는 것이 일반적입니다. 여기서는 단순화를 위해 부분 피봇팅을 적용하겠습니다. 피봇팅을 반영하기 위해 **치환 벡터(Permutation Vector)**를 함께 관리합니다.
src/linalg.rs라는 모듈을 만들어 선형대수 관련 함수들을 모아봅시다.
// src/linalg.rs
use crate::matrix::Matrix;
/// LU 분해 결과를 담는 구조체
pub struct LUDecomposition {
pub lu: Matrix<f64>, // LU를 합친 행렬 (L과 U를 하나의 행렬에 저장)
pub pivot: Vec<usize>, // 피봇팅 정보
pub det_sign: f64, // 행 교환 횟수에 따른 부호 (-1 또는 1)
}
impl LUDecomposition {
/// 행렬 A를 LU분해
pub fn decompose(a: &Matrix<f64>) -> Result<Self, String> {
let n = a.rows();
assert_eq!(n, a.cols(), "정방 행렬만 LU분해 가능");
let mut lu_data = a.data.clone();
let mut pivot = (0..n).collect::<Vec<_>>();
let mut det_sign = 1.0;
for k in 0..n {
// 피봇 선택: 절댓값이 가장 큰 행 선택
let mut max_row = k;
let mut max_val = lu_data[k*a.cols()+k].abs();
for i in (k+1)..n {
let val = lu_data[i*a.cols()+k].abs();
if val > max_val {
max_val = val;
max_row = i;
}
}
if max_val < 1e-14 {
return Err("행렬이 특이(singular)하거나 수치적으로 불안정".to_string());
}
if max_row != k {
// 행 교환
for c in 0..n {
lu_data.swap(k*n+c, max_row*n+c);
}
pivot.swap(k, max_row);
det_sign = -det_sign;
}
// LU 분해 진행
for i in (k+1)..n {
lu_data[i*n+k] /= lu_data[k*n+k];
let factor = lu_data[i*n+k];
for j in (k+1)..n {
lu_data[i*n+j] -= factor * lu_data[k*n+j];
}
}
}
Ok(LUDecomposition {
lu: Matrix::new(n, n, lu_data),
pivot,
det_sign,
})
}
/// 행렬식(determinant) 계산: U의 대각원소 곱 * det_sign
pub fn det(&self) -> f64 {
let n = self.lu.rows();
let mut d = self.det_sign;
for i in 0..n {
d *= self.lu[(i,i)];
}
d
}
}
여기서는 LU 분해 로직을 단순히 구현했지만, 실제로는 수치안정성을 고려한 다양한 최적화 기법을 적용할 수 있습니다.
역행렬 계산
LU 분해를 통해 A를 LU로 쪼갰다면, 역행렬 A⁻¹을 구하는 것은 각 표준 기저 벡터 e_i에 대해 Ax = e_i를 풀어 x를 구한 뒤, 이 x를 A⁻¹의 열로 구성하는 방식으로 할 수 있습니다.
// src/linalg.rs (추가)
use crate::vector::Vector;
impl LUDecomposition {
/// LU분해 결과를 이용해 Ax = b 형태의 방정식 해를 구한다.
fn solve(&self, b: &Vector<f64>) -> Vector<f64> {
let n = self.lu.rows();
assert_eq!(n, b.len());
let mut x = vec![0.0; n];
// 피봇팅 적용
for i in 0..n {
x[i] = b[self.pivot[i]];
}
let lu_data = &self.lu.data;
// 전방대치(forward substitution)
for i in 0..n {
let sum = x[i];
let row_start = i*n;
let val = (i+1..n).fold(sum, |acc, j| acc - lu_data[j+row_start]*x[j]);
x[i] = val;
}
// 후방대치(back substitution)
for i in (0..n).rev() {
let row_start = i*n;
let sum = x[i];
x[i] = sum / lu_data[i*n+i];
for j in 0..i {
x[j] -= lu_data[j+row_start] * x[i];
}
}
Vector::new(x)
}
/// 역행렬 계산
pub fn inverse(&self) -> Result<Matrix<f64>, String> {
let n = self.lu.rows();
// det이 0이면 역행렬 없음
if self.det().abs() < 1e-14 {
return Err("비가역 행렬".to_string());
}
let mut inv_data = vec![0.0; n*n];
// e_i 벡터 생성: i번째 열에 해당
for i in 0..n {
let mut e = vec![0.0; n];
e[i] = 1.0;
let x = self.solve(&Vector::new(e));
for r in 0..n {
inv_data[r*n+i] = x[r];
}
}
Ok(Matrix::new(n, n, inv_data))
}
}
선형 방정식 Ax = b 풀이
이미 solve 메서드가 구현되었으니, LU 분해만 하면 Ax = b 형태를 간단히 풀 수 있습니다.
// 사용 예제
use mathlib::{matrix::Matrix, linalg::LUDecomposition, vector::Vector};
fn main() {
let a = Matrix::new(3,3, vec![
2.0,1.0,1.0,
1.0,3.0,2.0,
1.0,2.0,3.0
]);
let b = Vector::new(vec![1.0,2.0,3.0]);
let lu = LUDecomposition::decompose(&a).unwrap();
let x = lu.solve(&b);
// Ax ≈ b 가 성립하는지 확인
let ax = a.mat_vec_mul(&x);
for i in 0..b.len() {
assert!((ax[i] - b[i]).abs() < 1e-12);
}
// 역행렬 확인
let inv_a = lu.inverse().unwrap();
let a_inv_a = a.mat_mat_mul(&inv_a);
// a_inv_a는 단위행렬에 가까워야 함
for i in 0..3 {
for j in 0..3 {
let expected = if i == j {1.0} else {0.0};
assert!((a_inv_a[(i,j)] - expected).abs() < 1e-12);
}
}
}
성능 개선 포인트
- 메모리 접근 패턴 최적화: LU 분해, 역행렬 계산 등에서는 캐시 친화적 접근 패턴 중요
- SIMD 활용: 루프 내 벡터 연산을 std::arch나 packed_simd등을 활용해 SIMD화 가능
- 멀티스레딩: 큰 행렬 연산에서 rayon 등을 사용해 병렬 연산
- 디버깅 도구 활용: 벤치마크(cargo bench), 프로파일링(perf, callgrind)을 통해 병목 식별 후 최적화
앞으로의 학습 방향
2편에서는 LU 분해와 이를 활용한 역행렬, 선형 방정식 풀이를 구현해보았습니다. 다음 글에서는 QR 분해, SVD(특이값 분해), 고유값/고유벡터(Eigenvalue/Eigenvector) 문제 등 더 발전된 선형대수 기법을 다루어볼 예정입니다. 또한, 복소수 지원, 심볼릭 연산, 최적화 문제 해결용 함수 구현 등으로 라이브러리를 확장할 수 있습니다.
이러한 과정을 통해 차츰 러스트로 구현한 수학 라이브러리가 고성능, 고품질, 유지보수 가능한 수준으로 성장하게 될 것입니다.
결론
이번 글에서 우리는 러스트를 사용해 LU 분해, 역행렬 계산, 선형 방정식 풀이라는 핵심 선형대수 기능을 구현했습니다. 이로써 전편에서 구축한 벡터/행렬 기반 위에 고급 연산 기능을 쌓아올리며, 점진적으로 로우레벨 BLAS/LAPACK와 유사한 기능을 러스트 생태계 안에 구현하고 있음을 확인할 수 있습니다.
유용한 링크와 리소스
- Numerical Recipes 책(개념 참고): LU 분해, 역행렬 계산 알고리즘에 대한 수학적 배경
- LAPACK 문서: http://www.netlib.org/lapack/
- Rust Numeric 라이브러리 예시: nalgebra, ndarray, linfa 등 참조