Vectors and Matrices

Linear algebra is the mathematics of vectors and matrices. It's fundamental to computer graphics, machine learning, scientific computing, data compression, and quantum computing. ## 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

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

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

  3. Implement matrix exponentiation using repeated squaring.

Advanced

  1. Implement blocked matrix multiplication for better cache efficiency.

  2. Create a sparse matrix class with efficient multiplication.

  3. 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 →