NumPy is a fundamental Python library for scientific computing and data analysis. It provides efficient implementation of multidimensional arrays and matrices along with a vast collection of high-level mathematical functions to operate on these arrays. Linear algebra is one of the key strengths of NumPy and is widely used for data science, machine learning, computer vision, signal processing, and more. This guide will provide an in-depth overview of performing linear algebra operations in NumPy.
Table of Contents
Open Table of Contents
- Introduction
- Installing and Importing NumPy
- Creating Arrays and Matrices
- Array Attributes and Properties
- Array Indexing and Slicing
- Basic Matrix Operations
- The Dot Product
- Transposing and Reshaping Arrays
- Computing the Matrix Determinant
- Finding the Matrix Inverse
- Solving Linear Systems
- Eigenvalues and Eigenvectors
- Matrix Decomposition
- Solving Least Squares Problems
- Practical Examples and Uses
- Summary
Introduction
Linear algebra is a branch of mathematics focused on vector spaces and linear transformations between them. It involves mathematical concepts like matrices, determinants, eigenvalues, eigenvectors, vector products, and more. NumPy provides a n-dimensional array data structure called ndarray
along with vectorized implementations of most common linear algebra operations which allow us to efficiently work with large matrix datasets.
In this guide, we will cover the following core topics:
- Creating arrays and matrices in NumPy
- Basic properties and attributes of NumPy arrays
- Matrix math operations like addition, subtraction, multiplication etc.
- Computing the dot product between arrays
- Transposing and changing dimensions of arrays
- Computing matrix determinants and inverses
- Solving linear systems of equations
- Computing eigenvalues and eigenvectors
- Decomposing matrices through SVD, LU, Cholesky etc.
We will look at relevant example codes in NumPy to illustrate the usage of various linear algebra routines. But first, let’s start by setting up NumPy.
Installing and Importing NumPy
NumPy can be easily installed using pip:
pip install numpy
Import NumPy into any Python script as follows:
import numpy as np
We will use the standard alias np
for NumPy throughout this guide.
Creating Arrays and Matrices
The fundamental data structure in NumPy is the ndarray
class. We can create new ndarrays from lists or tuples in Python:
vector = np.array([1, 2, 3])
matrix = np.array([(1, 2, 3), (4, 5, 6)])
The ndarray
constructor converts nested sequences like lists of lists to multidimensional arrays.
We can also create specific matrix shapes with the zeros
and ones
functions:
# 3x3 matrix with all zeros
zero_matrix = np.zeros((3,3))
# 4x2 matrix with all ones
one_matrix = np.ones((4,2))
NumPy offers arange
and linspace
to generate numeric ranges:
# Array from 0 to 9
arange_array = np.arange(10)
# 10 numbers evenly spaced between 0 and 5
linspace_array = np.linspace(0, 5, 10)
The shape
attribute returns the array dimensions while dtype
returns the data type:
matrix = np.array([(1, 2, 3), (4, 5, 6)])
print(matrix.shape) # (2, 3)
print(matrix.dtype) # int64
Array Attributes and Properties
NumPy arrays have useful attributes that provide information about the data:
ndim
- Number of array dimensionssize
- Total number of elements in the arrayitemsize
- Storage size of each element in bytes
For example:
import numpy as np
arr = np.arange(10)
print(arr.ndim) # 1
print(arr.size) # 10
print(arr.itemsize) # 8 (64-bit system)
The actual data is stored in the data
buffer attribute as a multidimensional chunk of bytes.
We can also get basic statistical properties like the mean, standard deviation, variance etc. easily for NumPy arrays:
stats = np.random.normal(size=100)
print(stats.mean())
print(stats.std()) # standard deviation
print(stats.var()) # variance
Array Indexing and Slicing
NumPy arrays support familiar Python slicing syntax to extract subarrays:
arr = np.arange(10)
# First 5 elements
print(arr[:5])
# Last 3 elements
print(arr[-3:])
Specific elements can be accessed via integer indices:
matrix = np.array([(1, 2, 3), (4, 5, 6)])
# Get first row
print(matrix[0])
# Get element at row 1, column 2
print(matrix[1, 2]) # 6
Boolean index masks can also extract elements based on conditions:
rand_ints = np.random.randint(0, 10, 10)
# Elements greater than 5
mask = rand_ints > 5
print(rand_ints[mask])
Basic Matrix Operations
Let’s look at some basic mathematical operations between arrays and scalars.
Addition
arr = np.arange(5)
print(arr + 5) # Add 5 to each element
mat1 = np.ones((2, 3))
mat2 = np.ones((2, 3))
print(mat1 + mat2) # Element-wise addition
Subtraction
arr = np.arange(5)
print(arr - 2) # Subtract 2 from each element
Matrix subtraction is similarly element-wise.
Multiplication
Scalar multiplication:
arr = np.full((5,), 5)
print(arr * 10) # Multiply each element by 10
Matrix multiplication using dot()
:
mat1 = [[1, 2], [3, 4]]
mat2 = [[5, 6], [7, 8]]
mat1 = np.array(mat1)
mat2 = np.array(mat2)
print(np.dot(mat1, mat2))
Element-wise multiplication using *
:
mat1 = np.ones((2, 3))
mat2 = np.full((2, 3), 2)
print(mat1 * mat2)
Division
For scalar division:
arr = np.full((5,), 10)
print(arr / 2) # Divide each element by 2
Element-wise matrix division:
mat1 = np.full((3, 3), 6)
mat2 = np.full((3, 3), 3)
print(mat1 / mat2)
Exponentiation
Raising scalars to array powers:
import numpy as np
arr = np.array([1, 2, 3])
print(2**arr) # [2, 4, 8]
Element-wise exponentiation is done with **
operator:
mat = np.full((2, 2), 2)
print(mat**3)
The Dot Product
The dot product is a fundamental linear algebra operation that projects two vectors onto each other. In NumPy, we can compute dot products using the dot()
function:
import numpy as np
v1 = np.array([1, 2])
v2 = np.array([3, 4])
print(np.dot(v1, v2)) # 11
For matrices, dot()
performs standard matrix multiplication:
mat1 = [[1, 2],
[3, 4]]
mat2 = [[5, 6],
[7, 8]]
mat1 = np.array(mat1)
mat2 = np.array(mat2)
print(np.dot(mat1, mat2))
# [[19 22]
# [43 50]]
The dot product is associative so a.dot(b).dot(c)
is equivalent to a.dot(b.dot(c))
.
We can also compute dot products across specific axes using np.tensordot()
:
mat1 = np.arange(8).reshape(2, 2, 2)
mat2 = np.arange(8).reshape(2, 2, 2)
print(np.tensordot(mat1, mat2, axes=((1, 0), (0, 1))))
Here we took dot products over the given axes combinations.
Transposing and Reshaping Arrays
The transpose()
method returns a view of the array with axes swapped:
mat = np.random.rand(3, 5)
print(mat.transpose()) # Transposed view
To reshape an array without changing data, use reshape()
:
vector = np.arange(10)
print(vector.reshape(5, 2)) # Reshape to 5x2 matrix
We can also flatten to 1D using flatten()
or ravel to continuous memory using ravel()
.
To add a new axis, use np.expand_dims()
:
vector = np.arange(5)
col_vector = np.expand_dims(vector, axis=1)
# [[0],
# [1],
# [2],
# [3],
# [4]]
Computing the Matrix Determinant
The determinant is a scalar value that provides important information about a matrix like its invertibility. In NumPy, we can find the determinant using np.linalg.det()
:
mat = np.array([[1, 2],
[3, 4]])
print(np.linalg.det(mat)) # -2.0
For non-square matrices, det()
raises an exception. The determinant value is useful for matrix inversion and computing eigenvalues.
Finding the Matrix Inverse
The inverse of a square matrix A
is denoted as A^-1
and satisfies:
A * A^-1 = I
Where I
is the identity matrix.
In NumPy, we can compute the inverse as follows:
from numpy.linalg import inv
mat = np.array([[1, 2],
[3, 4]])
mat_inv = inv(mat)
print(mat_inv)
# [[-2. , 1. ],
# [ 1.5, -0.5]]
This finds the unique inverse for the matrix. An exception is raised if the inverse does not exist.
Solving Linear Systems
We can use matrix inverses to solve systems of linear equations.
For example, consider:
3x + 5y = 14
2x + 4y = 8
This can be represented in matrix form as:
A * x = b
Where:
A = np.array([[3, 5], [2, 4]])
x = np.array([x, y])
b = np.array([14, 8])
Now to solve for x:
import numpy as np
A = np.array([[3, 5], [2, 4]])
b = np.array([14, 8])
x = np.linalg.inv(A).dot(b)
print(x)
# [2. 2.] # x = 2, y = 2
Here we computed the inverse of the coefficients matrix A
and multiplied it by b
to find x
.
Eigenvalues and Eigenvectors
Eigenvalues and eigenvectors are useful for matrix decomposition, principal component analysis, and more.
The eigenvalues can be found using np.linalg.eigvals()
:
mat = np.array([[1, 2],
[3, 4]])
print(np.linalg.eigvals(mat))
# [5.37228404 -0.37228404]
It returns the eigenvalues sorted in ascending order.
For eigenvectors, use np.linalg.eig()
:
evals, evecs = np.linalg.eig(mat)
print(evals) # Eigenvalues
print(evecs) # Matrix with eigenvectors as columns
The columns of the eigenvector matrix correspond to the eigenvectors for the sorted eigenvalues.
Matrix Decomposition
NumPy provides various matrix decomposition methods that are useful for machine learning and data analysis.
LU Decomposition
LU decomposition splits a matrix into lower and upper triangular matrices. We can perform LU decomposition in NumPy using np.linalg.lu()
:
A = np.array([[1, 2], [3, 4]])
P, L, U = np.linalg.lu(A)
print(L)
# [[1. 0.]
# [3. 1.]]
print(U)
# [[1. 2.]
# [0. 4.]]
This also returns a permutation matrix P
to apply row swaps.
Cholesky Decomposition
Cholesky decomposition is used for positive definite matrices. The np.linalg.cholesky()
method returns the Cholesky lower triangular factor:
A = np.array([[4, -2, 1], [-2, 6, -3], [1, -3, 7]])
L = np.linalg.cholesky(A)
print(L)
# [[2. 0. 0.]
# [-1. 1. 0.]
# [ 1. -1. 1.]]
Where L * L.T = A
.
Singular Value Decomposition (SVD)
SVD factorizes a matrix into singular vectors and singular values.
In NumPy:
A = np.array([[1, 0, 0, 0, 2],
[0, 0, 3, 0, 0],
[0, 0, 0, 0, 0],
[0, 2, 0, 0, 0]])
U, s, Vh = np.linalg.svd(A)
print(s)
# [4. 2. 1. 0.]
s
contains the singular values sorted by size. U
and Vh
contain the left and right singular vectors for each singular value.
SVD has many applications in signal processing, statistics, and machine learning.
Solving Least Squares Problems
We can use NumPy to solve least squares problems with linear equality constraints. The np.linalg.lstsq()
method computes the vector x
that minimizes the equation:
||Ax - b||_2
For example:
A = np.array([[1, 1], [2, 1], [3, 1]])
b = np.array([5, 7, 9])
x, residuals, rank, s = np.linalg.lstsq(A, b)
print(x)
# [2. 1.]
print(residuals)
# 1.1102230246251565e-16
# (near zero)
This finds the optimal x
while minimizing the residuals for the overdetermined system.
The return value contains the residuals sum, matrix rank, and singular values.
Practical Examples and Uses
Let’s go through some practical linear algebra examples and applications in data science and machine learning.
Principal Component Analysis (PCA)
PCA is used for dimensionality reduction and feature extraction. We’ll use SVD to perform PCA on a sample dataset:
from sklearn.datasets import load_iris
X, y = load_iris(return_X_y=True)
U, s, Vh = np.linalg.svd(X)
# Components in U with largest singular values
pca = U[:, :2]
print(pca.shape) # (150, 2)
This derives the two primary principal components from the Iris flower dataset. SVR reduced it from 4D to 2D while preserving the main latent features.
Pseudoinverse for Least Squares
The pseudoinverse found via SVD can solve least squares problems:
A = np.array([[1, 1], [3, 2], [4, 1]])
b = np.array([5, 6, 7])
A_pinv = np.linalg.pinv(A)
x = A_pinv.dot(b)
print(x)
# [2. 1.]
So SVD provides an efficient way to compute the pseudoinverse.
Image Compression with SVD
We can use SVD for image compression by reconstructing using only top singular values and vectors:
import matplotlib.pyplot as plt
from skimage import io
img = io.imread('image.jpg')
U, s, Vh = np.linalg.svd(img)
# Reconstruct using top 10 singular values
s10 = np.zeros(s.shape)
s10[:10] = s[:10]
img10 = U.dot(np.diag(s10).dot(Vh))
plt.imshow(img10, cmap='gray')
plt.show()
This provides an optimized lower rank approximation for the image.
Summary
To summarize, NumPy provides efficient Pythonic implementations of:
- Array creation routines like
array()
,zeros()
,ones()
,linspace()
etc. - Array attributes like
shape
,size
,ndim
and methods likereshape()
,flatten()
etc. - Mathematical operations like addition, subtraction, multiplication,