An In-Depth Guide to Linear Algebra in NumPy

Updated: at 04:28 AM

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

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:

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:

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: