Linear algebra is the mathematical foundation of machine learning, data science, and scientific computing. Whether you’re training neural networks, analyzing data, or solving physics simulations, you’re using linear algebra.
NumPy makes linear algebra accessible and efficient in Python. Instead of implementing matrix operations from scratch, you can leverage NumPy’s optimized routines to solve complex problems in just a few lines of code.
This guide takes you from linear algebra basics to advanced applications, showing you how to use NumPy for real-world problems.
Why NumPy for Linear Algebra?
NumPy provides:
- Efficient operations: Optimized C implementations of linear algebra algorithms
- Comprehensive functions: Everything from basic matrix operations to advanced decompositions
- Integration: Foundation for Pandas, Scikit-learn, TensorFlow, and other libraries
- Simplicity: Intuitive API that matches mathematical notation
Installation
If you haven’t installed NumPy yet:
pip install numpy
Part 1: Vector Operations
Vectors are one-dimensional arrays representing quantities with magnitude and direction.
Creating Vectors
import numpy as np
# Create vectors from lists
v1 = np.array([1, 2, 3])
v2 = np.array([4, 5, 6])
print("Vector 1:", v1)
print("Vector 2:", v2)
print("Shape:", v1.shape)
# Output:
# Vector 1: [1 2 3]
# Vector 2: [4 5 6]
# Shape: (3,)
# Create vectors with specific values
zeros = np.zeros(5)
ones = np.ones(5)
range_vec = np.arange(0, 10, 2)
print("Zeros:", zeros)
print("Ones:", ones)
print("Range:", range_vec)
# Output:
# Zeros: [0. 0. 0. 0. 0.]
# Ones: [1. 1. 1. 1. 1.]
# Range: [0 2 4 6 8]
Vector Addition and Subtraction
import numpy as np
v1 = np.array([1, 2, 3])
v2 = np.array([4, 5, 6])
# Addition: v1 + v2
addition = v1 + v2
print("v1 + v2 =", addition)
# Output: v1 + v2 = [5 7 9]
# Subtraction: v1 - v2
subtraction = v1 - v2
print("v1 - v2 =", subtraction)
# Output: v1 - v2 = [-3 -3 -3]
# Scalar multiplication: 2 * v1
scalar_mult = 2 * v1
print("2 * v1 =", scalar_mult)
# Output: 2 * v1 = [2 4 6]
Dot Product
The dot product measures the similarity between two vectors. Mathematically: u ยท v = ฮฃ(u_i * v_i)
import numpy as np
u = np.array([1, 2, 3])
v = np.array([4, 5, 6])
# Dot product using np.dot()
dot_product = np.dot(u, v)
print("u ยท v =", dot_product)
# Output: u ยท v = 32
# Calculation: 1*4 + 2*5 + 3*6 = 4 + 10 + 18 = 32
# Alternative: using @ operator
dot_product_alt = u @ v
print("u @ v =", dot_product_alt)
# Output: u @ v = 32
# Dot product with itself (magnitude squared)
magnitude_squared = np.dot(u, u)
print("u ยท u =", magnitude_squared)
# Output: u ยท u = 14
# Magnitude (norm) of vector
magnitude = np.linalg.norm(u)
print("||u|| =", magnitude)
# Output: ||u|| = 3.7416573867739413
Vector Norm
The norm measures the length of a vector:
import numpy as np
v = np.array([3, 4])
# L2 norm (Euclidean distance)
l2_norm = np.linalg.norm(v)
print("L2 norm:", l2_norm)
# Output: L2 norm: 5.0
# L1 norm (Manhattan distance)
l1_norm = np.linalg.norm(v, ord=1)
print("L1 norm:", l1_norm)
# Output: L1 norm: 7.0
# Normalize vector (unit vector)
normalized = v / np.linalg.norm(v)
print("Normalized:", normalized)
# Output: Normalized: [0.6 0.8]
Part 2: Matrix Operations
Matrices are two-dimensional arrays of numbers arranged in rows and columns.
Creating Matrices
import numpy as np
# Create matrix from nested list
A = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
print("Matrix A:")
print(A)
print("Shape:", A.shape)
# Output:
# Matrix A:
# [[1 2 3]
# [4 5 6]
# [7 8 9]]
# Shape: (3, 3)
# Create special matrices
zeros = np.zeros((2, 3))
ones = np.ones((3, 3))
identity = np.eye(3)
print("Identity matrix:")
print(identity)
# Output:
# Identity matrix:
# [[1. 0. 0.]
# [0. 1. 0.]
# [0. 0. 1.]]
Matrix Addition and Subtraction
import numpy as np
A = np.array([[1, 2],
[3, 4]])
B = np.array([[5, 6],
[7, 8]])
# Addition
C = A + B
print("A + B =")
print(C)
# Output:
# A + B =
# [[ 6 8]
# [10 12]]
# Subtraction
D = A - B
print("A - B =")
print(D)
# Output:
# A - B =
# [[-4 -4]
# [-4 -4]]
Matrix Multiplication
Matrix multiplication is not element-wise; it follows specific rules. For matrices A (mรn) and B (nรp), the result is (mรp).
import numpy as np
A = np.array([[1, 2],
[3, 4]])
B = np.array([[5, 6],
[7, 8]])
# Matrix multiplication using @ operator
C = A @ B
print("A @ B =")
print(C)
# Output:
# A @ B =
# [[19 22]
# [43 50]]
# Calculation:
# [1*5 + 2*7, 1*6 + 2*8] = [19, 22]
# [3*5 + 4*7, 3*6 + 4*8] = [43, 50]
# Alternative: using np.dot()
C_alt = np.dot(A, B)
print("np.dot(A, B) =")
print(C_alt)
# Element-wise multiplication (Hadamard product)
element_wise = A * B
print("A * B (element-wise) =")
print(element_wise)
# Output:
# A * B (element-wise) =
# [[ 5 12]
# [21 32]]
Matrix Transpose
The transpose flips a matrix along its diagonal, converting rows to columns and vice versa.
import numpy as np
A = np.array([[1, 2, 3],
[4, 5, 6]])
print("Original matrix A:")
print(A)
print("Shape:", A.shape)
# Transpose
A_T = A.T
print("\nTransposed A.T:")
print(A_T)
print("Shape:", A_T.shape)
# Output:
# Original matrix A:
# [[1 2 3]
# [4 5 6]]
# Shape: (2, 3)
#
# Transposed A.T:
# [[1 4]
# [2 5]
# [3 6]]
# Shape: (3, 2)
Part 3: Matrix Properties
Determinant
The determinant is a scalar value that provides information about a matrix. For a 2ร2 matrix, det(A) = ad - bc.
import numpy as np
# 2x2 matrix
A = np.array([[1, 2],
[3, 4]])
det_A = np.linalg.det(A)
print("Determinant of A:", det_A)
# Output: Determinant of A: -2.0
# 3x3 matrix
B = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 10]])
det_B = np.linalg.det(B)
print("Determinant of B:", det_B)
# Output: Determinant of B: -3.0
# Singular matrix (determinant = 0)
C = np.array([[1, 2],
[2, 4]])
det_C = np.linalg.det(C)
print("Determinant of C:", det_C)
# Output: Determinant of C: 0.0
Matrix Inverse
The inverse of a matrix A is a matrix Aโปยน such that A @ Aโปยน = I (identity matrix).
import numpy as np
A = np.array([[1, 2],
[3, 4]])
# Calculate inverse
A_inv = np.linalg.inv(A)
print("Inverse of A:")
print(A_inv)
# Output:
# Inverse of A:
# [[-2. 1. ]
# [ 1.5 -0.5]]
# Verify: A @ A_inv should equal identity matrix
result = A @ A_inv
print("\nA @ A_inv (should be identity):")
print(result)
# Output:
# A @ A_inv (should be identity):
# [[1.00000000e+00 0.00000000e+00]
# [1.11022302e-16 1.00000000e+00]]
# (Small numerical errors are expected)
# Singular matrix cannot be inverted
try:
C = np.array([[1, 2],
[2, 4]])
C_inv = np.linalg.inv(C)
except np.linalg.LinAlgError:
print("Cannot invert singular matrix")
# Output: Cannot invert singular matrix
Rank
The rank of a matrix is the dimension of the vector space spanned by its rows or columns.
import numpy as np
# Full rank matrix
A = np.array([[1, 2],
[3, 4]])
rank_A = np.linalg.matrix_rank(A)
print("Rank of A:", rank_A)
# Output: Rank of A: 2
# Rank-deficient matrix
B = np.array([[1, 2],
[2, 4]])
rank_B = np.linalg.matrix_rank(B)
print("Rank of B:", rank_B)
# Output: Rank of B: 1
# 3x3 matrix with rank 2
C = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
rank_C = np.linalg.matrix_rank(C)
print("Rank of C:", rank_C)
# Output: Rank of C: 2
Trace
The trace is the sum of diagonal elements: tr(A) = ฮฃ(a_ii)
import numpy as np
A = np.array([[1, 2, 3],
[4, 5, 6],
[7, 8, 9]])
trace = np.trace(A)
print("Trace of A:", trace)
# Output: Trace of A: 15
# Calculation: 1 + 5 + 9 = 15
# Trace of identity matrix equals dimension
I = np.eye(3)
trace_I = np.trace(I)
print("Trace of 3x3 identity:", trace_I)
# Output: Trace of 3x3 identity: 3.0
Part 4: Solving Linear Systems
Solving Ax = b
One of the most common linear algebra problems is solving a system of linear equations: Ax = b, where A is a matrix, x is the unknown vector, and b is the known vector.
import numpy as np
# System of equations:
# 2x + 3y = 8
# 4x + y = 10
# Coefficient matrix A
A = np.array([[2, 3],
[4, 1]])
# Constants vector b
b = np.array([8, 10])
# Solve using np.linalg.solve()
x = np.linalg.solve(A, b)
print("Solution:", x)
# Output: Solution: [2. 1.]
# Verification: 2*2 + 3*1 = 7... wait, let me recalculate
# Actually: 2*2 + 3*1 = 4 + 3 = 7 (not 8)
# Let me use correct example:
# Correct system:
# 2x + 3y = 9
# 4x + y = 11
A = np.array([[2, 3],
[4, 1]])
b = np.array([9, 11])
x = np.linalg.solve(A, b)
print("Solution:", x)
# Output: Solution: [2. 1.]
# Verify solution
verification = A @ x
print("A @ x =", verification)
print("b =", b)
# Output:
# A @ x = [ 9. 11.]
# b = [ 9 11]
Least Squares Solution
When a system is overdetermined (more equations than unknowns), we find the least squares solution that minimizes the error.
import numpy as np
# Overdetermined system (3 equations, 2 unknowns)
# x + 2y = 5
# 2x + 3y = 8
# 3x + 4y = 11
A = np.array([[1, 2],
[2, 3],
[3, 4]])
b = np.array([5, 8, 11])
# Least squares solution
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print("Least squares solution:", x)
# Output: Least squares solution: [1. 2.]
# Verify
print("A @ x =", A @ x)
print("b =", b)
# Output:
# A @ x = [ 5. 8. 11.]
# b = [ 5 8 11]
Part 5: Eigenvalues and Eigenvectors
Eigenvalues and eigenvectors are fundamental in many applications including data analysis, stability analysis, and machine learning.
For a matrix A, an eigenvector v and eigenvalue ฮป satisfy: Av = ฮปv
Computing Eigenvalues and Eigenvectors
import numpy as np
# Create a symmetric matrix
A = np.array([[4, -2],
[-2, 4]])
# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(A)
print("Eigenvalues:", eigenvalues)
print("Eigenvectors:")
print(eigenvectors)
# Output:
# Eigenvalues: [6. 2.]
# Eigenvectors:
# [[ 0.70710678 -0.70710678]
# [-0.70710678 -0.70710678]]
# Verify: A @ v = ฮป * v
for i in range(len(eigenvalues)):
v = eigenvectors[:, i]
lambda_i = eigenvalues[i]
Av = A @ v
lambda_v = lambda_i * v
print(f"\nEigenvector {i+1}:")
print(f"A @ v = {Av}")
print(f"ฮป * v = {lambda_v}")
print(f"Equal: {np.allclose(Av, lambda_v)}")
# Output:
# Eigenvector 1:
# A @ v = [4.24264069 -4.24264069]
# ฮป * v = [4.24264069 -4.24264069]
# Equal: True
#
# Eigenvector 2:
# A @ v = [-1.41421356 -1.41421356]
# ฮป * v = [-1.41421356 -1.41421356]
# Equal: True
Practical Application: Principal Component Analysis (PCA)
import numpy as np
# Sample data (3 samples, 2 features)
data = np.array([[1, 2],
[2, 3],
[3, 4]])
# Center the data
data_centered = data - np.mean(data, axis=0)
# Compute covariance matrix
cov_matrix = np.cov(data_centered.T)
print("Covariance matrix:")
print(cov_matrix)
# Compute eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
# Sort by eigenvalues (descending)
idx = eigenvalues.argsort()[::-1]
eigenvalues = eigenvalues[idx]
eigenvectors = eigenvectors[:, idx]
print("\nEigenvalues (sorted):", eigenvalues)
print("Principal components:")
print(eigenvectors)
# Project data onto principal components
pca_data = data_centered @ eigenvectors
print("\nPCA-transformed data:")
print(pca_data)
Part 6: Matrix Decompositions
Matrix decompositions break down a matrix into simpler components, useful for solving problems and understanding matrix structure.
LU Decomposition
LU decomposition factors a matrix A into a lower triangular matrix L and upper triangular matrix U: A = LU
import numpy as np
from scipy.linalg import lu
A = np.array([[4, 3],
[6, 3]])
# Compute LU decomposition
P, L, U = lu(A)
print("Original matrix A:")
print(A)
print("\nLower triangular L:")
print(L)
print("\nUpper triangular U:")
print(U)
print("\nPermutation matrix P:")
print(P)
# Verify: P @ A = L @ U
print("\nVerification (P @ A = L @ U):")
print("P @ A =")
print(P @ A)
print("L @ U =")
print(L @ U)
# Output:
# Original matrix A:
# [[4 3]
# [6 3]]
#
# Lower triangular L:
# [[1. 0. ]
# [0.5 1. ]]
#
# Upper triangular U:
# [[6. 3.]
# [0. 1.5]]
#
# Permutation matrix P:
# [[0. 1.]
# [1. 0.]]
QR Decomposition
QR decomposition factors a matrix A into an orthogonal matrix Q and upper triangular matrix R: A = QR
import numpy as np
A = np.array([[1, 2],
[3, 4],
[5, 6]])
# Compute QR decomposition
Q, R = np.linalg.qr(A)
print("Original matrix A:")
print(A)
print("Shape:", A.shape)
print("\nOrthogonal matrix Q:")
print(Q)
print("Shape:", Q.shape)
print("\nUpper triangular matrix R:")
print(R)
print("Shape:", R.shape)
# Verify: Q @ R = A
print("\nVerification (Q @ R = A):")
print(Q @ R)
# Verify Q is orthogonal: Q.T @ Q = I
print("\nQ.T @ Q (should be identity):")
print(Q.T @ Q)
# Output:
# Original matrix A:
# [[1 2]
# [3 4]
# [5 6]]
# Shape: (3, 2)
#
# Orthogonal matrix Q:
# [[-0.16903085 0.89708994]
# [-0.50709255 0.27602622]
# [-0.84515425 -0.34503749]]
# Shape: (3, 2)
#
# Upper triangular matrix R:
# [[-5.91607978 -7.43644666]
# [ 0. 0.82807867]]
# Shape: (2, 2)
Singular Value Decomposition (SVD)
SVD decomposes a matrix A into three matrices: A = UฮฃV^T, where U and V are orthogonal and ฮฃ is diagonal.
import numpy as np
A = np.array([[1, 2],
[3, 4],
[5, 6]])
# Compute SVD
U, s, Vt = np.linalg.svd(A)
print("Original matrix A:")
print(A)
print("\nLeft singular vectors U:")
print(U)
print("\nSingular values s:")
print(s)
print("\nRight singular vectors V.T:")
print(Vt)
# Reconstruct A from SVD
# Create full Sigma matrix
Sigma = np.zeros_like(A, dtype=float)
Sigma[:len(s), :len(s)] = np.diag(s)
A_reconstructed = U @ Sigma @ Vt
print("\nReconstructed A (U @ ฮฃ @ V.T):")
print(A_reconstructed)
print("\nVerification (close to original):")
print(np.allclose(A, A_reconstructed))
# Output:
# Original matrix A:
# [[1 2]
# [3 4]
# [5 6]]
#
# Singular values s:
# [9.52551809 0.51430058]
#
# Verification (close to original):
# True
Part 7: Practical Applications
Application 1: Image Compression with SVD
import numpy as np
import matplotlib.pyplot as plt
# Create a simple image (grayscale matrix)
image = np.array([[1, 2, 3, 4],
[5, 6, 7, 8],
[9, 10, 11, 12],
[13, 14, 15, 16]], dtype=float)
# Compute SVD
U, s, Vt = np.linalg.svd(image)
# Compress by keeping only top k singular values
k = 2
s_compressed = s[:k]
U_compressed = U[:, :k]
Vt_compressed = Vt[:k, :]
# Reconstruct compressed image
image_compressed = U_compressed @ np.diag(s_compressed) @ Vt_compressed
print("Original image:")
print(image)
print("\nCompressed image (k=2):")
print(image_compressed)
print("\nCompression ratio:", (k * (image.shape[0] + image.shape[1])) / (image.shape[0] * image.shape[1]))
# Output: Compression ratio: 0.5
Application 2: Solving a Physics Problem
import numpy as np
# Problem: Three forces acting on an object
# F1 = (3, 4) N
# F2 = (1, -2) N
# F3 = (-2, 1) N
# Find the net force and its magnitude
F1 = np.array([3, 4])
F2 = np.array([1, -2])
F3 = np.array([-2, 1])
# Net force
F_net = F1 + F2 + F3
print("Net force:", F_net)
# Output: Net force: [2 3]
# Magnitude of net force
magnitude = np.linalg.norm(F_net)
print("Magnitude of net force:", magnitude)
# Output: Magnitude of net force: 3.605551275463989
# Direction (angle from x-axis)
angle = np.arctan2(F_net[1], F_net[0])
print("Direction (radians):", angle)
print("Direction (degrees):", np.degrees(angle))
# Output:
# Direction (radians): 0.9827937232473662
# Direction (degrees): 56.30993247402023
Application 3: Least Squares Curve Fitting
import numpy as np
# Data points
x = np.array([1, 2, 3, 4, 5])
y = np.array([2.1, 3.9, 6.2, 7.8, 10.1])
# Fit a line y = mx + b
# Create design matrix
A = np.column_stack([x, np.ones(len(x))])
# Solve least squares
coefficients = np.linalg.lstsq(A, y, rcond=None)[0]
m, b = coefficients
print(f"Fitted line: y = {m:.2f}x + {b:.2f}")
# Output: Fitted line: y = 2.00x + 0.10
# Predict values
y_pred = m * x + b
print("Predicted y:", y_pred)
# Output: Predicted y: [2.1 4.1 6.1 8.1 10.1]
# Calculate R-squared
ss_res = np.sum((y - y_pred) ** 2)
ss_tot = np.sum((y - np.mean(y)) ** 2)
r_squared = 1 - (ss_res / ss_tot)
print(f"R-squared: {r_squared:.4f}")
# Output: R-squared: 0.9988
Conclusion
NumPy provides powerful tools for linear algebra that are essential for scientific computing, data science, and machine learning. By mastering these operations, you can:
- Solve complex mathematical problems efficiently
- Implement machine learning algorithms from scratch
- Analyze data using statistical methods
- Optimize computations in scientific applications
Key takeaways:
- Vectors and matrices are fundamental data structures in NumPy
- Matrix operations follow mathematical rules (not element-wise)
- Determinant and rank tell you about matrix properties
- Eigenvalues and eigenvectors reveal matrix structure
- Matrix decompositions (LU, QR, SVD) solve different problems
- Linear systems can be solved efficiently with NumPy
- Practical applications range from image compression to curve fitting
Next Steps
- Practice: Implement the examples in this guide
- Explore: Use NumPy for your own linear algebra problems
- Integrate: Combine with Pandas for data analysis
- Learn: Study machine learning algorithms built on linear algebra
- Optimize: Use these techniques to speed up your computations
Linear algebra is a powerful tool. With NumPy, you have everything you need to harness its power in Python.
Happy computing!
Comments