지난 글에서 우리는 LU 분해를 통해 역행렬, 선형 방정식 풀이를 구현하며 선형대수의 기반 기능을 갖추었습니다. 이제 한 단계 더 나아가 QR 분해(Orthogonal Decomposition), SVD(특이값 분해, Singular Value Decomposition), 고유값/고유벡터(Eigenvalue/Eigenvector) 문제 해결 등의 고급 연산에 도전하겠습니다.
이러한 연산은 머신러닝(주성분 분석, 차원 축소), 신호처리(잡음 제거, 이미지 압축), 물리/공학 시뮬레이션(시스템 동작 분석) 등 매우 광범위한 영역에서 핵심 도구로 사용됩니다. 이번 글에서는 기초 구현과 개념 위주로 접근하며, 최적화와 수치안정성을 최대한 고려하려고 합니다.
목표 연산 정리
- QR 분해(QR Decomposition): MxN 행렬 A를 정규직교(Orthogonal) 행렬 Q와 상삼각 행렬 R로 분해하는 과정. 직교화 기법(그람-슈미트, Householder 반사)을 통해 구현할 수 있습니다. 여기서는 수치안정성이 좋은 Householder 반사 기반 QR 분해를 간단히 구현해 봅시다.
- SVD(특이값 분해): A를 U Σ V^T 형태로 분해하는 것. 여기서는 완전한 구현 대신, 개념과 간단한 방법(예: Golub–Kahan SVD 알고리즘 스케치)을 통해 SVD 구조를 잡아봅니다. 실전에서는 LAPACK 바인딩이나 Crate를 활용할 수도 있습니다.
- 고유값/고유벡터 문제(Eigenvalue/Eigenvector): 대칭 행렬에 대해 파워 메서드(Power Method)를 사용한 근사적 고유값 계산 기법을 구현해봅니다. 이후 Rayleigh Quotient Iteration, QR 알고리즘 등을 통해 더 정확한 고유분해를 시도할 수 있습니다.
QR 분해 (Householder 반사 기반)
Householder 반사를 이용하면, 수치안정성 높은 QR 분해를 구현할 수 있습니다. 기본 아이디어는 A를 순차적으로 Householder 행렬 H_i로 좌측에서 곱해, 상삼각 형태 R을 만들어내고, Q = H₁H₂...H_k 를 얻는 것입니다.
src/decomp.rs 파일을 만들어 QR 관련 코드를 작성해봅시다.
// src/decomp.rs
use crate::matrix::Matrix;
/// QR 분해 결과
pub struct QRDecomposition {
pub q: Matrix<f64>,
pub r: Matrix<f64>,
}
impl QRDecomposition {
pub fn decompose(a: &Matrix<f64>) -> Self {
let m = a.rows();
let n = a.cols();
let mut mat = a.data.clone();
// Householder 반사를 순차 적용
let min_dim = m.min(n);
let mut reflectors = Vec::new();
for k in 0..min_dim {
let col_start = k*n + k;
let col_len = m - k;
// 현재 열에서 벡터 x 추출
let mut x = vec![0.0; col_len];
for i in 0..col_len {
x[i] = mat[(k+i)*n + k];
}
// Householder 벡터 v 계산
let alpha = x[0].copysign(1.0) * x.iter().map(|v| v*v).sum::<f64>().sqrt();
x[0] -= alpha;
let norm = x.iter().map(|v| v*v).sum::<f64>().sqrt();
for v in &mut x { *v /= norm; }
// v를 이용해 H = I - 2 vv^T 형성
// mat에 직접 H를 적용해 R 형성
for j in k..n {
let mut dot = 0.0;
for i in 0..col_len {
dot += x[i]*mat[(k+i)*n + j];
}
dot *= 2.0;
for i in 0..col_len {
mat[(k+i)*n + j] -= dot * x[i];
}
}
reflectors.push((k, x));
}
// R은 mat의 상삼각 부분
let mut r_data = vec![0.0; m*n];
for i in 0..m {
for j in i..n {
if i <= j { r_data[i*n+j] = mat[i*n+j]; }
}
}
let r = Matrix::new(m, n, r_data);
// Q 계산: reflectors를 역순으로 적용
let mut q_data = vec![0.0; m*m];
for i in 0..m {
q_data[i*m+i] = 1.0; // 단위행렬
}
let mut q = Matrix::new(m,m,q_data);
for &(k, ref v) in reflectors.iter().rev() {
// Q <- Q * (I - 2 vv^T)
// vv^T * Q 부분을 효율적으로 계산
for j in 0..m {
let mut dot = 0.0;
for i in 0..v.len() {
dot += v[i]*q[(k+i,j)];
}
dot *= 2.0;
for i in 0..v.len() {
q[(k+i,j)] -= dot * v[i];
}
}
}
Self { q, r }
}
}
위 구현은 개념 시연용으로, 수치안정성이나 최적화 측면에서 단순화했습니다. 실제로는 더 복잡한 관리가 필요합니다.
예제
use mathlib::{matrix::Matrix, decomp::QRDecomposition};
fn main() {
let a = Matrix::new(3,2, vec![
1.0,2.0,
2.0,3.0,
3.0,4.0
]);
let qr = QRDecomposition::decompose(&a);
let q = &qr.q;
let r = &qr.r;
// A ≈ Q R 검증
let qr_a = q.mat_mat_mul(r);
for i in 0..a.rows() {
for j in 0..a.cols() {
assert!((qr_a[(i,j)] - a[(i,j)]).abs() < 1e-10);
}
}
}
SVD(특이값 분해)
SVD 구현은 매우 복잡합니다. 여기서는 골루브-레인절(Golub-Reinsch) 알고리즘 등 고급 알고리즘을 전부 구현하기는 어렵고, 개념적 스케치만 제시합니다. 실전에서는 ndarray-linalg 같은 크레이트나 LAPACK 바인딩을 사용해 SVD를 얻는 경우가 많습니다.
아래 코드는 완전한 구현이 아닌, SVD 개념 전달용 의사코드에 가깝습니다.
// src/decomp.rs (추가)
pub struct SVD {
pub u: Matrix<f64>,
pub s: Vector<f64>,
pub vt: Matrix<f64>,
}
impl SVD {
pub fn decompose(_a: &Matrix<f64>) -> Result<Self,String> {
// 여기서는 "not implemented"를 반환하거나,
// 실제 구현 대신 외부 라이브러리 호출 예시를 남김
Err("SVD 구현은 복잡하며, LAPACK 바인딩 활용을 권장".to_string())
}
}
실제 SVD는 Householder 반사를 통해 A^T A를 대각화한 뒤 특이값을 얻거나, bidiagonalization 후 QR 알고리즘을 적용하는 등의 과정이 필요합니다. 이 시리즈에서는 시간 관계상 스킵하거나 단순히 외부 라이브러리 의존을 권장합니다.
고유값/고유벡터(Eigenvalue) 계산 (파워 메서드)
파워 메서드는 대칭 행렬 A에 대해 최대 고유값 및 대응 고유벡터를 근사하는 단순한 방법입니다.
// src/decomp.rs (추가)
use crate::vector::Vector;
impl Matrix<f64> {
pub fn power_method(&self, max_iter: usize, tol: f64) -> (f64, Vector<f64>) {
let n = self.rows();
assert_eq!(n, self.cols());
let mut x = Vector::new(vec![1.0; n]);
let mut lambda_old = 0.0;
for _ in 0..max_iter {
let y = self.mat_vec_mul(&x);
let lambda = y.data.iter().map(|v| v.abs()).fold(0./0., f64::max); // max norm
for i in 0..n {
x[i] = y[i]/lambda;
}
if (lambda - lambda_old).abs() < tol {
return (lambda, x);
}
lambda_old = lambda;
}
(lambda_old, x)
}
}
파워 메서드는 최대 절대값 고유값을 근사할 때 사용 가능하며, 보다 정교한 고유값 알고리즘(QR 알고리즘, Jacobi 알고리즘 등)은 후속 확장을 통해 구현할 수 있습니다.
앞으로의 학습 방향
3편에서는 QR 분해를 통해 정규직교 기저를 얻고 상삼각 형태로 행렬을 분해하는 과정을 살펴보았고, SVD와 고유값 문제를 개략적으로나마 다루어보았습니다. 다음 글에서는 다음과 같은 주제를 고려해볼 수 있습니다.
- 복소수 지원: Complex 타입을 지원하여 FFT나 고유값 문제에 더 광범위하게 대응
- 최적화 문제 해결용 라이브러리: 경사하강법, 뉴턴법, L-BFGS 등 최적화 알고리즘 구현
- 심볼릭 연산이나 자동미분(AD) 기반 기능 추가
- GPU 연동(CUDA, OpenCL, wgpu): 대규모 계산을 GPU 가속으로 처리
결론
이번 글에서 QR 분해, SVD 개념, 고유값 문제 등 보다 고급 연산 방향을 살펴보았습니다. 실제 구현은 매우 복잡하고 최적화가 필요하지만, 러스트의 타입 안전성, 모듈성, 성능 잠재력은 이런 고급 연산에 잘 어울립니다. 추후 외부 라이브러리나 BLAS/LAPACK 바인딩을 조합하면, 안정적이며 고성능인 수학 계산 라이브러리를 완성할 수 있을 것입니다.
유용한 링크와 리소스
- QR 분해 개념: Golub & Van Loan "Matrix Computations" 참고
- SVD 개념: Trefethen & Bau "Numerical Linear Algebra"
- Eigenvalue 문제: Wikipedia "Power Method", "QR algorithm"
- Rust Linear Algebra Crates: nalgebra, ndarray-linalg (LAPACK 바인딩)