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.

Understanding the Core Types ​

Before diving into the exercises, let's examine the three fundamental types that form the foundation of our linear algebra library:

The Scalar Trait ​

The Scalar trait defines what operations a numeric type must support to 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)
    PartialEq +            // a == b
    PartialOrd             // a < b
{
    fn zero() -> Self;      // Identity for addition
    fn one() -> Self;       // Identity for multiplication
    fn to_f32(&self) -> f32; // Convert for normalization
    fn fma(a: Self, b: Self, c: Self) -> Self; // (a * b) + c
}

The trait includes:

• Basic arithmetic operations
• Assignment operations
• Comparison operations
• Identity elements (zero and one)
• Conversion to f32 for normalization
• FMA (Fused Multiply-Add) for better numerical precision

The Vector Type ​

pub struct Vector<K: Scalar> {
    pub data: Vec<K>,
}

Represents a mathematical vector as a single-dimensional array of scalar values. The generic type K must implement the Scalar trait, allowing the vector to work with any compatible numeric type.

The Matrix Type ​

pub struct Matrix<K: Scalar> {
    pub data: Vec<Vec<K>>,
}

Represents a matrix as a two-dimensional array of scalar values. Like Vector, it uses a generic type K that implements Scalar.

These three types work together to provide a foundation for all linear algebra operations in the library:

• Scalar provides the numeric operations
• Vector handles one-dimensional collections
• Matrix manages two-dimensional collections

Each exercise builds on these types, implementing increasingly complex linear algebra operations while maintaining type safety and numerical stability.

In place or new instances? ​

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 the fundamental operations: vector/matrix addition, subtraction, and scalar multiplication. Each operation follows the principle of modifying the data in place to maintain optimal memory usage.

The key insight here is keeping things simple but efficient. Both time and space complexity stay within O(n) - we make a single pass through the elements with only the memory needed for the operation result. For example, when adding vectors, we're just going element by element, adding the corresponding values. The same applies for matrices, just with one more dimension to consider.

Addition creates a new sum by combining corresponding elements. Picture two 2D vectors [1,2] and [3,4] - adding them gives us [4,6]. For matrices, we do the same thing position by position. What's neat is this operation preserves the structure - add two 2x2 matrices and you get another 2x2 matrix.

Subtraction works similarly but takes away values instead of combining them. That same [1,2] and [3,4] example under subtraction gives [-2,-2]. The operation lets us find the difference between vectors/matrices, which has applications in everything from error calculations to finding relative positions.

Scaling multiplies every element by a single number (scalar). If you scale [1,2] by 2, you get [2,4]. This operation stretches or shrinks the vector/matrix uniformly. A scale factor of 0 zeros out everything, while 1 leaves it unchanged.

The beauty of these operations is their simplicity - they're building blocks that combine to enable far more complex transformations later. Each maintains clear mathematical properties while being computationally straightforward.

All this is achieved through direct arithmetic - no complex math functions needed. Just addition, subtraction, and multiplication operating systematically on the data structures.

Exercise 01 - Linear combination ​

A linear combination takes multiple vectors and their corresponding coefficients to produce a new vector through multiplication and addition. For example, given vectors v1=[1,0] and v2=[0,1] with coefficients 2 and 3, we multiply each vector by its coefficient (getting [2,0] and [0,3]) then add the results to get [2,3].

The implementation maintains O(n) complexity by making a single pass through the vectors, using Fused Multiply-Add (FMA) operations for both better performance and numerical precision. FMA combines the coefficient multiplication and addition in one step, reducing rounding errors that could accumulate with separate operations.

The function validates matching vector sizes and coefficient counts before proceeding - for example, you can't combine a 2D vector with a 3D vector, or have more coefficients than vectors. This size checking ensures mathematical correctness.

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.

Exercise 02 - Linear interpolation ​

Linear interpolation (lerp) finds values along a straight line between two endpoints. Think of it like a slider - when t=0 you're at the start, t=1 you're at the end, and t=0.5 puts you exactly halfway. This works for both single numbers and entire vectors.

The implementation stays efficient with O(n) time complexity by touching each element just once. The formula a + t(b-a) requires only basic arithmetic, no complex math functions. To ensure numeric stability, we use a single scalar multiplication rather than the weighted average form (1-t)a + tb.

What makes this function powerful is its generic nature - it works with any type that supports basic arithmetic, from simple numbers to matrices. The only requirement is that the type can be scaled and added.

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.

Exercise 03 - Dot product ​

The dot product of two vectors reveals how aligned they are. By multiplying corresponding elements and summing the results, we get a single number that tells us a lot about the vectors' relationship.

The function maintains O(n) complexity by processing each pair of elements just once. Using Fused Multiply-Add (FMA) operations keeps things numerically stable - instead of separate multiplications and additions that might lose precision, FMA handles both in one step.

The dot product gives us some interesting insights: perpendicular vectors (like [1,0] and [0,1]) have a dot product of 0, while vectors pointing in the same direction give a positive result proportional to their lengths. This makes it perfect for checking angles between vectors - an operation we'll need in later exercises.

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

Exercise 04 - Norms ​

A vector norm measures "length" in different ways. Our implementation covers three key norms, each with O(n) complexity as we only need one pass through the vector:

The Manhattan norm (1-norm) sums absolute values - imagine walking through a city grid, where you can only move along streets. For a vector [3,4], you walk 7 blocks total (3 + 4).

The Euclidean norm (2-norm) is our familiar straight-line distance - think "as the crow flies". For that same vector [3,4], Pythagoras tells us it's 5 units long.

The Supremum norm (infinity-norm) takes the largest absolute component. In [3,4], that's 4. Think of it as the minimum box size needed to contain the vector.

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

An interesting property these share: for any vector, infinity-norm ≤ 2-norm ≤ 1-norm. Our city walker always travels further than the crow flies!

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.

Exercise 05 - Cosine ​

The angle cosine function reveals how aligned two vectors are by computing their dot product divided by the product of their lengths. The result is beautifully bounded between -1 and 1, where:

1 means vectors point the same way, 0 means they're perpendicular, and -1 means they point opposite directions. For example, if you have [1,0] and [0,1], they're perpendicular, so their angle cosine is 0.

The implementation maintains O(n) complexity by reusing our dot product and norm functions. We carefully avoid dividing by zero - if either vector has zero length, the angle is undefined. Using FMA operations throughout helps maintain numerical precision.

Mathematical Foundation ​

The cosine formula between vectors is:

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

Where:

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

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

Exercise 06 - Cross product ​

The cross product creates a new vector perpendicular to two input vectors, essential for 3D geometry. Unlike our previous operations that worked in any dimension, this is specifically for 3D vectors.

What makes it special is that the result tells us two things: its direction is perpendicular to both input vectors (following the right-hand rule), and its length represents the area of the parallelogram they form.

For instance, crossing unit vectors follows a neat pattern: [1,0,0] × [0,1,0] = [0,0,1]. This is why i × j = k in physics. The operation is also anticommutative - swap the inputs and you get the opposite result.

The implementation uses FMA for precision in the component calculations, each being a determinant of a 2×2 slice of the input vectors. We always check for 3D vectors first - trying to take a cross product in 2D or 4D wouldn't make geometric sense.

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

Exercise 07 - Linear map, Matrix multiplication ​

A matrix can act as a transformation machine - feed it a vector, and it maps that vector to a new position. This forms the foundation of computer graphics and physics simulations. We handle both matrix-vector multiplication (applying the transform) and matrix-matrix multiplication (combining transforms).

For example, feeding vector [4,2] through the matrix [[2,0],[0,2]] doubles each component, giving [8,4]. The identity matrix [[1,0],[0,1]] leaves vectors unchanged. These operations require careful dimension matching - the matrix's column count must match the input vector's size.

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

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

Exercise 08 - Trace ​

The trace is one of the simplest yet most useful matrix operations - just sum up the diagonal elements. For a 2×2 matrix [[a,b],[c,d]], the trace is a+d. It requires the matrix to be square since we need matching diagonal positions.

The operation maintains O(n) time complexity since we only need to visit n diagonal elements in an n×n matrix. Some handy properties emerge: the trace of an identity matrix equals the matrix dimension (each diagonal is 1), while a zero matrix gives trace zero.

A neat mathematical insight: the trace tells us the sum of a matrix's eigenvalues without having to calculate them individually. This makes it valuable for checking matrix properties quickly.

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ₙₙ

Exercise 09 - Transpose ​

Transposing a matrix flips it across its main diagonal - rows become columns and columns become rows. For example, [[1,2],[3,4]] becomes [[1,3],[2,4]]. Unlike many previous operations, transpose works on rectangular matrices, not just square ones.

The implementation maintains O(nm) complexity by visiting each element exactly once for both time and space. What makes transpose interesting is how some matrices behave: identity matrices stay unchanged (symmetric), while for any matrix A, transposing twice returns the original (A^T^T = A).

This operation is particularly important for converting between column-major and row-major matrix formats, which matters in graphics programming and numerical computing.

Mathematical Definition ​

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

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 ​

Reduced Row Echelon Form (RREF) takes row echelon form one step further. Besides the stair-step pattern, it ensures:

1. Each pivot is exactly 1
2. Pivots are the only non-zero entry in their column

The implementation uses Gaussian-Jordan elimination with partial pivoting - choosing the largest value in each column as pivot for stability. After getting a pivot, we scale its row to get a leading 1, then eliminate entries both below AND above it. This gives us O(n³) time complexity.

RREF is particularly useful because it gives us a unique canonical form for each matrix. When solving systems of equations, this means we can read the solutions directly. For example, if a variable's column reduces to [1,0,0], that variable equals the value in the rightmost column of that row.

Implementation Strategy ​

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

Exercise 11 - Determinant ​

The determinant measures how a matrix transformation affects area (2D) or volume (3D). A determinant of 2 means the transformation doubles the area/volume, while 0 means the transformation squishes space into a lower dimension (like flattening a cube into a square).

For efficiency up to 4x4 matrices, we use specialized formulas per dimension:

• 2x2: ad - bc for matrix [[a,b],[c,d]]
• 3x3: Cofactor expansion along first row
• 4x4: Cofactor expansion with 3x3 minors

A key insight - when det(A) = 0, the matrix is singular (non-invertible). This means some vector gets transformed to zero, effectively losing information. Geometrically, this collapses space onto a line or point.

The implementation achieves O(n³) complexity through recursion on cofactors, using FMA operations for numerical stability particularly in the base 2x2 case.

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

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 - Inverse ​

A matrix inverse undoes a transformation - multiplying a matrix by its inverse gives you the identity matrix. Not every matrix has an inverse though. When det(A)=0, the matrix is singular (non-invertible) because the transformation loses information that can't be recovered.

The implementation uses the classical adjugate method:

1. Calculate determinant - if too close to zero, matrix is singular
2. Compute cofactors (signed minors) for each element
3. Transpose to get adjugate matrix
4. Divide adjugate by determinant

Since we work with floating point numbers, we use a small tolerance value to detect "effectively zero" determinants. The operation requires O(n³) time due to cofactor calculations, and O(n²) space for the result matrix.

Key test: multiplying a matrix by its computed inverse should give identity matrix (within floating point precision).

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.

Mathematical Foundation ​

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).

Properties ​

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

Exercise 13 - Rank ​

The rank represents the number of linearly independent rows or columns in a matrix. In practical terms, it reveals three crucial things:

1. The dimension of the image/output space - a rank 2 matrix maps vectors into a 2D space, even if the matrix is larger
2. How much information the matrix preserves - full rank means no information loss, lower rank means some dimensions get collapsed
3. Whether we have redundant equations in a linear system - rank less than the number of equations means some are combinations of others

Our implementation leverages RREF to calculate rank - after getting reduced row echelon form, we count non-zero rows. This works because RREF makes dependencies obvious - any linearly dependent rows become all zeros during reduction.

For example, in a system:

• Rank equals rows (max rank): Every equation gives new information
• Rank less than rows: Some equations are redundant
• Rank = 0: All equations are trivial (0 = 0)

Exercise 14 - Bonus: Projection matrix ​

The perspective projection matrix is fundamental to 3D graphics, transforming 3D points into 2D screen coordinates while maintaining realistic perspective effects. Understanding each component reveals how it creates the illusion of depth.

Resources ​

• Projecting 3D Points into a 2D Screen
• The Perspective and Orthographic Projection Matrix

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)

The View Frustum ​

Think of the frustum as a pyramid with its top cut off:

• Near plane: Closest visible distance (anything closer is clipped)
• Far plane: Farthest visible distance
• FOV (Field of View): Angular extent of scene visible to camera
• Aspect ratio: Screen width/height ratio (like 16:9)

Understanding Projection Matrix Mathematics ​

1. Homogeneous Coordinates ​

We use 4D vectors (x,y,z,w) to represent 3D points because:

• Enables perspective division (x/w, y/w, z/w)
• Makes translation linear (impossible with 3D vectors)
• Final w component stores -z for perspective division

2. Basic Trigonometry Relations ​

tan(fov/2) = (height/2) / near
height = 2 * near * tan(fov/2)
width = height * aspect_ratio

3. Matrix Components Explained ​

[2n/(r-l)    0           (r+l)/(r-l)      0          ]
[0           2n/(t-b)    (t+b)/(t-b)      0          ]
[0           0           -(f+n)/(f-n)     -2fn/(f-n) ]
[0           0           -1               0          ]

Each row serves a specific purpose:

1. X-scaling & offset: Maps view width to normalized device coordinates
2. Y-scaling & offset: Maps view height while maintaining aspect ratio
3. Z-mapping: Converts depth to [0,1] range for z-buffer
4. W-component: Creates perspective divide effect (far objects appear smaller)

4. Key Transformations ​

• World → Camera space (ModelView)
• Camera → NDC space (Projection)
• NDC → Screen space (Viewport)

The beauty is how this single matrix combines frustum mapping, perspective effects, and depth handling in one operation that hardware can efficiently process.

FOV Effects ​

• Wide FOV (100°): Fish-eye effect, good for tight spaces
• Normal FOV (70°): Approximates human vision
• Narrow FOV (40°): Telephoto effect, compresses depth

Numerical Pipeline ​

1. World coordinates → Camera space (by view matrix)
2. Camera space → Clip space (by projection matrix)
3. Clip space → NDC (by perspective divide)
4. NDC → Screen coordinates (by viewport transform)

The projection matrix works by first creating a perspective divide effect through the z coordinate, then mapping the resulting points into a canonical view volume ([-1,1]³).

Optimized Projection Matrix ​

For symmetric view frustums (which is almost always the case), the calculations can be greatly simplified:

1. Scale Optimization

   Original:
   top = near * tan(fov/2)
   right = top * ratio
   x_scale = 2n/(r-l) = 1/(tan(fov/2))
   
   Simplified:
   x_scale = 1/tan(fov/2) directly

   This works because in a symmetric frustum, left=-right and bottom=-top.

2. Depth Mapping The depth calculations remain unchanged as they can't be further simplified without losing precision:

   depth_norm = (f+n)/(n-f)
   depth_map = 2fn/(n-f)

This optimization:

• Reduces floating point operations
• Maintains numerical stability
• Works for all valid parameter combinations
• Produces identical results to the full calculation

The only requirement is that the frustum remains symmetric (centered on the view axis), which is standard for most 3D graphics applications.

Further Reading ​

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

Exercise 15 - Bonus: Complex vector spaces ​

Our journey through linear algebra culminates in an interesting extension: implementing the entire library for complex numbers. This 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 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 { Self::new(0.0, 0.0) }
    fn one() -> Self { Self::new(1.0, 0.0) }

    /// Implement conversion from Complex to f32 by returning the modulus.
    fn to_f32(&self) -> f32 {
        self.modulus()
    }
}

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 ​

Well, we've made it through building a linear algebra library from scratch! It's been quite a journey from simple vector math to handling complex numbers and 3D projections. Let's break down what we've accomplished and what we learned along the way.

We managed to put together a solid set of tools that handle everything from basic vector operations to matrix transformations. The code isn't just theoretical - it's practical stuff you could use in graphics engines, physics simulations, or data analysis. Using Rust helped keep things safe and fast without getting in the way too much.

The neat thing about this project was seeing how the math actually turns into useful code. Those dot products we implemented? They're what you'd use to figure out angles in a game engine. The matrix operations? They're pushing pixels around in 3D graphics. Even the complex number support opens doors to signal processing and quantum mechanics applications.

There's still plenty we could add - things like eigenvalues, better matrix factorization, or parallel processing for bigger matrices. We could also make the error handling more Rust-like with proper Results and Options instead of just panicking when things go wrong.

But the big takeaway here is that we've built something that actually works and makes sense. The code is clean enough that you can come back to it later and still understand what's going on, and it's flexible enough to use in real projects.

This wasn't just about memorizing formulas or following patterns - it was about understanding how these pieces fit together to solve real problems. Whether you're making games, crunching numbers, or just learning the ropes of linear algebra, having code that you understand from the ground up is pretty valuable.