This article covers orthogonality and projections, part eight of the series. Orthogonality is fundamental to many algorithms including least squares regression, QR decomposition, and machine learning techniques like PCA.

Linear Algebra: Python Series - View all articles in this series.

Previous articles in this series:

1. Linear Algebra: Vectors
2. Linear Algebra: Matrices
3. Linear Algebra: Systems of Linear Equations
4. Linear Algebra: Matrix Inverses and Determinants
5. Linear Algebra: Vector Spaces and Subspaces
6. Linear Algebra: Eigenvalues and Eigenvectors Part 1
7. Linear Algebra: Eigenvalues and Eigenvectors Part 2

Previous articles covered vectors, matrices, systems, inverses, vector spaces, and eigenvalues, all linked above.

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

§2026 Update

The code here runs as written, and np.linalg.qr is still the workhorse. Two things worth knowing in practice.

The classical Gram-Schmidt shown here is the clearest way to watch orthogonalization happen, but it is numerically shaky: in floating point it gradually loses orthogonality as the vectors pile up. When you actually need an orthonormal basis, get it from a QR decomposition. np.linalg.qr uses Householder reflections under the hood and stays orthogonal where hand-rolled Gram-Schmidt drifts.

The projection matrix here is written as $\boldsymbol{A}(\boldsymbol{A}^T\boldsymbol{A})^{-1}\boldsymbol{A}^T$, which is the clean way to state it. For real computation, skip forming $(\boldsymbol{A}^T\boldsymbol{A})^{-1}$, which squares the conditioning, and project with the $\boldsymbol{Q}$ from a QR instead: $\boldsymbol{Q}\boldsymbol{Q}^T\vec{b}$. Same projection, computed more stably. The next article puts all of this to work in least squares.

The article continues below. The 2026 Update above covers what’s worth knowing today.

§Orthogonal Vectors

Two vectors $\vec{u}$ and $\vec{v}$ are orthogonal (perpendicular) if their dot product is zero:

# Orthogonal vectors
u = np.array([1, 0, 1])
v = np.array([0, 1, 0])

dot_product = np.dot(u, v)
print(f"u = {u}")
print(f"v = {v}")
print(f"u . v = {dot_product}")
print(f"Orthogonal: {dot_product == 0}")

u = [1 0 1]
v = [0 1 0]
u . v = 0
Orthogonal: True

§Vector Projection

The projection of vector $\vec{b}$ onto vector $\vec{a}$ is:

def vector_projection(b, a):
    """Project vector b onto vector a."""
    return (np.dot(a, b) / np.dot(a, a)) * a

a = np.array([2, 0])
b = np.array([1, 2])

proj = vector_projection(b, a)
residual = b - proj

print(f"a = {a}")
print(f"b = {b}")
print(f"proj_a(b) = {proj}")
print(f"residual = b - proj = {residual}")
print(f"proj . residual = {np.dot(proj, residual):.10f} (should be 0)")

a = [2 0]
b = [1 2]
proj_a(b) = [1. 0.]
residual = b - proj = [0. 2.]
proj . residual = 0.0000000000 (should be 0)

# Visualize the projection
plt.figure(figsize=(8, 6))
origin = np.array([0, 0])

plt.quiver(*origin, *a, color='blue', scale=1, scale_units='xy', angles='xy', label='a')
plt.quiver(*origin, *b, color='green', scale=1, scale_units='xy', angles='xy', label='b')
plt.quiver(*origin, *proj, color='red', scale=1, scale_units='xy', angles='xy', label='proj_a(b)')
plt.plot([b[0], proj[0]], [b[1], proj[1]], 'k--', label='residual')

plt.xlim(-1, 3)
plt.ylim(-1, 3)
plt.grid(True, alpha=0.3)
plt.legend()
plt.axis('equal')
plt.title('Vector Projection')
plt.show()

§The Gram-Schmidt Process

The Gram-Schmidt process transforms a set of linearly independent vectors into an orthonormal set (orthogonal and unit length).

Given vectors ${\vec{v}_1, \vec{v}_2, \ldots, \vec{v}_n}$:

Then normalize: $\vec{e}_k = \frac{\vec{u}_k}{|\vec{u}_k|}$

def gram_schmidt(vectors):
    """Apply Gram-Schmidt orthonormalization."""
    vectors = np.array(vectors, dtype=float)
    n = len(vectors)
    orthonormal = []

    for i in range(n):
        v = vectors[i].copy()

        # Subtract projections onto previous vectors
        for u in orthonormal:
            v -= np.dot(v, u) * u

        # Normalize
        norm = np.linalg.norm(v)
        if norm > 1e-10:  # Check for linear dependence
            orthonormal.append(v / norm)

    return np.array(orthonormal)

# Apply Gram-Schmidt
v1 = np.array([1, 1, 0])
v2 = np.array([1, 0, 1])
v3 = np.array([0, 1, 1])

orthonormal = gram_schmidt([v1, v2, v3])

print("Original vectors:")
print(f"v1 = {v1}")
print(f"v2 = {v2}")
print(f"v3 = {v3}")

print("\nOrthonormal vectors:")
for i, e in enumerate(orthonormal):
    print(f"e{i+1} = {e}")

Original vectors:
v1 = [1 1 0]
v2 = [1 0 1]
v3 = [0 1 1]

Orthonormal vectors:
e1 = [0.70710678 0.70710678 0.        ]
e2 = [ 0.40824829 -0.40824829  0.81649658]
e3 = [-0.57735027  0.57735027  0.57735027]

# Verify orthonormality
print("Verification:")
print("Dot products (should be 0 for i != j):")
for i in range(3):
    for j in range(i, 3):
        dot = np.dot(orthonormal[i], orthonormal[j])
        print(f"  e{i+1} . e{j+1} = {dot:.6f}")

Verification:
Dot products (should be 0 for i != j):
  e1 . e1 = 1.000000
  e1 . e2 = 0.000000
  e1 . e3 = 0.000000
  e2 . e2 = 1.000000
  e2 . e3 = -0.000000
  e3 . e3 = 1.000000

§QR Decomposition

QR decomposition factors a matrix into an orthogonal matrix $\boldsymbol{Q}$ and an upper triangular matrix $\boldsymbol{R}$:

This is essentially Gram-Schmidt applied to the columns of $\boldsymbol{A}$.

# Create a matrix
A = np.array([[1, 1, 0],
              [1, 0, 1],
              [0, 1, 1]], dtype=float)

# QR decomposition
Q, R = np.linalg.qr(A)

print(f"Matrix A:\n{A}\n")
print(f"Orthogonal matrix Q:\n{Q}\n")
print(f"Upper triangular R:\n{R}\n")

# Verify
print(f"Q @ R:\n{Q @ R}")
print(f"\nEquals A: {np.allclose(A, Q @ R)}")

Matrix A:
[[1. 1. 0.]
 [1. 0. 1.]
 [0. 1. 1.]]

Orthogonal matrix Q:
[[-0.70710678  0.40824829 -0.57735027]
 [-0.70710678 -0.40824829  0.57735027]
 [-0.          0.81649658  0.57735027]]

Upper triangular R:
[[-1.41421356 -0.70710678 -0.70710678]
 [ 0.          1.22474487  0.40824829]
 [ 0.          0.          1.15470054]]

Q @ R:
[[1. 1. 0.]
 [1. 0. 1.]
 [0. 1. 1.]]

Equals A: True

# Verify Q is orthogonal: Q^T @ Q = I
print("Q^T @ Q (should be identity):")
print(Q.T @ Q)
print(f"\nIs orthogonal: {np.allclose(Q.T @ Q, np.eye(3))}")

Q^T @ Q (should be identity):
[[ 1.00000000e+00  1.02766874e-16 -1.28818495e-16]
 [ 1.02766874e-16  1.00000000e+00 -1.03593182e-16]
 [-1.28818495e-16 -1.03593182e-16  1.00000000e+00]]

Is orthogonal: True

§Projection onto a Subspace

To project a vector $\vec{b}$ onto a subspace spanned by columns of $\boldsymbol{A}$:

The matrix $\boldsymbol{P} = \boldsymbol{A}(\boldsymbol{A}^T\boldsymbol{A})^{-1}\boldsymbol{A}^T$ is called the projection matrix.

def projection_matrix(A):
    """Compute the projection matrix onto the column space of A."""
    return A @ np.linalg.inv(A.T @ A) @ A.T

# Project onto a plane in R^3
# Plane spanned by [1,0,1] and [0,1,1]
A = np.array([[1, 0],
              [0, 1],
              [1, 1]])

P = projection_matrix(A)
print(f"Projection matrix P:\n{P}\n")

# Project a point onto the plane
b = np.array([1, 2, 5])
proj = P @ b

print(f"Original vector b = {b}")
print(f"Projection onto plane: {proj}")
print(f"Residual: {b - proj}")

Projection matrix P:
[[ 0.66666667 -0.33333333  0.33333333]
 [-0.33333333  0.66666667  0.33333333]
 [ 0.33333333  0.33333333  0.66666667]]

Original vector b = [1 2 5]
Projection onto plane: [1.66666667 2.66666667 4.33333333]
Residual: [-0.66666667 -0.66666667  0.66666667]

§Properties of Projection Matrices

Projection matrices satisfy:

1. $\boldsymbol{P}^2 = \boldsymbol{P}$ (idempotent)
2. $\boldsymbol{P}^T = \boldsymbol{P}$ (symmetric)

# Verify properties
print("Projection matrix properties:")
print(f"P^2 = P: {np.allclose(P @ P, P)}")
print(f"P^T = P: {np.allclose(P.T, P)}")

Projection matrix properties:
P^2 = P: True
P^T = P: True

§Orthogonal Matrices

An orthogonal matrix $\boldsymbol{Q}$ satisfies:

Properties:

• $\boldsymbol{Q}^{-1} = \boldsymbol{Q}^T$
• Columns (and rows) are orthonormal
• Preserves lengths and angles
• $|\det(\boldsymbol{Q})| = 1$

# Rotation matrix is orthogonal
theta = np.pi / 4  # 45 degrees
R = np.array([[np.cos(theta), -np.sin(theta)],
              [np.sin(theta),  np.cos(theta)]])

print(f"45-degree rotation matrix:\n{R}\n")
print(f"R^T @ R:\n{R.T @ R}\n")
print(f"Is orthogonal: {np.allclose(R.T @ R, np.eye(2))}")
print(f"Determinant: {np.linalg.det(R):.6f}")

45-degree rotation matrix:
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]

R^T @ R:
[[1. 0.]
 [0. 1.]]

Is orthogonal: True
Determinant: 1.000000

§Summary

In this article, we covered:

• Orthogonal vectors: $\vec{u} \cdot \vec{v} = 0$
• Vector projection: $\text{proj}_{\vec{a}}\vec{b}$
• Gram-Schmidt process: Creating orthonormal bases
• QR decomposition: $\boldsymbol{A} = \boldsymbol{Q}\boldsymbol{R}$
• Projection matrices: $\boldsymbol{P} = \boldsymbol{A}(\boldsymbol{A}^T\boldsymbol{A})^{-1}\boldsymbol{A}^T$
• Orthogonal matrices: $\boldsymbol{Q}^T\boldsymbol{Q} = \boldsymbol{I}$

§Resources

• Orthogonality - Khan Academy
• Gram-Schmidt - Wikipedia
• NumPy QR

Linear Algebra: Python Series - View all articles in this series.