러스트로 만드는 수학 및 과학 계산 라이브러리 시리즈 - 1편: 벡터와 행렬 연산 기초

이번 시리즈에서는 러스트(Rust)를 활용해 수학 및 과학 계산 라이브러리를 만드는 과정을 단계별로 살펴보겠습니다. 데이터 사이언스, 시뮬레이션, 머신러닝 등 다양한 분야에서 필수적인 벡터 연산, 행렬 연산, 선형대수, 확률/통계 연산, 고성능 FFT(Fast Fourier Transform)까지 확장 가능한 라이브러리를 직접 구현해볼 예정입니다.

이 시리즈는 다음과 같은 특징을 갖습니다.

• 안전한 메모리 관리: 러스트의 소유권과 빌림 규칙을 통해 메모리 안전성을 확보하면서 수학 연산을 처리합니다.
• 고성능 지향: 과학 계산에서 성능은 중요합니다. SIMD 활용, 메모리 레이아웃 최적화, 멀티스레딩 등을 통해 C++ BLAS/LAPACK 수준의 성능에 도전합니다.
• 추상화와 제네릭 활용: 러스트의 제네릭과 트레이트를 활용하여 실수, 복소수, 부동소수점 정밀도(float, double), 정수 타입 등에 유연하게 대응할 수 있는 라이브러리 구조를 모색합니다.
• 테스트 주도 개발: 각 단계별로 유닛 테스트와 벤치마크를 작성해 신뢰성 있는 수학 연산 코드를 완성합니다.

이번 1편에서는 벡터(Vector)와 행렬(Matrix)의 기본적인 자료구조 및 연산을 정의하는 것으로 시작하겠습니다. 이를 통해 향후 선형대수 연산, 고급 알고리즘 구현을 위한 토대를 마련합니다.

프로젝트 초기화

새로운 라이브러리를 위한 Cargo 프로젝트를 만듭니다.

cargo new mathlib --lib
cd mathlib

Cargo.toml:

[package]
name = "mathlib"
version = "0.1.0"
edition = "2021"

[dev-dependencies]
criterion = "0.4" # 벤치마크용

우리는 벤치마크를 위해 criterion 크레이트를 사용할 예정입니다. 아직 벤치마크 코드는 없지만 미리 설치해두겠습니다.

벡터 구조 정의

src/vector.rs 파일을 만들어, 벡터를 표현하는 타입을 정의해봅시다. 여기서는 일차원 배열 기반 벡터를 상정합니다.

// src/vector.rs
#[derive(Debug, Clone, PartialEq)]
pub struct Vector<T> {
    data: Vec<T>,
}

impl<T> Vector<T> {
    pub fn new(data: Vec<T>) -> Self {
        Self { data }
    }

    pub fn len(&self) -> usize {
        self.data.len()
    }
}

impl<T> std::ops::Index<usize> for Vector<T> {
    type Output = T;
    fn index(&self, idx: usize) -> &Self::Output {
        &self.data[idx]
    }
}

impl<T> std::ops::IndexMut<usize> for Vector<T> {
    fn index_mut(&mut self, idx: usize) -> &mut Self::Output {
        &mut self.data[idx]
    }
}

이렇게 하면 Vector<T>는 내부적으로 Vec<T>를 보유하며, 인덱싱으로 요소 접근이 가능합니다.

예제

use mathlib::vector::Vector;

fn main() {
    let v = Vector::new(vec![1.0, 2.0, 3.0]);
    assert_eq!(v.len(), 3);
    assert_eq!(v[0], 1.0);
}

벡터 연산: 스칼라 곱, 벡터 합

수학 라이브러리에서 기본적인 벡터 연산으로 스칼라 곱, 벡터 덧셈 등을 구현해봅시다. Vector<f64> 형태를 가정하고, 추후 제네릭으로 확장할 수 있도록 고려하겠습니다.

// src/vector.rs (추가)
impl Vector<f64> {
    pub fn scalar_mul(&mut self, scalar: f64) {
        for x in &mut self.data {
            *x *= scalar;
        }
    }

    pub fn add(&mut self, other: &Vector<f64>) {
        assert_eq!(self.len(), other.len());
        for (a, b) in self.data.iter_mut().zip(other.data.iter()) {
            *a += b;
        }
    }

    pub fn dot(&self, other: &Vector<f64>) -> f64 {
        assert_eq!(self.len(), other.len());
        self.data.iter().zip(other.data.iter())
            .map(|(a, b)| a * b)
            .sum()
    }
}

예제

let mut v1 = Vector::new(vec![1.0, 2.0, 3.0]);
let v2 = Vector::new(vec![4.0, 5.0, 6.0]);

v1.scalar_mul(2.0); // v1: [2.0, 4.0, 6.0]
v1.add(&v2);        // v1: [6.0, 9.0, 12.0]
let d = v1.dot(&v2); // dot product: 6*4 + 9*5 + 12*6 = 24 + 45 + 72 = 141

행렬 구조 정의

이제 2차원 배열로 행렬을 표현해봅시다. **행렬(Matrix)**은 보통 Row-major order로 저장하며, (rows, cols) 형태를 기억해야 합니다.

// src/matrix.rs
#[derive(Debug, Clone, PartialEq)]
pub struct Matrix<T> {
    rows: usize,
    cols: usize,
    data: Vec<T>,
}

impl<T> Matrix<T> {
    pub fn new(rows: usize, cols: usize, data: Vec<T>) -> Self {
        assert_eq!(data.len(), rows * cols);
        Self { rows, cols, data }
    }

    pub fn rows(&self) -> usize { self.rows }
    pub fn cols(&self) -> usize { self.cols }

    fn index_1d(&self, r: usize, c: usize) -> usize {
        r * self.cols + c
    }
}

impl<T> std::ops::Index<(usize, usize)> for Matrix<T> {
    type Output = T;
    fn index(&self, idx: (usize, usize)) -> &Self::Output {
        &self.data[self.index_1d(idx.0, idx.1)]
    }
}

impl<T> std::ops::IndexMut<(usize, usize)> for Matrix<T> {
    fn index_mut(&mut self, idx: (usize, usize)) -> &mut Self::Output {
        let pos = self.index_1d(idx.0, idx.1);
        &mut self.data[pos]
    }
}

이로써 Matrix<T>를 (row, col) 튜플 인덱스로 접근할 수 있습니다.

예제

use mathlib::matrix::Matrix;

fn main() {
    let m = Matrix::new(2, 2, vec![1.0, 2.0, 3.0, 4.0]);
    assert_eq!(m[(0, 0)], 1.0);
    assert_eq!(m[(1, 1)], 4.0);
}

행렬 연산: 행렬-벡터 곱, 행렬-행렬 곱

행렬 연산은 수학 라이브러리의 핵심입니다. 먼저 행렬-벡터 곱을 구현해봅시다 (f64 전용).

// src/matrix.rs (추가)
use crate::vector::Vector;

impl Matrix<f64> {
    pub fn mat_vec_mul(&self, v: &Vector<f64>) -> Vector<f64> {
        assert_eq!(self.cols, v.len());
        let mut result = vec![0.0; self.rows];
        for r in 0..self.rows {
            let mut sum = 0.0;
            for c in 0..self.cols {
                sum += self[(r, c)] * v[c];
            }
            result[r] = sum;
        }
        Vector::new(result)
    }

    pub fn mat_mat_mul(&self, other: &Matrix<f64>) -> Matrix<f64> {
        assert_eq!(self.cols, other.rows);
        let mut result = vec![0.0; self.rows * other.cols];
        for r in 0..self.rows {
            for c in 0..other.cols {
                let mut sum = 0.0;
                for k in 0..self.cols {
                    sum += self[(r, k)] * other[(k, c)];
                }
                result[r * other.cols + c] = sum;
            }
        }
        Matrix::new(self.rows, other.cols, result)
    }
}

예제

let m = Matrix::new(2, 3, vec![1.0, 2.0, 3.0,
                               4.0, 5.0, 6.0]);
let v = Vector::new(vec![1.0, 0.0, -1.0]); 
// m * v = [1*1 + 2*0 + 3*(-1), 4*1 + 5*0 + 6*(-1)] 
//       = [1 - 3, 4 - 6] = [-2.0, -2.0]

let mv = m.mat_vec_mul(&v);
assert_eq!(mv, Vector::new(vec![-2.0, -2.0]));

let m2 = Matrix::new(3, 2, vec![1.0,4.0,
                                2.0,5.0,
                                3.0,6.0]);
// m * m2 = 2x3 * 3x2 = 2x2
// result = [[1*1+2*2+3*3, 1*4+2*5+3*6],
//           [4*1+5*2+6*3, 4*4+5*5+6*6]]
//        = [[1+4+9,4+10+18],[4+10+18,16+25+36]]
//        = [[14,32],[32,77]]

let mm = m.mat_mat_mul(&m2);
assert_eq!(mm, Matrix::new(2,2, vec![14.0,32.0,32.0,77.0]));

테스트와 벤치마크

tests/vector_test.rs:

use mathlib::vector::Vector;

#[test]
fn test_vector_operations() {
    let mut v1 = Vector::new(vec![1.0, 2.0, 3.0]);
    let v2 = Vector::new(vec![4.0, 5.0, 6.0]);

    v1.scalar_mul(2.0); 
    assert_eq!(v1, Vector::new(vec![2.0, 4.0, 6.0]));

    v1.add(&v2);
    assert_eq!(v1, Vector::new(vec![6.0, 9.0, 12.0]));

    let d = v1.dot(&v2);
    assert_eq!(d, 141.0);
}

cargo test로 테스트를 검증할 수 있습니다.

벤치마크는 benches/bench.rs 생성 후 cargo bench로 성능 측정.

// benches/bench.rs
use criterion::{black_box, Criterion, criterion_group, criterion_main};
use mathlib::{vector::Vector, matrix::Matrix};

fn bench_mat_mul(c: &mut Criterion) {
    let m = Matrix::new(100, 100, (0..10_000).map(|x| x as f64).collect());
    let m2 = Matrix::new(100, 100, (0..10_000).map(|x| (x*2) as f64).collect());

    c.bench_function("mat_mul_100x100", |b| {
        b.iter(|| black_box(&m).mat_mat_mul(black_box(&m2)))
    });
}

criterion_group!(benches, bench_mat_mul);
criterion_main!(benches);

cargo bench를 통해 행렬 곱 성능을 측정하고, 최적화 포인트를 탐색할 수 있습니다.

앞으로의 학습 방향

이번 1편에서는 벡터와 행렬의 기본 자료구조 및 기본 연산을 구현해보았습니다. 다음 편에서는 선형대수 연산(LU 분해, QR 분해), 특이값 분해(SVD), 고성능 연산을 위한 SIMD 활용, 복소수/다른 수 타입 지원 등을 다루며 라이브러리를 더욱 풍부하게 만들 예정입니다.

 

또한, C++ BLAS/LAPACK와 성능 비교나, 메모리 풀링, 스레드 병렬화, GPU 연동(CUDA, OpenCL), FFT 등 고급 주제도 후속 편에서 차근차근 시도해볼 계획입니다.

결론

러스트로 수학 및 과학 계산 라이브러리를 구현하기 위한 첫 걸음으로 벡터와 행렬 구조, 기본 연산을 살펴보았습니다. 안정적인 메모리 관리, 타입 안전성, 테스트와 벤치마크 지원 덕분에 처음부터 신뢰성 있는 코드 베이스를 구축할 수 있음을 확인할 수 있습니다.

앞으로 이 시리즈를 통해 점진적으로 복잡한 계산 기능을 구현하며, 러스트 기반 수학/과학 계산 생태계를 확장하는 경험을 쌓아봅시다.

유용한 링크와 리소스

• The Rust Book: https://doc.rust-lang.org/book/
• Criterion Benchmarks: https://bheisler.github.io/criterion.rs/
• Rust Numeric Ecosystem 프로젝트들: nalgebra, ndarray, linfa 등 오픈소스 라이브러리 참고