Part 2: Linear Algebra Fundamentals

Part of the Mathematics for Programming 101 Series

Why Linear Algebra Clicked for Me

I was implementing image preprocessing for a computer vision model. The code was a mess—nested loops iterating over pixels, converting RGB to grayscale, normalizing values, rotating images.

Hundreds of lines. Slow execution. Hard to maintain.

Then I learned about matrix operations:

# Before: Nested loops
for i in range(height):
    for j in range(width):
        gray[i][j] = 0.299*r[i][j] + 0.587*g[i][j] + 0.114*b[i][j]

# After: One line with matrix multiplication
gray = image @ np.array([0.299, 0.587, 0.114])

10x faster. One line instead of hundreds. Easier to understand.

That's when I realized: linear algebra isn't abstract theory—it's the language of data transformation.

What is Linear Algebra?

Linear algebra is the mathematics of vectors and matrices. In programming terms:

• Vectors = Arrays of numbers (like lists)

• Matrices = 2D arrays (like nested lists)

• Operations = Transforming data efficiently

Every time you:

• Process images or video

• Train machine learning models

• Transform 3D graphics

• Manipulate datasets

• Solve systems of equations

You're using linear algebra.

Vectors: The Foundation

Understanding Vectors

A vector is an ordered collection of numbers. Think of it as:

• A point in space

• A direction with magnitude

• A list of features

import numpy as np

# Different ways to think about vectors
position = np.array([3, 4])  # Point at x=3, y=4
direction = np.array([1, 1])  # Direction: northeast at 45°
features = np.array([age, income, score])  # ML features

# Vector in 3D space
point_3d = np.array([1, 2, 3])

Vector Operations

# Vector addition: combine movements/forces
v1 = np.array([1, 2])
v2 = np.array([3, 1])
result = v1 + v2  # [4, 3]

print(f"v1 + v2 = {result}")

# Scalar multiplication: scaling
v = np.array([1, 2, 3])
scaled = 3 * v  # [3, 6, 9]

print(f"3 * v = {scaled}")

# Dot product: measure similarity/projection
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
dot = np.dot(a, b)  # 1*4 + 2*5 + 3*6 = 32

print(f"a · b = {dot}")

# Magnitude (length): Euclidean distance
magnitude = np.linalg.norm(a)  # sqrt(1² + 2² + 3²) = sqrt(14)

print(f"||a|| = {magnitude:.3f}")

# Unit vector: direction without magnitude
unit = a / np.linalg.norm(a)
print(f"Unit vector of a: {unit}")
print(f"Magnitude of unit vector: {np.linalg.norm(unit):.3f}")

Real Use Case: Cosine Similarity

Finding similar documents or products:

def cosine_similarity(a, b):
    """Measure similarity between vectors (0 to 1)"""
    dot_product = np.dot(a, b)
    magnitude_a = np.linalg.norm(a)
    magnitude_b = np.linalg.norm(b)
    return dot_product / (magnitude_a * magnitude_b)

# Document vectors (simplified TF-IDF representation)
doc1 = np.array([0.5, 0.8, 0.0, 0.3])  # Features: [python, java, ruby, go]
doc2 = np.array([0.6, 0.7, 0.0, 0.2])  # Similar to doc1
doc3 = np.array([0.0, 0.1, 0.9, 0.0])  # Different (focuses on ruby)

similarity_1_2 = cosine_similarity(doc1, doc2)
similarity_1_3 = cosine_similarity(doc1, doc3)

print(f"Doc1 vs Doc2 similarity: {similarity_1_2:.3f}")  # High similarity
print(f"Doc1 vs Doc3 similarity: {similarity_1_3:.3f}")  # Low similarity

Used in:

• Recommendation systems

• Document search

• Image similarity

• Face recognition

Matrices: Data Transformations

Understanding Matrices

A matrix is a 2D array of numbers. Think of it as:

• A transformation function

• A dataset (rows = samples, columns = features)

• A collection of vectors

# Matrix as dataset
data = np.array([
    [25, 50000, 5.0],  # Person 1: age, income, rating
    [30, 60000, 4.5],  # Person 2
    [35, 75000, 4.8],  # Person 3
])

print("Dataset shape:", data.shape)  # (3, 3) = 3 people, 3 features

# Matrix as transformation
rotation = np.array([
    [0, -1],  # 90-degree rotation matrix
    [1,  0]
])

Matrix Operations

# Matrix addition
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
C = A + B

print("A + B =\n", C)

# Scalar multiplication
D = 2 * A
print("2 * A =\n", D)

# Matrix multiplication (THE IMPORTANT ONE)
E = np.array([[1, 2], [3, 4]])
F = np.array([[5, 6], [7, 8]])
G = E @ F  # Python 3.5+ matrix multiplication operator

print("E @ F =\n", G)

# Transpose: flip rows and columns
H = np.array([[1, 2, 3], [4, 5, 6]])
H_T = H.T

print("H =\n", H)
print("H^T =\n", H_T)

Matrix Multiplication Intuition

Matrix multiplication transforms vectors:

# Transformation matrix
transform = np.array([
    [2, 0],  # Scale x by 2
    [0, 3]   # Scale y by 3
])

# Original point
point = np.array([1, 1])

# Transform the point
transformed = transform @ point  # [2, 3]

print(f"Original point: {point}")
print(f"Transformed point: {transformed}")

# Multiple transformations: compose by multiplying matrices
rotation_90 = np.array([[0, -1], [1, 0]])
scale_2x = np.array([[2, 0], [0, 2]])

# Apply both: first scale, then rotate
combined = rotation_90 @ scale_2x
result = combined @ point

print(f"Scaled then rotated: {result}")

Real-World Application: Image Processing

Grayscale Conversion

# RGB to grayscale using matrix multiplication
def rgb_to_grayscale(image):
    """Convert RGB image to grayscale using standard weights"""
    # Standard weights: R=0.299, G=0.587, B=0.114
    weights = np.array([0.299, 0.587, 0.114])
    
    # Matrix multiplication: (H×W×3) @ (3×1) = (H×W×1)
    grayscale = image @ weights
    
    return grayscale

# Example: 2x2 RGB image
rgb_image = np.array([
    [[255, 0, 0], [0, 255, 0]],    # Red, Green
    [[0, 0, 255], [255, 255, 255]]  # Blue, White
], dtype=float)

gray = rgb_to_grayscale(rgb_image)
print("Grayscale values:\n", gray)

Image Rotation

def rotate_image(image, angle):
    """Rotate image by angle (in degrees)"""
    # Convert to radians
    theta = np.radians(angle)
    
    # Rotation matrix
    cos_theta = np.cos(theta)
    sin_theta = np.sin(theta)
    rotation_matrix = np.array([
        [cos_theta, -sin_theta],
        [sin_theta, cos_theta]
    ])
    
    # Apply rotation to each pixel coordinate
    height, width = image.shape[:2]
    center = np.array([width//2, height//2])
    
    # This is simplified; real implementation uses interpolation
    # and handles image boundaries properly
    return rotation_matrix, center

# Example usage
rotation_matrix, center = rotate_image(rgb_image, 45)
print(f"45-degree rotation matrix:\n{rotation_matrix}")

Machine Learning Applications

Linear Regression

def linear_regression_matrix_form(X, y):
    """
    Solve linear regression using normal equation: θ = (X^T X)^(-1) X^T y
    
    This is pure linear algebra!
    """
    # Add bias term (column of ones)
    X_b = np.c_[np.ones((X.shape[0], 1)), X]
    
    # Normal equation: closed-form solution
    # θ = (X^T X)^(-1) X^T y
    theta = np.linalg.inv(X_b.T @ X_b) @ X_b.T @ y
    
    return theta

# Generate sample data
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 4 + 3 * X.flatten() + np.random.randn(100)

# Solve using linear algebra
theta = linear_regression_matrix_form(X, y)

print(f"Intercept: {theta[0]:.3f}")  # Should be ~4
print(f"Slope: {theta[1]:.3f}")      # Should be ~3

Neural Network Forward Pass

def neural_network_forward(X, W1, b1, W2, b2):
    """
    Two-layer neural network using matrix operations
    
    X: input matrix (batch_size × input_dim)
    W1: first layer weights (input_dim × hidden_dim)
    b1: first layer bias (hidden_dim,)
    W2: second layer weights (hidden_dim × output_dim)
    b2: second layer bias (output_dim,)
    """
    # First layer: Z1 = X @ W1 + b1
    Z1 = X @ W1 + b1
    
    # Activation: ReLU(Z1)
    A1 = np.maximum(0, Z1)
    
    # Second layer: Z2 = A1 @ W2 + b2
    Z2 = A1 @ W2 + b2
    
    # Output activation: Softmax
    exp_Z2 = np.exp(Z2 - np.max(Z2, axis=1, keepdims=True))
    A2 = exp_Z2 / np.sum(exp_Z2, axis=1, keepdims=True)
    
    return A2

# Example: 4 samples, 3 features, 5 hidden units, 2 classes
X = np.random.randn(4, 3)
W1 = np.random.randn(3, 5) * 0.01
b1 = np.zeros(5)
W2 = np.random.randn(5, 2) * 0.01
b2 = np.zeros(2)

predictions = neural_network_forward(X, W1, b1, W2, b2)
print("Predictions shape:", predictions.shape)  # (4, 2)
print("First sample probabilities:", predictions[0])
print("Sum to 1?", np.sum(predictions[0]))  # Should be 1.0

Eigenvalues and Eigenvectors

The Concept

An eigenvector of a matrix is a vector that only gets scaled (not rotated) when the matrix is applied to it.

# Matrix that scales in different directions
A = np.array([
    [3, 1],
    [0, 2]
])

# Find eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)

print("Eigenvalues:", eigenvalues)
print("Eigenvectors:\n", eigenvectors)

# Verify: A @ v = λ @ v
for i in range(len(eigenvalues)):
    v = eigenvectors[:, i]
    lambda_i = eigenvalues[i]
    
    left = A @ v
    right = lambda_i * v
    
    print(f"\nEigenvector {i+1}:")
    print(f"A @ v = {left}")
    print(f"λ * v = {right}")
    print(f"Close? {np.allclose(left, right)}")

Real Use Case: Principal Component Analysis (PCA)

Dimensionality reduction for machine learning:

def pca(X, n_components):
    """
    Perform PCA using eigendecomposition
    
    Returns the transformed data with reduced dimensions
    """
    # Center the data
    X_centered = X - np.mean(X, axis=0)
    
    # Compute covariance matrix
    cov_matrix = np.cov(X_centered.T)
    
    # Eigendecomposition
    eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
    
    # Sort by eigenvalue (descending)
    idx = eigenvalues.argsort()[::-1]
    eigenvalues = eigenvalues[idx]
    eigenvectors = eigenvectors[:, idx]
    
    # Select top n_components
    principal_components = eigenvectors[:, :n_components]
    
    # Transform data
    X_transformed = X_centered @ principal_components
    
    # Explained variance
    explained_variance = eigenvalues[:n_components] / np.sum(eigenvalues)
    
    return X_transformed, explained_variance

# Example: Reduce 4D data to 2D
np.random.seed(42)
X = np.random.randn(100, 4)

X_reduced, variance = pca(X, n_components=2)

print("Original shape:", X.shape)       # (100, 4)
print("Reduced shape:", X_reduced.shape) # (100, 2)
print("Explained variance:", variance)
print("Total variance captured:", np.sum(variance))

Why this matters:

• Reduce features for faster ML training

• Visualize high-dimensional data

• Remove noise and redundancy

• Compress data efficiently

Practical Tips

1. Think in Transformations

Don't memorize matrix multiplication rules. Think: "This matrix transforms my data from space A to space B."

# Transformation thinking
embedding_matrix = np.random.randn(vocab_size, embedding_dim)

# word_index → embedded_vector
# This matrix transforms discrete word indices to continuous vectors
word_embedding = embedding_matrix[word_index]

2. Use Broadcasting

NumPy broadcasting avoids explicit loops:

# Add bias to all samples at once
batch = np.random.randn(32, 128)  # 32 samples, 128 features
bias = np.random.randn(128)       # One bias per feature

# Broadcasting: automatically repeats bias for each sample
result = batch + bias  # Works! Shape: (32, 128)

3. Vectorize Operations

Replace loops with matrix operations:

# Slow: loop version
def compute_distances_slow(points1, points2):
    distances = []
    for p1 in points1:
        row = []
        for p2 in points2:
            dist = np.sqrt(np.sum((p1 - p2)**2))
            row.append(dist)
        distances.append(row)
    return np.array(distances)

# Fast: vectorized version
def compute_distances_fast(points1, points2):
    # ||a - b||² = ||a||² + ||b||² - 2(a·b)
    p1_sq = np.sum(points1**2, axis=1, keepdims=True)
    p2_sq = np.sum(points2**2, axis=1, keepdims=True)
    cross = points1 @ points2.T
    
    distances_sq = p1_sq + p2_sq.T - 2 * cross
    return np.sqrt(np.maximum(distances_sq, 0))

# Benchmark
points1 = np.random.randn(100, 50)
points2 = np.random.randn(80, 50)

import time

start = time.time()
slow_result = compute_distances_slow(points1, points2)
slow_time = time.time() - start

start = time.time()
fast_result = compute_distances_fast(points1, points2)
fast_time = time.time() - start

print(f"Slow version: {slow_time:.4f}s")
print(f"Fast version: {fast_time:.4f}s")
print(f"Speedup: {slow_time/fast_time:.1f}x")

When You Need Linear Algebra

You'll use these concepts when:

1. Machine Learning: Every ML algorithm uses matrix operations

2. Computer Graphics: Transformations, projections, rendering

3. Data Processing: Feature engineering, normalization, PCA

4. Computer Vision: Image filters, transformations, CNNs

5. Natural Language Processing: Word embeddings, transformers

6. Scientific Computing: Simulations, numerical methods

7. Game Development: Physics, transformations, camera systems

Key Takeaways

• Vectors represent data points, directions, or features

• Matrices represent transformations or datasets

• Matrix multiplication transforms data efficiently

• Eigendecomposition reveals the fundamental directions of transformation

• Vectorization replaces slow loops with fast matrix operations

• NumPy makes linear algebra practical and fast in Python

What's Next

In the next article, we'll explore calculus and optimization—the mathematics behind training machine learning models, understanding gradients, and building neural networks from scratch.

You'll learn:

• How derivatives power gradient descent

• Building backpropagation from scratch

• Optimization algorithms explained

• Real debugging with calculus

Continue to Part 3: Calculus and Optimization →

---

Navigation

• ← Part 1: Why Mathematics Matters

• Back to Series Overview

• Part 3: Calculus and Optimization →