⬡ Hub
Skip to content

SciPy: Linear Algebra

While NumPy provides excellent basic linear algebra capabilities, SciPy's scipy.linalg module extends this with more advanced and specialized routines. It includes functions for matrix decompositions, eigenvalues, matrix exponentials, and more, often built on highly optimized Fortran libraries like BLAS and LAPACK.

1. Basic Operations (Overlap with NumPy)

SciPy's linalg module includes most of the functions found in numpy.linalg, and often uses the same underlying optimized libraries.

import numpy as np
from scipy import linalg

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

# Determinant
print("Determinant of A:", linalg.det(A))

# Inverse
try:
    inv_A = linalg.inv(A)
    print("\nInverse of A:\n", inv_A)
    print("\nA * A_inv (should be identity):\n", A @ inv_A)
except linalg.LinAlgError as e:
    print(f"Cannot compute inverse: {e}")

# Solve linear system A @ x = B
# For A @ x = B where B is a vector or matrix
b_vec = np.array([1, 0])
x_vec = linalg.solve(A, b_vec)
print("\nSolution to A @ x = [1, 0]:", x_vec)
print("Check: A @ x_vec:", A @ x_vec)

2. Matrix Decompositions

Decompositions are fundamental to many linear algebra algorithms and are heavily used in numerical methods.

a. LU Decomposition (lu, lu_factor, lu_solve)

Decomposes a matrix A into a lower triangular matrix L, an upper triangular matrix U, and a permutation matrix P such that A = P L U. Used to solve linear systems efficiently.

import numpy as np
from scipy import linalg

A = np.array([[2, 5, 8],
              [3, 0, 1],
              [7, 9, 10]])
b = np.array([1, 2, 3])

# Perform LU decomposition
P, L, U = linalg.lu(A)
print("A:\n", A)
print("\nP:\n", P)
print("\nL:\n", L)
print("\nU:\n", U)
print("\nCheck P @ L @ U:\n", P @ L @ U) # Should be equal to A

# Solve A @ x = b using LU factors (more efficient for multiple 'b' vectors)
lu, piv = linalg.lu_factor(A) # lu_factor returns (LU matrix, pivots)
x = linalg.lu_solve((lu, piv), b)
print("\nSolution to A @ x = b using LU decomposition:", x)
print("Check A @ x:", A @ x)

b. Cholesky Decomposition (cholesky)

Decomposes a symmetric positive-definite matrix A into L L^T, where L is a lower triangular matrix. Often used in Monte Carlo simulations, Kalman filters, and to solve linear systems.

import numpy as np
from scipy import linalg

# A symmetric positive-definite matrix
A = np.array([[10, 2, 3],
              [2, 12, 4],
              [3, 4, 15]])

# Perform Cholesky decomposition
L = linalg.cholesky(A, lower=True)
print("Symmetric Positive-Definite Matrix A:\n", A)
print("\nCholesky Decomposition L:\n", L)
print("\nCheck L @ L.T:\n", L @ L.T) # Should be equal to A

c. QR Decomposition (qr)

Decomposes a matrix A into an orthogonal matrix Q and an upper triangular matrix R such that A = Q R. Used in solving linear least squares problems, eigenvalue problems, etc.

import numpy as np
from scipy import linalg

A = np.array([[1, 2],
              [3, 4],
              [5, 6]])

Q, R = linalg.qr(A)
print("Matrix A:\n", A)
print("\nQ (orthogonal matrix):\n", Q)
print("\nR (upper triangular matrix):\n", R)
print("\nCheck Q @ R:\n", Q @ R) # Should be equal to A
print("\nQ.T @ Q (should be identity):\n", Q.T @ Q) # Check orthogonality

3. Eigenvalue Problems (eig, eigh)

Finding eigenvalues and eigenvectors is crucial in principal component analysis (PCA), quantum mechanics, and many other fields.

  • eig(A): Computes eigenvalues and right eigenvectors for a general square matrix.
  • eigh(A): Computes eigenvalues and eigenvectors for a symmetric or Hermitian matrix (more stable and faster).
import numpy as np
from scipy import linalg

A = np.array([[2, 1],
              [1, 2]]) # A symmetric matrix

# Compute eigenvalues (w) and eigenvectors (v)
# w: eigenvalues, v: eigenvectors where v[:,i] is the eigenvector corresponding to w[i]
w, v = linalg.eigh(A) # Use eigh for symmetric matrices
print("Symmetric Matrix A:\n", A)
print("\nEigenvalues:", w)
print("Eigenvectors:\n", v)

# Check A @ v[:,0] = w[0] * v[:,0]
print("\nCheck A @ v[:,0]:", A @ v[:,0])
print("Check w[0] * v[:,0]:", w[0] * v[:,0])

4. Matrix Exponentials and Logarithms (expm, logm)

These functions are used in control theory, differential equations, and signal processing.

import numpy as np
from scipy import linalg

A = np.array([[1, 1],
              [0, 1]])

# Matrix exponential (e^A)
exp_A = linalg.expm(A)
print("Matrix A:\n", A)
print("\nMatrix Exponential (e^A):\n", exp_A)

# Matrix logarithm (log(A)) - inverse of expm
# Note: logm can return complex numbers for certain matrices
log_exp_A = linalg.logm(exp_A)
print("\nMatrix Logarithm of exp(A):\n", log_exp_A)

5. Other Advanced Functions:

  • Singular Value Decomposition (SVD): linalg.svd - fundamental for dimensionality reduction (like PCA) and solving linear least squares.
  • Pseudoinverse: linalg.pinv - computes the (Moore-Penrose) pseudo-inverse, useful for non-square matrices or singular matrices.
  • Norms: linalg.norm - computes various matrix or vector norms.

Further Topics:

  • Sparse linear algebra (scipy.sparse)
  • Iterative solvers for large linear systems
  • BLAS and LAPACK interfaces
  • Applications in signal processing and control theory.

SciPy's linalg module provides the robust and highly optimized tools necessary for tackling complex linear algebra problems, making it an indispensable part of any scientific computing toolkit in Python.