SciPy - Linalg

SciPy - Linear Algebra (scipy.linalg)

The scipy.linalg module provides a comprehensive set of linear algebra operations. It is built upon optimized libraries like ATLAS, LAPACK (Linear Algebra PACKage), and BLAS (Basic Linear Algebra Subprograms), ensuring high performance for matrix computations. SciPy’s linear algebra routines expect NumPy arrays (or objects convertible to them) as input and typically return results as NumPy arrays.

scipy.linalg vs. numpy.linalg

scipy.linalg contains all functions available in numpy.linalg. However, scipy.linalg offers additional, more advanced functions. A key advantage is that scipy.linalg is always compiled with BLAS/LAPACK support, guaranteeing optimized performance. For NumPy, BLAS/LAPACK support is optional during installation, so scipy.linalg functions can sometimes be faster or more robust.

It’s generally recommended to use scipy.linalg for linear algebra tasks in SciPy-based projects.

Solving Linear Systems

A common task is solving a system of linear equations, represented as $A \mathbf{x} = \mathbf{b}$, where $A$ is a matrix of coefficients, $\mathbf{x}$ is the vector of unknowns, and $\mathbf{b}$ is the result vector.

The scipy.linalg.solve(a, b) function computes $\mathbf{x}$.

Consider the system:

• $3x_1 + 2x_2 = 2$
• $1x_1 - 1x_2 = 4$
• $5x_2 + 1x_3 = -1$

This can be written as $A \mathbf{x} = \mathbf{b}$ with:

$$ A = \begin{bmatrix} 3 & 2 & 0 \ 1 & -1 & 0 \ 0 & 5 & 1 \end{bmatrix}, \quad \mathbf{x} = \begin{bmatrix} x_1 \ x_2 \ x_3 \end{bmatrix}, \quad \mathbf{b} = \begin{bmatrix} 2 \ 4 \ -1 \end{bmatrix} $$

import numpy as np
from scipy import linalg

# Coefficient matrix A
A = np.array([
    [3, 2, 0],
    [1, -1, 0],
    [0, 5, 1]
])

# Result vector b
b_vec = np.array([2, 4, -1])

# Solve for x
x_solution = linalg.solve(A, b_vec)

print("Solution vector x:")
print(x_solution)

# Verify the solution: A @ x should be close to b
print("\nVerification (A @ x):")
print(A @ x_solution)
# Check if A @ x_solution is close to b_vec using np.allclose
print(f"Is solution correct? {np.allclose(A @ x_solution, b_vec)}")

Output:

Solution vector x:
[ 2. -2.  9.]

Verification (A @ x):
[ 2.  4. -1.]
Is solution correct? True

Finding the Determinant

The determinant is a scalar value associated with a square matrix, denoted as $\det(A)$ or $|A|$. It’s computed using scipy.linalg.det(a).

import numpy as np
from scipy import linalg

# Define a square matrix
M = np.array([
    [1, 2],
    [3, 4]
])

# Calculate the determinant
det_M = linalg.det(M)

print(f"Matrix M:\n{M}")
print(f"Determinant of M: {det_M}") # Expected: 1*4 - 2*3 = -2

Output:

Matrix M:
[[1 2]
 [3 4]]
Determinant of M: -2.0

Eigenvalues and Eigenvectors

The eigenvalue-eigenvector problem involves finding scalars $\lambda$ (eigenvalues) and corresponding non-zero vectors $\mathbf{v}$ (eigenvectors) for a square matrix $A$ such that $A\mathbf{v} = \lambda\mathbf{v}$.

scipy.linalg.eig(a) computes eigenvalues and eigenvectors. It returns a tuple (eigenvalues, eigenvectors), where eigenvectors[:, i] corresponds to eigenvalues[i].

import numpy as np
from scipy import linalg

# Define a square matrix
A_eig = np.array([
    [1, 5],
    [3, 4]
])

# Calculate eigenvalues and eigenvectors
eigenvalues, eigenvectors = linalg.eig(A_eig)

print(f"Matrix A:\n{A_eig}")
print(f"\nEigenvalues (lambda):\n{eigenvalues}")
print(f"\nEigenvectors (v columns):\n{eigenvectors}")

# Verification for the first eigenvalue/eigenvector pair
lambda1 = eigenvalues[0]
v1 = eigenvectors[:, 0]
print(f"\nVerification for first pair (A @ v1):\n{A_eig @ v1}")
print(f"Verification for first pair (lambda1 * v1):\n{lambda1 * v1}")
print(f"Are they close? {np.allclose(A_eig @ v1, lambda1 * v1)}")

Output (eigenvalues/vectors can be complex for non-symmetric matrices):

Matrix A:
[[1 5]
 [3 4]]

Eigenvalues (lambda):
[-1.40512484+0.j  6.40512484+0.j]

Eigenvectors (v columns):
[[-0.84797312 -0.65708331]
 [ 0.53004982 -0.75381815]]

Verification for first pair (A @ v1):
[ 1.19165751+0.j -0.74484869+0.j]
Verification for first pair (lambda1 * v1):
[ 1.19165751+0.j -0.74484869+0.j]
Are they close? True

Singular Value Decomposition (SVD)

Singular Value Decomposition is a factorization of a real or complex matrix $M$ into $M = U \Sigma V^H$, where:

• $U$ is a unitary matrix having left singular vectors as columns.
• $\Sigma$ is a rectangular diagonal matrix with non-negative real numbers (singular values) on the diagonal.
• $V^H$ (the conjugate transpose of $V$) is a unitary matrix having right singular vectors as columns.

SVD is a powerful tool used in dimensionality reduction, solving linear systems, and more. scipy.linalg.svd(a) computes this.

import numpy as np
from scipy import linalg

# Define a matrix (can be non-square)
# Using a random matrix for demonstration
np.random.seed(0) # for reproducibility
A_svd = np.random.randn(3, 2) + 1j * np.random.randn(3, 2) # A 3x2 complex matrix

# Perform SVD
U, s, Vh = linalg.svd(A_svd)

# s contains the singular values (1D array)
# To reconstruct Sigma, create a diagonal matrix of appropriate shape
Sigma = np.zeros(A_svd.shape, dtype=complex) # A_svd.shape is (3,2)
Sigma[:A_svd.shape[1], :A_svd.shape[1]] = np.diag(s) # s has min(3,2)=2 elements

# Reconstruct the original matrix
A_reconstructed = U @ Sigma @ Vh

print(f"Original Matrix A (shape {A_svd.shape}):\n{A_svd}")
print(f"\nU matrix (shape {U.shape}):\n{U}")
print(f"\nSingular values (s - 1D array, shape {s.shape}):\n{s}")
print(f"\nSigma matrix (shape {Sigma.shape}):\n{Sigma}")
print(f"\nVh matrix (conjugate transpose of V, shape {Vh.shape}):\n{Vh}")
print(f"\nReconstructed Matrix (U @ Sigma @ Vh):")
print(A_reconstructed)
print(f"\nIs reconstruction close to original? {np.allclose(A_svd, A_reconstructed)}")

The output will show the matrices U, Vh, the singular values s, the constructed Sigma, and the reconstructed A, which should be very close to the original A_svd.

scipy.linalg offers many other functions, including matrix factorizations (LU, QR, Cholesky), matrix exponentials, and special matrix constructions. Always consult the official SciPy linalg documentation for a full list and advanced usage.