Vectors and Matrices

Introduction

Linear algebra is the mathematics of vectors and matrices. It's fundamental to computer graphics, machine learning, scientific computing, data compression, and quantum computing. This reading introduces vectors, matrices, and their operations.

Learning Objectives

By the end of this reading, you will be able to:

• Represent and manipulate vectors in code
• Perform matrix operations (addition, multiplication, transpose)
• Understand geometric interpretations of linear algebra
• Implement efficient matrix algorithms

1. Vectors

Definition

A vector is an ordered list of numbers. In n-dimensional space:

v = [v₁, v₂, ..., vₙ]

import math

class Vector:
    """n-dimensional vector"""

    def __init__(self, components):
        self.components = list(components)
        self.n = len(self.components)

    def __repr__(self):
        return f"Vector({self.components})"

    def __len__(self):
        return self.n

    def __getitem__(self, i):
        return self.components[i]

    def __setitem__(self, i, value):
        self.components[i] = value

    def __eq__(self, other):
        return self.components == other.components

    # Vector addition
    def __add__(self, other):
        if self.n != other.n:
            raise ValueError("Vectors must have same dimension")
        return Vector([a + b for a, b in zip(self.components, other.components)])

    # Vector subtraction
    def __sub__(self, other):
        if self.n != other.n:
            raise ValueError("Vectors must have same dimension")
        return Vector([a - b for a, b in zip(self.components, other.components)])

    # Scalar multiplication
    def __mul__(self, scalar):
        return Vector([scalar * x for x in self.components])

    def __rmul__(self, scalar):
        return self * scalar

    # Negation
    def __neg__(self):
        return Vector([-x for x in self.components])

# Examples
v = Vector([1, 2, 3])
w = Vector([4, 5, 6])
print(f"v + w = {v + w}")       # Vector([5, 7, 9])
print(f"2 * v = {2 * v}")       # Vector([2, 4, 6])
print(f"v - w = {v - w}")       # Vector([-3, -3, -3])

Dot Product

The dot product (inner product) of two vectors:

v · w = Σᵢ vᵢwᵢ = v₁w₁ + v₂w₂ + ... + vₙwₙ

def dot(self, other):
    """Dot product of two vectors"""
    if self.n != other.n:
        raise ValueError("Vectors must have same dimension")
    return sum(a * b for a, b in zip(self.components, other.components))

# Add to Vector class
Vector.dot = dot

v = Vector([1, 2, 3])
w = Vector([4, 5, 6])
print(f"v · w = {v.dot(w)}")  # 1*4 + 2*5 + 3*6 = 32

Vector Magnitude (Norm)

The Euclidean norm (length) of a vector:

||v|| = √(v · v) = √(v₁² + v₂² + ... + vₙ²)

def magnitude(self):
    """Euclidean norm (length) of vector"""
    return math.sqrt(sum(x**2 for x in self.components))

def normalize(self):
    """Return unit vector in same direction"""
    mag = self.magnitude()
    if mag == 0:
        raise ValueError("Cannot normalize zero vector")
    return Vector([x / mag for x in self.components])

Vector.magnitude = magnitude
Vector.normalize = normalize

v = Vector([3, 4])
print(f"||v|| = {v.magnitude()}")  # 5.0
print(f"unit v = {v.normalize()}")  # Vector([0.6, 0.8])

Geometric Interpretation

def angle_between(v, w):
    """Angle between two vectors in radians"""
    # cos(θ) = (v · w) / (||v|| ||w||)
    cos_theta = v.dot(w) / (v.magnitude() * w.magnitude())
    # Clamp to handle numerical errors
    cos_theta = max(-1, min(1, cos_theta))
    return math.acos(cos_theta)

def projection(v, w):
    """Project v onto w"""
    # proj_w(v) = ((v · w) / (w · w)) * w
    scalar = v.dot(w) / w.dot(w)
    return scalar * w

v = Vector([1, 2])
w = Vector([3, 0])
print(f"Angle: {math.degrees(angle_between(v, w)):.1f}°")  # 63.4°
print(f"Projection of v onto w: {projection(v, w)}")  # Vector([1, 0])

Cross Product (3D only)

def cross(v, w):
    """Cross product of 3D vectors"""
    if len(v) != 3 or len(w) != 3:
        raise ValueError("Cross product only defined for 3D vectors")

    return Vector([
        v[1] * w[2] - v[2] * w[1],
        v[2] * w[0] - v[0] * w[2],
        v[0] * w[1] - v[1] * w[0]
    ])

v = Vector([1, 0, 0])
w = Vector([0, 1, 0])
print(f"v × w = {cross(v, w)}")  # Vector([0, 0, 1])

2. Matrices

Definition

A matrix is a 2D array of numbers. An m×n matrix has m rows and n columns.

class Matrix:
    """m × n matrix"""

    def __init__(self, data):
        """Create matrix from 2D list"""
        self.data = [list(row) for row in data]
        self.rows = len(self.data)
        self.cols = len(self.data[0]) if self.data else 0

    def __repr__(self):
        rows = [f"  {row}" for row in self.data]
        return "Matrix([\n" + ",\n".join(rows) + "\n])"

    def __getitem__(self, idx):
        if isinstance(idx, tuple):
            i, j = idx
            return self.data[i][j]
        return self.data[idx]

    def __setitem__(self, idx, value):
        if isinstance(idx, tuple):
            i, j = idx
            self.data[i][j] = value
        else:
            self.data[idx] = value

    @property
    def shape(self):
        return (self.rows, self.cols)

    @staticmethod
    def zeros(rows, cols):
        """Create zero matrix"""
        return Matrix([[0] * cols for _ in range(rows)])

    @staticmethod
    def identity(n):
        """Create n×n identity matrix"""
        data = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
        return Matrix(data)

    def transpose(self):
        """Return transpose of matrix"""
        data = [[self.data[i][j] for i in range(self.rows)] for j in range(self.cols)]
        return Matrix(data)

# Examples
A = Matrix([
    [1, 2, 3],
    [4, 5, 6]
])
print(f"Shape: {A.shape}")  # (2, 3)
print(f"A[0,1] = {A[0,1]}")  # 2
print(f"A^T =\n{A.transpose()}")

I = Matrix.identity(3)
print(f"I =\n{I}")

Matrix Addition

def matrix_add(self, other):
    """Add two matrices"""
    if self.shape != other.shape:
        raise ValueError(f"Shape mismatch: {self.shape} vs {other.shape}")

    data = [[self[i,j] + other[i,j] for j in range(self.cols)]
            for i in range(self.rows)]
    return Matrix(data)

def scalar_multiply(self, scalar):
    """Multiply matrix by scalar"""
    data = [[scalar * self[i,j] for j in range(self.cols)]
            for i in range(self.rows)]
    return Matrix(data)

Matrix.__add__ = matrix_add
Matrix.__mul__ = scalar_multiply
Matrix.__rmul__ = scalar_multiply

Matrix Multiplication

Key rule: A (m×n) × B (n×p) = C (m×p)

C[i,j] = Σₖ A[i,k] × B[k,j]

def matrix_multiply(self, other):
    """Multiply two matrices"""
    if self.cols != other.rows:
        raise ValueError(f"Cannot multiply {self.shape} by {other.shape}")

    result = Matrix.zeros(self.rows, other.cols)

    for i in range(self.rows):
        for j in range(other.cols):
            total = 0
            for k in range(self.cols):
                total += self[i,k] * other[k,j]
            result[i,j] = total

    return result

Matrix.__matmul__ = matrix_multiply  # Use @ operator

# Example
A = Matrix([
    [1, 2],
    [3, 4]
])
B = Matrix([
    [5, 6],
    [7, 8]
])
C = A @ B
print(f"A @ B =\n{C}")
# [[1*5+2*7, 1*6+2*8], [3*5+4*7, 3*6+4*8]]
# [[19, 22], [43, 50]]

Matrix-Vector Multiplication

def matvec(self, v):
    """Multiply matrix by vector"""
    if self.cols != len(v):
        raise ValueError(f"Cannot multiply {self.shape} matrix by {len(v)}-vector")

    result = []
    for i in range(self.rows):
        total = sum(self[i,j] * v[j] for j in range(self.cols))
        result.append(total)

    return Vector(result)

Matrix.matvec = matvec

A = Matrix([
    [1, 2],
    [3, 4]
])
v = Vector([5, 6])
print(f"A × v = {A.matvec(v)}")  # Vector([17, 39])

3. Special Matrices

Symmetric Matrices

A matrix A is symmetric if A = A^T.

def is_symmetric(self, tolerance=1e-9):
    """Check if matrix is symmetric"""
    if self.rows != self.cols:
        return False
    for i in range(self.rows):
        for j in range(i+1, self.cols):
            if abs(self[i,j] - self[j,i]) > tolerance:
                return False
    return True

Matrix.is_symmetric = is_symmetric

Diagonal Matrices

Only diagonal elements are non-zero.

@staticmethod
def diagonal(values):
    """Create diagonal matrix from values"""
    n = len(values)
    data = [[values[i] if i == j else 0 for j in range(n)] for i in range(n)]
    return Matrix(data)

Matrix.diagonal = diagonal

D = Matrix.diagonal([1, 2, 3])
print(f"Diagonal matrix:\n{D}")

Sparse Matrices

Most elements are zero - use special representation.

class SparseMatrix:
    """Sparse matrix using dictionary of keys"""

    def __init__(self, rows, cols):
        self.rows = rows
        self.cols = cols
        self.data = {}  # (i,j) -> value

    def __getitem__(self, idx):
        return self.data.get(idx, 0)

    def __setitem__(self, idx, value):
        if value != 0:
            self.data[idx] = value
        elif idx in self.data:
            del self.data[idx]

    def to_dense(self):
        """Convert to dense matrix"""
        data = [[self[i,j] for j in range(self.cols)]
                for i in range(self.rows)]
        return Matrix(data)

    def density(self):
        """Fraction of non-zero elements"""
        return len(self.data) / (self.rows * self.cols)

4. Linear Transformations

Matrices represent linear transformations.

2D Transformations

def rotation_matrix_2d(theta):
    """2D rotation matrix by angle theta (radians)"""
    c, s = math.cos(theta), math.sin(theta)
    return Matrix([
        [c, -s],
        [s, c]
    ])

def scaling_matrix_2d(sx, sy):
    """2D scaling matrix"""
    return Matrix([
        [sx, 0],
        [0, sy]
    ])

def reflection_matrix_x():
    """Reflect across x-axis"""
    return Matrix([
        [1, 0],
        [0, -1]
    ])

def shear_matrix_x(k):
    """Shear in x direction"""
    return Matrix([
        [1, k],
        [0, 1]
    ])

# Example: Rotate point by 90 degrees
v = Vector([1, 0])
R = rotation_matrix_2d(math.pi / 2)  # 90 degrees
rotated = R.matvec(v)
print(f"Rotated: {rotated}")  # Approximately Vector([0, 1])

3D Transformations

def rotation_matrix_z(theta):
    """3D rotation around z-axis"""
    c, s = math.cos(theta), math.sin(theta)
    return Matrix([
        [c, -s, 0],
        [s, c, 0],
        [0, 0, 1]
    ])

def rotation_matrix_x(theta):
    """3D rotation around x-axis"""
    c, s = math.cos(theta), math.sin(theta)
    return Matrix([
        [1, 0, 0],
        [0, c, -s],
        [0, s, c]
    ])

def rotation_matrix_y(theta):
    """3D rotation around y-axis"""
    c, s = math.cos(theta), math.sin(theta)
    return Matrix([
        [c, 0, s],
        [0, 1, 0],
        [-s, 0, c]
    ])

Homogeneous Coordinates

Include translation with matrix multiplication.

def translation_matrix_2d(tx, ty):
    """2D translation using homogeneous coordinates"""
    return Matrix([
        [1, 0, tx],
        [0, 1, ty],
        [0, 0, 1]
    ])

def transform_point_2d(M, point):
    """Apply 3x3 transformation to 2D point"""
    homogeneous = Vector([point[0], point[1], 1])
    result = M.matvec(homogeneous)
    return Vector([result[0] / result[2], result[1] / result[2]])

# Translate (1, 2) by (3, 4)
T = translation_matrix_2d(3, 4)
point = Vector([1, 2])
translated = transform_point_2d(T, point)
print(f"Translated: {translated}")  # Vector([4, 6])

5. Efficient Matrix Operations

Strassen's Algorithm

Matrix multiplication in O(n^2.807) instead of O(n³).

def strassen_multiply(A, B):
    """
    Strassen's algorithm for matrix multiplication

    Reduces multiplications from 8 to 7 for 2x2 blocks
    Complexity: O(n^2.807) vs O(n³) for standard

    For simplicity, assumes square matrices of power-of-2 size
    """
    n = A.rows

    # Base case
    if n == 1:
        return Matrix([[A[0,0] * B[0,0]]])

    # Split into quadrants
    mid = n // 2

    def submatrix(M, row_start, col_start, size):
        data = [[M[row_start+i, col_start+j] for j in range(size)]
                for i in range(size)]
        return Matrix(data)

    def add_matrices(X, Y):
        return Matrix([[X[i,j] + Y[i,j] for j in range(X.cols)]
                       for i in range(X.rows)])

    def sub_matrices(X, Y):
        return Matrix([[X[i,j] - Y[i,j] for j in range(X.cols)]
                       for i in range(X.rows)])

    A11 = submatrix(A, 0, 0, mid)
    A12 = submatrix(A, 0, mid, mid)
    A21 = submatrix(A, mid, 0, mid)
    A22 = submatrix(A, mid, mid, mid)

    B11 = submatrix(B, 0, 0, mid)
    B12 = submatrix(B, 0, mid, mid)
    B21 = submatrix(B, mid, 0, mid)
    B22 = submatrix(B, mid, mid, mid)

    # Strassen's 7 products
    M1 = strassen_multiply(add_matrices(A11, A22), add_matrices(B11, B22))
    M2 = strassen_multiply(add_matrices(A21, A22), B11)
    M3 = strassen_multiply(A11, sub_matrices(B12, B22))
    M4 = strassen_multiply(A22, sub_matrices(B21, B11))
    M5 = strassen_multiply(add_matrices(A11, A12), B22)
    M6 = strassen_multiply(sub_matrices(A21, A11), add_matrices(B11, B12))
    M7 = strassen_multiply(sub_matrices(A12, A22), add_matrices(B21, B22))

    # Combine results
    C11 = add_matrices(sub_matrices(add_matrices(M1, M4), M5), M7)
    C12 = add_matrices(M3, M5)
    C21 = add_matrices(M2, M4)
    C22 = add_matrices(sub_matrices(add_matrices(M1, M3), M2), M6)

    # Assemble result
    result = Matrix.zeros(n, n)
    for i in range(mid):
        for j in range(mid):
            result[i, j] = C11[i, j]
            result[i, j+mid] = C12[i, j]
            result[i+mid, j] = C21[i, j]
            result[i+mid, j+mid] = C22[i, j]

    return result

Exercises

Basic

1. Implement vector operations: addition, scalar multiplication, dot product, magnitude.

2. Write a function to check if two vectors are orthogonal.

3. Implement matrix transpose and verify (AB)^T = B^T A^T.

Intermediate

4. Implement a function to check if a matrix is orthogonal (A^T A = I).

5. Create a 2D transformation that rotates around an arbitrary point (not origin).

6. Implement matrix exponentiation using repeated squaring.

Advanced

7. Implement blocked matrix multiplication for better cache efficiency.

8. Create a sparse matrix class with efficient multiplication.

9. Implement the Coppersmith-Winograd style multiplication (conceptually).

Summary

• Vectors are ordered lists of numbers with operations: add, scale, dot product, norm
• Dot product relates to angles: v·w = ||v|| ||w|| cos(θ)
• Matrices are 2D arrays representing linear transformations
• Matrix multiplication combines transformations (not commutative!)
• Homogeneous coordinates allow translation via matrix multiplication
• Special matrices: identity, diagonal, symmetric, sparse

Next Reading

Linear Systems and Gaussian Elimination →