As part of my journey learning Rust and diving deeper into linear algebra, I've started working on implementing a linear algebra library for the 42 School curriculum. The project, called "matrix", aims to implement fundamental vector and matrix operations from scratch.

Building a Linear Algebra Library in Rust ​

Setting Up the Project Structure ​

One of the first things that struck me about Rust was its excellent built-in tooling. The project structure was easy to set up using Cargo, Rust's package manager. I organized the code into three main modules:

• scalar.rs: Defines the Scalar trait and its implementations
• vector.rs: Implements the Vector type and its operations
• matrix.rs: Implements the Matrix type and its operations

The modularity of Rust's module system made it natural to separate these concerns while maintaining clear interfaces between them.

The Scalar Trait ​

The Scalar trait defines the requirements for numeric types that can be used in vectors and matrices:

pub trait Scalar:
    Copy +                  // Values must be copyable
    Add<Output = Self> +    // a + b
    Sub<Output = Self> +    // a - b
    Mul<Output = Self> +    // a * b
    Div<Output = Self> +    // a / b
    AddAssign +             // a += b
    SubAssign +             // a -= b
    MulAssign +             // a *= b
    DivAssign +             // a /= b
    Debug +                 // println!("{:?}", a)
    Display                 // println!("{}", a)
{
    fn zero() -> Self;
    fn one() -> Self;
}

I love how Rust's trait system allows for such clear specification of requirements. The trait bounds ensure that any type implementing Scalar has all the necessary operations for linear algebra.

Vector and Matrix Types ​

For both Vector and Matrix types, I had to make a key design decision - should operations modify in place or return new instances? I chose to implement in-place modifications:

impl<K: Scalar> Vector<K> {
    pub fn add(&mut self, other: &Vector<K>) {
        if self.size() != other.size() {
            panic!("Addition requires vectors of the same size.");
        }
        for (i, val) in other.data.iter().enumerate() {
            self.data[i] += *val;
        }
    }
}

This decision led to using panic! for error handling rather than returning Result types. While this might change as the project evolves (Rust generally prefers explicit error handling), it simplified the initial implementation.

Documentation and Testing ​

Rust's documentation system is phenomenal. Using doc-comments and the cargo doc command, I can automatically generate comprehensive documentation:

/// A vector of scalar values.
///
/// # Example
/// ```
/// use matrix::Vector;
/// let v = Vector::from([1.0, 2.0, 3.0]);
/// ```
pub struct Vector<K: Scalar> {
    data: Vec<K>
}

I've also set up GitHub Actions to automatically build and deploy this documentation to GitHub Pages: https://mlrcbsousa.github.io/matrix/

Exercise 00 - Add, Subtract and Scale ​

The first exercise implements basic vector and matrix operations: addition, subtraction, and scaling. All operations maintain O(n) time and space complexity as required.

The implementations are straightforward but required careful consideration of edge cases and error handling. For example, vector addition checks that dimensions match before proceeding:

pub fn add(&mut self, other: &Vector<K>) {
    if self.size() != other.size() {
        panic!("Addition requires vectors of the same size.");
    }
    // ... implementation
}

This is just the beginning - there are 14 more exercises to implement! I'll continue documenting my progress as I work through implementing determinants, matrix multiplication, and more complex linear algebra operations.

Exercise 01 - Linear combination ​

While implementing the linear combination operation, I encountered an interesting optimization opportunity through Fused Multiply-Add (FMA) operations. The subject mentioned some x86 assembly instructions - vfmadd132ps, vfmadd213ps, and vfmadd231ps - which led me down a fascinating path of hardware-level optimizations.

Understanding Fused Multiply-Add (FMA) ​

A Fused Multiply-Add computes (a * b) + c in a single operation with a single rounding step. Compared to performing multiplication and addition separately, which requires two rounding steps, FMA provides better numerical precision and can be faster on hardware that supports it.

Finding the Right Implementation ​

The subject mentioned three x86 SIMD instructions: vfmadd132ps, vfmadd213ps, and vfmadd231ps. These are AVX instructions that perform FMA with different orderings of operands. At the lowest level, Rust provides access to these through the std::arch::x86_64 module's intrinsic fmaf32.

We could implement FMA using these low-level intrinsics:

unsafe {
    use std::arch::x86_64::{
        __m128, _mm_fmadd_ps, // SIMD intrinsics
        _mm_load_ps, _mm_store_ps
    };

    // Unsafe, requires CPU feature checking
    if is_x86_feature_detected!("fma") {
        // Maps to vfmadd213ps under the hood
        let result = _mm_fmadd_ps(a, b, c);
    }
}

However, Rust provides a much safer and more portable solution through the stabilized mul_add method on floating point types:

impl Scalar for f32 {
    fn fma(a: Self, b: Self, c: Self) -> Self {
        // Safe, portable FMA using std lib
        // Internally optimizes to appropriate FMA instruction
        a.mul_add(b, c)
    }
}

The mul_add method provides a safe interface that will:

1. Use hardware FMA instructions when available
2. Fall back to separate multiply-add when needed
3. Handle CPU feature detection automatically
4. Work across different architectures

This exemplifies Rust's philosophy of providing safe abstractions over low-level features while maintaining performance. The compiler can still optimize this to use the appropriate SIMD instructions where available, but we get the benefits of safe, portable code.

Why Not Use Vector Add and Scale? ​

A reasonable question arose - why not implement linear combination using the previously implemented vector addition and scaling operations? While it would work:

// Possible but less efficient implementation
let mut result = Vector::new(vec![K::zero(); size]);
for (vector, &coef) in u.iter().zip(coefs.iter()) {
    let mut temp = vector.clone();
    temp.scl(coef);
    result.add(&temp);
}

This approach has several disadvantages:

1. Space Complexity: Creating temporary vectors for scaling would make the space complexity O(nm) where n is vector length and m is number of vectors. Using FMA maintains O(n).

2. Time Complexity: Each iteration requires O(n) operations for scaling and O(n) for addition. While the total is still O(n), FMA reduces the constant factor.

3. Precision: Multiple rounding steps between operations can accumulate numerical errors. FMA's single rounding provides better numerical stability.

Implementation Results ​

The final implementation achieves both performance and precision goals:

pub fn linear_combination<K: Scalar>(u: &[Vector<K>], coefs: &[K]) -> Vector<K> {
    // ... input validation ...

    let mut result = vec![K::zero(); size];

    // Single pass using FMA
    for (vector, &coef) in u.iter().zip(coefs.iter()) {
        for (r, &el) in result.iter_mut().zip(vector.data.iter()) {
            *r = K::fma(coef, el, *r);
        }
    }

    Vector::new(result)
}

This implementation maintains:

• O(n) space complexity - only allocates result vector
• O(n) time complexity but with minimal operations
• Maximum precision through FMA
• Clean, safe code using Rust standard library features

The journey through implementing FMA taught me valuable lessons about numerical computing, hardware optimizations, and Rust's excellent support for both high-level safety and low-level performance.

Exercise 02 - Linear interpolation ​

Linear interpolation is a fundamental operation in graphics, animation, and scientific computing. While conceptually simple, implementing it generically in Rust raised some interesting design decisions.

Two Approaches to Lerp ​

There are two common formulas for linear interpolation:

1. lerp(u, v, t) = u + t(v - u)          // Direct formula
2. lerp(u, v, t) = (1-t)u + tv           // Weighted average

While mathematically equivalent, the weighted average form proved more suitable for generic implementation because:

1. Better numerical stability near t = 1
2. More intuitive geometric interpretation
3. Requires fewer operations when implemented with assignment operators

Generic Implementation with Assignment Operators ​

Rather than using standard operators (+, -, *), I chose to use assignment operators (+=, -=, *=) for a cleaner implementation:

pub fn lerp<T>(u: T, v: T, t: f32) -> T
where
    T: Clone + AddAssign + SubAssign + MulAssign<f32>
{
    let mut a = u.clone();
    let mut b = v.clone();
    a *= 1.0 - t;  // Scale first point by (1-t)
    b *= t;        // Scale second point by t
    a += b;        // Add together
    a              // Return result
}

Key design decisions:

1. Uses Clone to preserve input values
2. Assignment operators keep mutations clear
3. Works with any type implementing basic scaling and addition
4. Returns new value rather than modifying in place

Making Types Interpolatable ​

To make our existing Vector and Matrix types work with lerp, we added the necessary assignment operator traits:

impl<K: Scalar> MulAssign<f32> for Vector<K> {
    fn mul_assign(&mut self, rhs: f32) {
        self.scl(rhs);
    }
}

impl<K: Scalar> AddAssign for Vector<K> {
    fn add_assign(&mut self, rhs: Self) {
        self.add(&rhs);
    }
}

This pattern of implementing operator traits in terms of our existing methods:

1. Maintains consistent behavior
2. Reuses tested code
3. Keeps complexity O(n)
4. Makes our types more ergonomic to use

Understanding t Values ​

While t is typically used between 0 and 1 for interpolation, our implementation works for any t value:

• t = 0: Returns first input (u)
• 0 < t < 1: Returns weighted average
• t = 1: Returns second input (v)
• t < 0 or t > 1: Extrapolates beyond inputs

This flexibility can be useful for extrapolation and other mathematical operations.

The exercise demonstrated how Rust's trait system enables generic implementations while maintaining type safety and performance. The clean separation of traits also made it natural to extend our existing types to support new operations.

Exercise 03 - Dot product ​

The dot product (or inner product) is one of the most fundamental operations in linear algebra. It takes two vectors and produces a scalar that reveals important geometric relationships between them.

Implementation with FMA ​

Building on our previous work with FMA operations, we can implement the dot product efficiently and with good numerical precision:

impl<K: Scalar> Vector<K> {
    pub fn dot(&self, other: &Vector<K>) -> K {
        if self.size() != other.size() {
            panic!("Dot product requires vectors of the same size");
        }

        let mut result = K::zero();
        for (x, y) in self.data.iter().zip(other.data.iter()) {
            result = K::fma(*x, *y, result); // Use FMA for better precision
        }
        result
    }
}

Key features:

1. Uses iterator zip for clean element pairing
2. Leverages FMA for numerical stability
3. O(n) time complexity
4. O(1) space complexity (just the accumulator)

Geometric Interpretation ​

The dot product has beautiful geometric properties that make it essential in graphics, physics, and many other domains:

// Orthogonal vectors have dot product 0
let v1 = Vector::from([1.0, 0.0]);
let v2 = Vector::from([0.0, 1.0]);
assert_eq!(v1.dot(&v2), 0.0);

// Parallel vectors show magnitude relationship
let v1 = Vector::from([2.0, 0.0]);
let v2 = Vector::from([4.0, 0.0]);
assert_eq!(v1.dot(&v2), 8.0);

// Opposing vectors have negative dot product
let v1 = Vector::from([1.0, 0.0]);
let v2 = Vector::from([-1.0, 0.0]);
assert_eq!(v1.dot(&v2), -1.0);

The dot product reveals:

• Vector alignment (perpendicular = 0)
• Magnitude relationships (parallel vectors)
• Directional relationships (positive vs negative)
• Angle between vectors (through the cosine formula)

FMA Benefits in Dot Product ​

Using FMA for dot product calculation provides two key benefits:

1. Better numerical precision by avoiding intermediate rounding
2. Potential hardware acceleration on supported platforms

Consider calculating dot product for long vectors:

// Without FMA:
result = (x₁ * y₁) + (x₂ * y₂) + ... + (xₙ * yₙ)
// Each parenthesis represents a rounding point

// With FMA:
result = FMA(x₁, y₁, FMA(x₂, y₂, ...))
// Single rounding at each step

This implementation shows how a deep understanding of hardware capabilities (FMA) and mathematical properties (dot product geometry) can lead to both efficient and numerically stable code.

The exercise also demonstrates Rust's strengths:

• Generic implementation through traits
• Zero-cost abstractions for hardware features
• Clear error handling
• Expressive iterator combinators

Once again, Rust's type system helps ensure correctness while the resulting code remains performant and readable.

Exercise 04 - Norm ​

A norm is a function that assigns a positive length or size to a vector. Different norms capture different aspects of a vector's "magnitude" and are useful in different contexts.

Type Conversion Strategy ​

While implementing norms, I encountered an interesting challenge: our generic Vector<K> needed to compute norms as f32 values, but our K type could be any numeric type. Initially, I considered adding an abs method to our Scalar trait, but realized we could simplify by just converting to f32 first:

pub trait Scalar:
    // ... other trait bounds ...
    Into<f32>  // Single conversion trait handles all our needs
{
    fn zero() -> Self;
    fn one() -> Self;
    fn fma(a: Self, b: Self, c: Self) -> Self;
}

This elegant solution leverages Rust's standard library conversions between numeric types. Instead of requiring each Scalar type to implement its own abs, we simply convert to f32 and use its abs method.

Three Norms Implementation ​

I implemented three different norms, each with its own geometric interpretation:

impl<K: Scalar> Vector<K> {
    // Manhattan norm - sum of absolute values
    pub fn norm_1(&self) -> f32 {
        let mut sum = 0.0;
        for val in &self.data {
            sum += Into::<f32>::into(*val).abs();  // Convert then abs
        }
        sum
    }

    // Euclidean norm - square root of sum of squares
    pub fn norm(&self) -> f32 {
        let mut sum_sq = 0.0;
        for val in &self.data {
            let val_f32: f32 = (*val).into();
            sum_sq = f32::fma(val_f32, val_f32, sum_sq);
        }
        sum_sq.sqrt()
    }

    // Supremum norm - maximum absolute value
    pub fn norm_inf(&self) -> f32 {
        self.data
            .iter()
            .map(|x| Into::<f32>::into(*x).abs())
            .fold(0.0, f32::max)
    }
}

Each norm has specific properties and use cases:

1. Manhattan Norm (1-norm)

   • ∥v∥₁ = Σ|vᵢ|
   • Like measuring distance by city blocks
   • Useful in error analysis and optimization

2. Euclidean Norm (2-norm)

   • ∥v∥₂ = √(Σvᵢ²)
   • Our familiar notion of distance
   • Uses FMA for better numerical precision

3. Supremum Norm (∞-norm)

   • ∥v∥∞ = max(|vᵢ|)
   • Maximum magnitude in any dimension
   • Important in convergence analysis

Norm Relationships ​

An interesting property of these norms is their relationship:

let v = Vector::from([3.0, -4.0]);
assert!(v.norm_inf() <= v.norm());  // 4.0 <= 5.0
assert!(v.norm() <= v.norm_1());    // 5.0 <= 7.0

For any vector v: ∥v∥∞ ≤ ∥v∥₂ ≤ ∥v∥₁

This implementation showcases several Rust strengths:

• Generic programming with minimal trait bounds
• Clean conversion between numeric types
• Iterators for expressive, efficient code
• FMA optimization for numerical stability

The exercise also demonstrated the value of simplifying trait requirements - by relying on conversion to f32 before applying operations like abs, we reduced the requirements on our generic type parameter while maintaining all functionality.

The journey from a more complex trait definition to this simpler version mirrors a common pattern in Rust development: start with the most general solution, then pare it down to the minimal requirements needed to solve the problem at hand.

Exercise 05 - Cosine ​

Calculating the angle between vectors is a fundamental operation in graphics, physics, and machine learning. The cosine of this angle provides a normalized measure of how aligned two vectors are.

Mathematical Foundation ​

The cosine formula between vectors is elegant:

cos(θ) = (u·v) / (∥u∥∥v∥)

Where:

• u·v is the dot product
• ∥u∥ is the Euclidean norm (2-norm)
• θ is the angle between vectors

Implementation Strategy ​

Having already implemented dot product and norm, we can combine them:

pub fn angle_cos<K: Scalar>(u: &Vector<K>, v: &Vector<K>) -> f32 {
    if u.size() != v.size() {
        panic!("Vectors must have same size");
    }

    let dot = u.dot(v);
    let norm_u = u.norm();
    let norm_v = v.norm();

    if norm_u == 0.0 || norm_v == 0.0 {
        panic!("Zero vectors have undefined angle");
    }

    // Use FMA for final multiplication
    dot / f32::fma(norm_u, norm_v, 0.0)
}

The implementation:

1. Validates input vectors
2. Uses previous dot product and norm implementations
3. Maintains numerical precision with FMA
4. Returns angle cosine in [-1, 1]

Geometric Interpretation ​

The cosine value tells us about vector alignment:

• cos(θ) = 1: Vectors point in same direction
• cos(θ) = 0: Vectors are perpendicular
• cos(θ) = -1: Vectors point in opposite directions
• cos(θ) ∈ (0,1): Acute angle
• cos(θ) ∈ (-1,0): Obtuse angle

// Perpendicular vectors (90°)
let v1 = Vector::from([1.0, 0.0]);
let v2 = Vector::from([0.0, 1.0]);
assert!(angle_cos(&v1, &v2).abs() < f32::EPSILON); // ~0.0

// Parallel vectors (0°)
let v3 = Vector::from([2.0, 0.0]);
let v4 = Vector::from([4.0, 0.0]);
assert_eq!(angle_cos(&v3, &v4), 1.0);

// Opposite vectors (180°)
let v5 = Vector::from([1.0, 0.0]);
let v6 = Vector::from([-1.0, 0.0]);
assert_eq!(angle_cos(&v5, &v6), -1.0);

Applications ​

This operation is crucial in many fields:

• Graphics: Camera angles and lighting
• Machine Learning: Cosine similarity for text analysis
• Physics: Force components and work calculations
• Data Science: Feature correlation

The implementation showed how previously built functions (dot product, norm) can be composed to create more complex operations while maintaining good numerical properties through FMA usage.

Exercise 06 - Cross product ​

The cross product is a fundamental operation in 3D geometry and physics - producing a vector perpendicular to two input vectors. Unlike our previous operations that worked in any dimension, cross product is specifically defined for 3D vectors.

Mathematical Foundation ​

For vectors u = [u₁, u₂, u₃] and v = [v₁, v₂, v₃], the cross product is:

u × v = [u₂v₃ - u₃v₂, u₃v₁ - u₁v₃, u₁v₂ - u₂v₁]

This formula yields a vector that is:

1. Perpendicular to both input vectors
2. Has magnitude equal to area of parallelogram spanned by inputs
3. Direction follows the right-hand rule

Implementation Strategy ​

pub fn cross_product<K: Scalar>(u: &Vector<K>, v: &Vector<K>) -> Vector<K> {
    if u.size() != 3 || v.size() != 3 {
        panic!("Cross product requires 3D vectors");
    }

    let (u1, u2, u3) = (u.data[0], u.data[1], u.data[2]);
    let (v1, v2, v3) = (v.data[0], v.data[1], v.data[2]);

    // Use FMA for better precision
    let x = K::fma(u2, v3, -(u3 * v2));
    let y = K::fma(u3, v1, -(u1 * v3));
    let z = K::fma(u1, v2, -(u2 * v1));

    Vector::from([x, y, z])
}

Key implementation features:

1. Generic over scalar type
2. Dimension checking
3. FMA for numerical stability
4. Clean component-wise calculation

Important Properties ​

The cross product has several key properties:

1. Right Hand Rule: Cross products with standard basis vectors follow circular pattern:

let i = Vector::from([1.0, 0.0, 0.0]);
let j = Vector::from([0.0, 1.0, 0.0]);
let k = Vector::from([0.0, 0.0, 1.0]);

assert_eq!(cross_product(&i, &j).data, k.data); // i × j = k
assert_eq!(cross_product(&j, &k).data, i.data); // j × k = i
assert_eq!(cross_product(&k, &i).data, j.data); // k × i = j

2. Anticommutative: Swapping inputs negates result:

// u × v = -(v × u)
let u = Vector::from([2.0, 3.0, 4.0]);
let v = Vector::from([5.0, 6.0, 7.0]);
let uv = cross_product(&u, &v);
let vu = cross_product(&v, &u);
// uv = -vu

3. Orthogonality: Result is perpendicular to both inputs:

let w = cross_product(&u, &v);
assert!(u.dot(&w).abs() < f32::EPSILON);  // u ⊥ w
assert!(v.dot(&w).abs() < f32::EPSILON);  // v ⊥ w

Applications ​

Cross product is essential in:

• Physics: Torque, angular momentum
• Computer Graphics: Surface normals, camera orientation
• Robotics: Robot arm movements
• Navigation: Course calculations

The implementation demonstrates how a seemingly complex operation can be expressed clearly in Rust while maintaining both safety (through dimension checks) and performance (through FMA usage).

Exercise 07 - Linear map, Matrix multiplication ​

Matrix multiplication is a fundamental operation that represents the composition of linear transformations. This exercise required implementing both matrix-vector multiplication (applying a linear transformation to a vector) and matrix-matrix multiplication (composing two linear transformations).

Mathematical Foundation ​

Matrix multiplication has two key forms:

1. Matrix-Vector multiplication (Mx = b):

bᵢ = Σⱼ(mᵢⱼxⱼ)  for i = 1..n

2. Matrix-Matrix multiplication (AB = C):

cᵢⱼ = Σₖ(aᵢₖbₖⱼ)  for i = 1..n, j = 1..p

Implementation Strategy ​

The implementation uses nested loops with FMA operations for better numerical precision:

pub fn mul_vec(&self, vec: &Vector<K>) -> Vector<K> {
    if self.cols() != vec.size() {
        panic!("Matrix columns must match vector size");
    }

    let mut result = vec![K::zero(); self.rows()];
    for i in 0..self.rows() {
        for (j, v_j) in vec.data.iter().enumerate() {
            result[i] = K::fma(self.data[i][j], *v_j, result[i]);
        }
    }
    Vector::new(result)
}

Key features:

1. Dimension validation
2. Efficient FMA usage
3. Clear row-by-row computation
4. Generic over scalar type

Geometric Interpretation ​

Matrix multiplication represents composition of linear transformations:

1. Matrix-Vector: Applying transformation to a point

// Identity matrix preserves vector
let m = Matrix::from([[1.0, 0.0], [0.0, 1.0]]);
let v = Vector::from([4.0, 2.0]);
assert_eq!(m.mul_vec(&v).data, v.data);

// Scaling matrix multiplies components
let scale = Matrix::from([[2.0, 0.0], [0.0, 2.0]]);
// Results in vector [8.0, 4.0]

2. Matrix-Matrix: Combining transformations

// Composing two transformations
let rotate = Matrix::from([[0.0, -1.0], [1.0, 0.0]]); // 90° rotation
let scale = Matrix::from([[2.0, 0.0], [0.0, 2.0]]);   // 2x scaling
// rotate.mul_mat(&scale) applies scaling then rotation

Computational Complexity ​

The implementations maintain the required complexity bounds:

• Matrix-Vector: O(nm) time, O(nm) space
• Matrix-Matrix: O(nmp) time, O(nm + mp + np) space

Where:

• n = rows of first matrix
• m = columns of first matrix/rows of second
• p = columns of second matrix

Applications ​

Matrix multiplication is crucial in:

1. Computer Graphics

   • View transformations
   • Perspective projection
   • Object transformations

2. Machine Learning

   • Neural network layers
   • Feature transformations
   • Batch operations

3. Physics Simulations

   • Coordinate transformations
   • Rigid body dynamics
   • Quantum state evolution

The implementation balances clarity, numerical stability (through FMA), and performance considerations, while maintaining type safety through Rust's type system.

Exercise 08 - Trace ​

The trace of a matrix is one of the simplest yet most important matrix operations. It has connections to eigenvalues, matrix similarity, and appears in many physical equations.

Mathematical Definition ​

The trace of a square matrix is the sum of elements on its main diagonal:

tr(A) = Σᵢaᵢᵢ = a₁₁ + a₂₂ + ... + aₙₙ

Implementation ​

The implementation is straightforward but requires checking that the matrix is square:

pub fn trace(&self) -> K {
    if !self.is_square() {
        panic!("Trace only defined for square matrices");
    }

    let mut sum = K::zero();
    for i in 0..self.rows() {
        sum += self.data[i][i];
    }
    sum
}

Key features:

1. Square matrix validation
2. Linear time complexity O(n)
3. Constant space complexity
4. Generic over scalar type

Important Properties ​

The trace has several interesting properties:

1. Linearity:

   tr(A + B) = tr(A) + tr(B)
   tr(cA) = c·tr(A)

2. Similarity Invariance:

   tr(AB) = tr(BA)

3. Identity Matrix:

   rust

   let identity = Matrix::from([
       [1.0, 0.0],
       [0.0, 1.0]
   ]);
   assert_eq!(identity.trace(), 2.0); // Equals dimension

4. Relation to Eigenvalues:

   • The trace equals the sum of matrix eigenvalues
   • This makes it useful in characteristic equations

Applications ​

The trace appears in many areas:

1. Physics

   • Quantum Mechanics: Trace of density matrix = 1
   • Stress tensors in mechanics
   • Conservation laws

2. Linear Algebra

   • Matrix similarity tests
   • Characteristic polynomial calculations
   • Matrix decompositions

3. Statistics

   • Covariance matrices
   • Principal Component Analysis
   • Fisher Information Matrix

The implementation shows how a simple operation (summing diagonal elements) can have deep mathematical significance and wide-ranging applications.

Exercise 09 - Transpose ​

The transpose operation is a fundamental matrix transformation defined as:

For a matrix A ∈ K^(m×n), the transpose A^T ∈ K^(n×m) is given by: (A^T)ᵢⱼ = Aⱼᵢ for all i,j

This "flips" the matrix across its main diagonal, swapping rows and columns. The operation has elegant mathematical properties:

1. Involution: (A^T)^T = A
2. Linear: (A + B)^T = A^T + B^T and (cA)^T = cA^T
3. Reverses multiplication: (AB)^T = B^T A^T

The implementation is straightforward but must handle dimensions carefully:

pub fn transpose(&self) -> Matrix<K> {
    let mut result = vec![vec![K::zero(); self.rows()]; self.cols()];
    for i in 0..self.rows() {
        for j in 0..self.cols() {
            result[j][i] = self.data[i][j];
        }
    }
    Matrix::new(result)
}

This achieves the required O(nm) time complexity for assignments and O(nm) space complexity. The transpose has important applications in normal equations (A^T A x = A^T b), matrix decompositions, and optimization algorithms.

Unlike operations like determinant or matrix multiplication, transpose is numerically stable as it only involves reordering values without arithmetic operations.

Exercise 10 - Row Echelon Form ​

The Reduced Row Echelon Form (RREF) transforms a matrix into a unique canonical form through Gaussian elimination. This standardized form makes matrix properties like rank and null space immediately visible.

Mathematical Foundation ​

A matrix R is in RREF if:

1. All zero rows are at the bottom
2. Leading coefficient (pivot) of each non-zero row is 1
3. Each pivot is the only nonzero entry in its column
4. Each pivot is right of the pivots in rows above

Formally, if matrix A reduces to RREF R through elementary row operations: A ~ R where R ∈ R^(m×n) is uniquely determined by:

• R_{i,j} = 1 if j = p_i (pivot positions)
• R_{k,p_i} = 0 for all k ≠ i
• p_i < p_{i+1} for all valid i

Implementation Strategy ​

Our implementation uses Gauss-Jordan elimination with two key enhancements:

1. Partial Pivoting: For each column, we choose the largest absolute value as pivot. This improves numerical stability by reducing error propagation, though at cost less than full pivoting would require.

2. Numerical Tolerance: Instead of checking x == 0 directly, we use |x| < TOLERANCE. This constant (1e-10) is:

   • Smaller than f32::EPSILON (≈1.19e-7) for stricter zero detection
   • Small enough to maintain precision for most practical applications
   • Critical for cleaning up "nearly zero" entries after elimination

The algorithm proceeds through these steps:

1. Forward elimination to create pivot columns
2. Each pivot operation:
   • Find largest value in current column
   • Scale row to make pivot 1
   • Eliminate entries in column
3. Clean up near-zero values

This requires O(n³) time but produces a unique, well-defined form useful for:

• Solving linear systems
• Computing matrix rank
• Finding null spaces
• Detecting linear dependencies

The implementation balances numerical robustness (through pivoting and tolerance) with computational efficiency while maintaining clarity of the underlying mathematical operations.

Exercise 11 - Matrix Determinants ​

The determinant is a scalar value that captures essential properties of a square matrix, with applications ranging from solving linear equations to computer graphics.

Mathematical Foundation ​

Definition and Properties ​

For a square matrix A, the determinant (denoted det(A) or |A|) has several equivalent definitions:

1. 2x2 Matrix

   |a b| = ad - bc
   |c d|

2. Larger Matrices For an n×n matrix, using cofactor expansion along row i:

   det(A) = Σⱼ (-1)ⁱ⁺ʲ aᵢⱼ det(Mᵢⱼ)

   where Mᵢⱼ is the minor matrix (remove row i and column j)

Key Properties ​

1. Invertibility

   • det(A) ≠ 0 ⟺ A is invertible
   • det(A) = 0 ⟺ A is singular (linear transformation "collapses" space)

2. Multiplicativity

   det(AB) = det(A)·det(B)

3. Determinant of Special Matrices

   • Identity matrix: det(I) = 1
   • Triangular matrices: product of diagonal elements
   • Transpose: det(Aᵀ) = det(A)

Geometric Interpretation ​

The determinant tells us about how a linear transformation affects space:

1. |det(A)| = scaling factor of area/volume
2. sign(det(A)) = whether orientation is preserved (+1) or reversed (-1)
3. det(A) = 0 means the transformation reduces dimension

Implementation Strategy ​

Our implementation uses a recursive approach with dimension-specific optimizations:

2×2 Determinant ​

fn det2x2(&self) -> K {
    let a = self.data[0][0];
    let b = self.data[0][1];
    let c = self.data[1][0];
    let d = self.data[1][1];

    // Use FMA: ad - bc
    K::fma(a, d, -(b * c))
}

The 2×2 case uses the direct formula with FMA for better numerical precision.

3×3 Determinant ​

fn det3x3(&self) -> K {
    let mut det = K::zero();

    // Expand along first row using 2x2 determinants
    for j in 0..3 {
        let minor = self.get_minor(0, j);
        let sign = if j % 2 == 0 { K::one() } else { -K::one() };
        det = K::fma(self.data[0][j] * sign, minor.det2x2(), det);
    }

    det
}

For 3×3 matrices, we use cofactor expansion along the first row, computing each term using the previously implemented 2×2 determinant.

4×4 Determinant ​

fn det4x4(&self) -> K {
    let mut det = K::zero();

    // Expand along first row using 3x3 determinants
    for j in 0..4 {
        let minor = self.get_minor(0, j);
        let sign = if j % 2 == 0 { K::one() } else { -K::one() };
        det = K::fma(self.data[0][j] * sign, minor.det3x3(), det);
    }

    det
}

The 4×4 case follows the same pattern but uses 3×3 determinants for the minors.

Numerical Considerations ​

Several techniques ensure numerical stability:

1. Fused Multiply-Add (FMA)

   • More precise than separate multiply and add
   • Hardware-accelerated on modern CPUs
   rust

   // Instead of a*d - b*c
   K::fma(a, d, -(b * c))

2. Recursive Structure

   • Each dimension reuses smaller determinant calculations
   • Ensures consistent handling of numerical issues
   • Makes code more maintainable

3. Memory Management

   • Pre-allocates vectors for minors
   • Minimizes temporary allocations
   rust

   let mut minor_data = Vec::with_capacity(n-1);

Complexity Analysis ​

Our implementation maintains:

• Time: O(n³)

  • Each n×n determinant requires n cofactor calculations
  • Each cofactor requires an (n-1)×(n-1) determinant
  • This recurrence solves to O(n³)

• Space: O(n²)

  • Storing minor matrices during recursion
  • Each level requires (n-1)×(n-1) space
  • Maximum depth of n-2 for our implementation

Exercise 12 - Matrix Inverse ​

Error Handling with Style ​

One elegant aspect of our implementation is the use of Rust's type system for error handling. Instead of using simple string errors or booleans, we created a dedicated error type:

pub enum MatrixError {
    Singular,    // Matrix isn't invertible (det ≈ 0)
    NotSquare,   // Matrix dimensions don't match
}

// Usage is clear and type-safe:
impl<K: Scalar> Matrix<K> {
    pub fn inverse(&self) -> Result<Matrix<K>, MatrixError> {
        if !self.is_square() {
            return Err(MatrixError::NotSquare);
        }
        // ... rest of implementation
    }
}

// Client code is expressive:
let m = Matrix::from([[1.0, 2.0], [2.0, 4.0]]);  // Singular matrix
match m.inverse() {
    Ok(inv) => println!("Inverse found: {}", inv),
    Err(MatrixError::Singular) => println!("Matrix is not invertible!"),
    Err(MatrixError::NotSquare) => println!("Matrix must be square!"),
}

This approach provides type safety, clear error cases, and elegant pattern matching for error handling.

The Classical Adjugate Method ​

While there are several ways to compute a matrix inverse (Gaussian elimination, LU decomposition), we implemented the classical adjugate method for its clarity and direct connection to earlier exercises. The inverse A⁻¹ of a matrix A is given by:

A⁻¹ = adj(A)/det(A)

where adj(A) is the adjugate matrix (transpose of cofactor matrix).

Implementation Strategy ​

Our implementation follows three key steps:

1. Early Validation:

   rust

   if !self.is_square() {
       return Err(MatrixError::NotSquare);
   }
   let det = self.determinant();
   if Into::<f32>::into(det).abs() < Self::TOLERANCE {
       return Err(MatrixError::Singular);
   }

2. Cofactor Calculation:

   rust

   for i in 0..n {
       for j in 0..n {
           let minor = self.get_minor(i, j);
           let cofactor = if (i + j) % 2 == 0 {
               minor.determinant()
           } else {
               -minor.determinant()
           };

3. Adjugate Construction and Scaling:

   rust

   // Transpose while building (j,i instead of i,j)
   // and divide by determinant
   inverse_data[j][i] = cofactor / det;

Numerical Considerations ​

Several techniques ensure numerical stability:

1. Tolerance-Based Singularity Check:

   rust

   const TOLERANCE: f32 = 1e-10;

   This constant helps identify "effectively zero" determinants.

2. Error Handling:

   rust

   pub enum MatrixError {
       Singular,
       NotSquare,
   }

   Clear error types guide users when inversion fails.

Performance Analysis ​

The implementation achieves the required complexity bounds:

• Time: O(n³) - dominated by determinant calculations
• Space: O(n²) - storage for result matrix

Properties and Verification ​

Key inverse properties our implementation maintains:

1. Identity Property: AA⁻¹ = A⁻¹A = I

   rust

   let prod = a.mul_mat(&a_inv);
   // Should approximate identity matrix

2. Scaling Property: (kA)⁻¹ = (1/k)A⁻¹

   rust

   let m = Matrix::from([[2.0, 0.0], [0.0, 2.0]]);
   let inv = m.inverse(); // Results in 0.5 * Identity

3. Error Cases:

   • Singular matrices (det ≈ 0) return Err(Singular)
   • Non-square matrices return Err(NotSquare)

This implementation combines clean code organization, proper error handling, and numerical stability while maintaining the required algorithmic efficiency.

Exercise 13 - Rank ​

The rank of a matrix is a fundamental concept in linear algebra. It reveals the maximum number of linearly independent rows or columns in a matrix, playing a crucial role in understanding its structure and properties.

Mathematical Foundation ​

The rank of a matrix is determined by reducing it to its row-echelon form (RREF) and counting the number of non-zero rows. Intuitively, the rank indicates how much "information" the matrix encodes and is vital in solving systems of linear equations.

1. Full Rank: A matrix is full rank if its rank equals the smaller of its row or column count. This implies all rows or columns are linearly independent.
2. Rank Deficiency: A rank-deficient matrix has rows or columns that are linear combinations of others.

Implementation ​

Building on the existing row_echelon method, we can efficiently calculate the rank by counting the non-zero rows in the RREF:

impl<K: Scalar> Matrix<K> {
    pub fn rank(&self) -> usize {
        let rref = self.row_echelon();
        rref.data
            .iter()
            .filter(|row| row.iter().any(|&val| Into::<f32>::into(val).abs() > Self::TOLERANCE))
            .count()
    }
}

Key features:

• Row-Echelon Form: Ensures each pivot column has leading 1, making rank determination straightforward.
• Numerical Tolerance: Avoids false negatives due to floating-point inaccuracies by treating values below 1e-10 as zero.
• Time Complexity: O(n³) due to RREF computation.
• Space Complexity: O(n²) for intermediate storage.

Examples ​

Full Rank Matrix ​

let m = Matrix::from([[1.0, 0.0], [0.0, 1.0]]);
assert_eq!(m.rank(), 2);

Rank Deficiency ​

let m = Matrix::from([[1.0, 2.0], [2.0, 4.0]]);
assert_eq!(m.rank(), 1);

Zero Matrix ​

let m = Matrix::from([[0.0, 0.0], [0.0, 0.0]]);
assert_eq!(m.rank(), 0);

Applications ​

1. Solving Linear Systems:

   • The rank determines the solvability of (Ax = b).
   • If rank(A) = rank([A|b]), the system is consistent.

2. Matrix Decompositions:

   • Decompositions like Singular Value Decomposition (SVD) rely on rank to determine the number of significant components.

3. Dimension Reduction:

   • In data science, rank helps in techniques like Principal Component Analysis (PCA).

4. Graph Theory:

   • Rank of adjacency matrices relates to graph properties like connectivity.

The rank computation encapsulates a critical aspect of matrix analysis, showcasing the power of combining robust numerical techniques with clear mathematical foundations.

Exercise 14: Building a 3D Projection Matrix ​

The projection matrix is a fundamental component in 3D graphics pipelines that transforms 3D points in camera space into 2D points on the screen. Let's dive deep into how we implemented this in Rust.

Understanding Projection Matrices ​

A projection matrix transforms 3D points using perspective division to create the illusion of depth on a 2D screen. The transformation pipeline looks like:

World Space → Camera Space → Screen Space
     ↓             ↓            ↓
  (Model)       (View)     (Projection)

Implementation Details ​

Let's break down our implementation:

pub fn projection(fov: f32, ratio: f32, near: f32, far: f32) -> Matrix<f32> {
    // Input validation
    if near <= 0.0 || far <= 0.0 || near >= far {
        panic!("Invalid near and far clipping planes.");
    }

    let tan_half_fov = (fov / 2.0).tan();
    let depth = far - near;

    // Calculate frustum dimensions
    let t = near * tan_half_fov;     // top
    let b = -t;                      // bottom
    let r = t * ratio;               // right
    let l = -r;                      // left

    // Build the matrix
    let mut data = vec![vec![0.0; 4]; 4];

    // Scale factors
    data[0][0] = 2.0 * near / (r - l);  // X scale
    data[1][1] = 2.0 * near / (t - b);  // Y scale

    // Perspective transformation
    data[2][0] = (r + l) / (r - l);     // X offset
    data[2][1] = (t + b) / (t - b);     // Y offset
    data[2][2] = -(far + near) / depth; // Z mapping
    data[2][3] = -2.0 * far * near / depth; // W factor

    data[3][2] = -1.0;  // Perspective division

    Matrix::new(data)
}

Key Components ​

1. Field of View (fov)

   • Controls the angle of view
   • Wider FOV = more perspective distortion
   • tan_half_fov scales vertical dimensions

2. Aspect Ratio

   • Maintains proper proportions
   • Typically matches screen dimensions (e.g., 16:9)
   • Scales horizontal dimensions

3. Clipping Planes

   • near: Closest visible distance
   • far: Farthest visible distance
   • Objects between these remain visible

Matrix Structure ​

The resulting 4×4 matrix looks like:

┌                                  ┐
│ 2n/(r-l)   0        (r+l)/(r-l)   0     │
│   0      2n/(t-b)   (t+b)/(t-b)   0     │
│   0        0       -(f+n)/d    -2fn/d   │
│   0        0           -1         0      │
└                                  ┘

Where:

• n = near plane
• f = far plane
• d = far - near
• r,l,t,b = frustum bounds

Testing Strategy ​

We test three key aspects:

1. Normal Cases: Valid parameters produce correct matrices
2. Edge Cases: Near/far plane validation
3. Field of View: Different FOVs affect perspective

#[test]
fn test_projection_valid() {
    let proj = projection(90.0_f32.to_radians(), 16.0/9.0, 0.1, 100.0);
    assert!(proj.data[0][0] > 0.0);  // X scale positive
    assert!(proj.data[1][1] > 0.0);  // Y scale positive
    assert!(proj.data[2][2] < 0.0);  // Z mapping negative
    assert!(proj.data[3][2] < 0.0);  // W factor negative
}

Visual Validation ​

To verify the projection visually:

1. Lower FOV = More zoom effect
2. Higher FOV = More "fish-eye" effect
3. Near/far ratios affect depth precision

The provided display utility can help visualize these effects.

Usage Examples ​

Setting up a basic perspective projection:

// 90° FOV, 16:9 aspect ratio, 0.1 near, 100.0 far
let proj = projection(
    90.0_f32.to_radians(),
    16.0 / 9.0,
    0.1,
    100.0
);

// For VR-style wide view
let wide_proj = projection(
    120.0_f32.to_radians(),
    16.0 / 9.0,
    0.1,
    100.0
);

// For telephoto style narrow view
let narrow_proj = projection(
    30.0_f32.to_radians(),
    16.0 / 9.0,
    0.1,
    100.0
);

Common Issues and Solutions ​

1. Z-Fighting

   • Symptom: Flickering surfaces
   • Cause: Near/far ratio too large
   • Solution: Keep near/far ratio under 1000:1

2. Distortion

   • Symptom: Objects look stretched
   • Cause: Extreme aspect ratios
   • Solution: Maintain reasonable ratios

3. Precision Loss

   • Symptom: Jittery far objects
   • Cause: Large depth range
   • Solution: Use appropriate near/far values

Further Reading ​

1. OpenGL Perspective Matrix
2. Normalized Device Coordinates (NDC)
3. Homogeneous Coordinates
4. View Frustum

Exercise 15 - Complex Vector Spaces ​

Our journey through linear algebra culminates in a fascinating extension: implementing the entire library for complex numbers. This exercise ties together everything we've learned while revealing deep connections between complex arithmetic and linear transformations.

Working with Complex Numbers ​

The foundation of our implementation is the Complex type, which represents complex numbers using their real and imaginary components:

pub struct Complex {
    pub r: f32,    // Real part
    pub i: f32,    // Imaginary part
}

I chose 'r' and 'i' as field names for their clarity and common usage in mathematical contexts. The implementation provides standard operations like addition and multiplication, along with special complex operations like conjugate and modulus calculations.

Extending Our Generic Code ​

The real test of our library's design came when we made it work with complex numbers. Our careful use of Rust's trait system paid off - most algorithms worked without modification. The Scalar trait provided a clean interface for implementing complex arithmetic:

impl Scalar for Complex {
    fn zero() -> Self { Complex::new(0.0, 0.0) }
    fn one() -> Self { Complex::new(1.0, 0.0) }
    fn fma(a: Self, b: Self, c: Self) -> Self {
        (a * b) + c // Complex FMA using operator overloading
    }
}

Numerical Stability in Complex Space ​

Complex arithmetic introduces new numerical stability concerns. For example, computing the norm of a complex vector requires careful handling of the modulus:

impl Into<f32> for Complex {
    fn into(self) -> f32 {
        self.modulus()  // For norm calculations
    }
}

This conversion enables our vector and matrix norms to work seamlessly with complex values while maintaining numerical precision.

Applications and Insights ​

Working with complex vector spaces opened up applications in:

• Signal Processing: Fourier transforms and filter design
• Quantum Mechanics: State vectors and operators
• Control Systems: Transfer functions and stability analysis

The implementation revealed how fundamental complex numbers are to many fields, and how a well-designed generic library can support them naturally.

Testing Complex Operations ​

Testing complex number operations required special attention to edge cases and known mathematical identities. Our test suite verifies:

• Basic arithmetic properties
• Special cases like division by zero
• Complex vector operations
• Matrix transformations with complex entries

Project Review and Conclusions ​

This linear algebra implementation journey taught valuable lessons about both mathematics and software design.

Technical Achievements ​

The library successfully implements:

• Fundamental vector and matrix operations
• Advanced linear algebra concepts
• Generic support for both real and complex numbers
• Numerically stable algorithms
• Clean, well-documented interfaces

Design Philosophy ​

Our design choices prioritized:

• Type safety through Rust's strong type system
• Generic programming using traits
• Clear error handling
• Performance through careful algorithm selection
• Maintainability through modular design

Future Directions ​

The library could be extended with:

• Eigenvalue calculations for complex matrices
• Matrix decompositions like SVD
• Parallel processing for large matrices
• Additional numeric types beyond real and complex
• More specialized complex number operations

Learning Outcomes ​

Building this library has provided deep insights into:

• The connection between abstract math and concrete code
• The power of generic programming
• The importance of numerical stability
• The value of comprehensive testing
• The benefits of good API design

Most importantly, we've seen how theoretical concepts in linear algebra translate into practical, reusable code. The exercise demonstrated that with careful design, we can create abstractions that are both mathematically correct and practically useful.

Final Thoughts ​

The journey from basic vector operations to complex vector spaces has been enlightening. We've built not just a library, but a deeper understanding of how linear algebra concepts manifest in code.

This project stands as a testament to Rust's ability to express mathematical concepts clearly and efficiently, while maintaining the performance and safety guarantees needed for numerical computing. Whether used for graphics, physics simulations, or machine learning, these fundamentals of linear algebra implementation will prove invaluable.

Resources for Further Learning ​

For those interested in diving deeper:

• "Linear Algebra and its Applications" by Gilbert Strang
• "Visual Complex Analysis" by Tristan Needham
• The Rust Programming Language book
• "Matrix Computations" by Golub and Van Loan

These resources can help expand on the concepts we've explored and open new avenues for further study and application.