Implementing Decompositions with SciPy/NumPy

The theory behind matrix decompositions like SVD, LU, and QR is essential, but performing these calculations by hand is impractical for matrices of any significant size. Fortunately, Python's scientific computing libraries, specifically NumPy and SciPy, provide efficient and numerically stable implementations of these algorithms.

We'll primarily use functions from scipy.linalg, which often offers more advanced features and optimizations compared to numpy.linalg, though NumPy provides foundational versions as well. Ensure you have NumPy and SciPy installed and import them:

import numpy as np
from scipy import linalg

Singular Value Decomposition (SVD)

SVD decomposes a matrix $A$ into $U \Sigma V^T$. The function scipy.linalg.svd handles this decomposition.

# Define a sample matrix A
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9],
              [10, 11, 12]])

print("Original Matrix A ({}):".format(A.shape))
print(A)

# Perform SVD
# full_matrices=False computes the economy SVD (more efficient for non-square)
U, s, Vh = linalg.svd(A, full_matrices=False)

print("\nU Matrix ({}):".format(U.shape))
print(U)

# s contains the singular values (1D array)
print("\nSingular Values (s):", s.shape)
print(s)

# Vh is the transpose of V (V.T)
print("\nVh Matrix (V transpose) ({}):".format(Vh.shape))
print(Vh)

Understanding the Output:

• U: The matrix of left singular vectors. Its shape depends on full_matrices. With full_matrices=False for an $m \times n$ matrix where $m \ge n$, $U$ will be $m \times n$.
• s: A 1D NumPy array containing the singular values ($\sigma_1, \sigma_2, ..., \sigma_k$) in descending order. It's not the diagonal matrix $\Sigma$.
• Vh: The matrix of right singular vectors, already transposed ($V^T$). Its shape depends on full_matrices. With full_matrices=False for $m \ge n$, $Vh$ will be $n \times n$.

Reconstructing the Original Matrix (Verification):

To verify the decomposition $A = U \Sigma V^T$, you need to construct the diagonal matrix $\Sigma$ from the singular values s.

# Construct the Sigma matrix (diagonal)
# Need to match dimensions for multiplication
Sigma = np.zeros(A.shape) # Start with zeros in the shape of A
# Populate the diagonal with singular values
# min(A.shape) gives the number of singular values
Sigma[:A.shape[1], :A.shape[1]] = np.diag(s) # Using shape of Vh (n x n) here

# Reconstruct A
A_reconstructed = U @ Sigma @ Vh # Using matrix multiplication operator @

print("\nSigma Matrix ({}):".format(Sigma.shape))
print(Sigma)

print("\nReconstructed Matrix (U @ Sigma @ Vh):")
print(A_reconstructed)

# Check if close to original (due to floating point precision)
print("\nIs reconstructed close to original?", np.allclose(A, A_reconstructed))

The np.allclose() function is useful for comparing floating-point arrays, accounting for small numerical inaccuracies. Using full_matrices=False (economy SVD) is generally recommended when $m \neq n$ as it avoids computing unnecessary zero vectors in $U$ or $V$.

LU Decomposition

LU decomposition factors a square matrix $A$ into $L$ (lower triangular) and $U$ (upper triangular), such that $A = LU$. Often, pivoting (row swapping) is needed for stability, resulting in $PA = LU$, where $P$ is a permutation matrix. scipy.linalg.lu handles this.

# Define a square matrix B
B = np.array([[2, 5, 8],
              [4, 2, 1],
              [9, 3, 7]])

print("Original Matrix B ({}):".format(B.shape))
print(B)

# Perform LU decomposition with permutation
P_mat, L, U = linalg.lu(B)

# Note: P_mat is the permutation matrix itself
print("\nPermutation Matrix P ({}):".format(P_mat.shape))
print(P_mat)

print("\nLower Triangular Matrix L ({}):".format(L.shape))
print(L)

print("\nUpper Triangular Matrix U ({}):".format(U.shape))
print(U)

# Verify: P @ B should be close to L @ U
print("\nIs P @ B close to L @ U?", np.allclose(P_mat @ B, L @ U))

Understanding the Output:

• P_mat: The permutation matrix $P$. Applying this to $A$ (as $P @ A$) represents the row swaps performed during Gaussian elimination.
• L: The lower triangular matrix with ones on the diagonal.
• U: The upper triangular matrix.

LU decomposition is a foundation for efficient direct solvers for systems of linear equations ($Ax=b$), although functions like scipy.linalg.solve often abstract this away.

QR Decomposition

QR decomposition factors any $m \times n$ matrix $A$ into $Q$ (an orthogonal matrix, $Q^T Q = I$) and $R$ (an upper/right triangular matrix), such that $A = QR$. scipy.linalg.qr is the function to use.

# Use the original matrix A again
# A = np.array([[1, 2, 3],
#               [4, 5, 6],
#               [7, 8, 9],
#               [10, 11, 12]])
print("Original Matrix A ({}):".format(A.shape))
print(A)

# Perform QR decomposition
# mode='economic' is efficient for non-square matrices
Q, R = linalg.qr(A, mode='economic')

print("\nOrthogonal Matrix Q ({}):".format(Q.shape))
print(Q)

print("\nUpper Triangular Matrix R ({}):".format(R.shape))
print(R)

# Verify Orthogonality of Q: Q.T @ Q should be close to Identity matrix
I = np.eye(Q.shape[1]) # Identity matrix of appropriate size
print("\nIs Q.T @ Q close to Identity?", np.allclose(Q.T @ Q, I))

# Verify Decomposition: A should be close to Q @ R
print("\nIs A close to Q @ R?", np.allclose(A, Q @ R))

Understanding the Output:

• Q: The orthogonal (or unitary for complex matrices) matrix. When mode='economic', if $m \ge n$, $Q$ is $m \times n$. If $m < n$, $Q$ is $m \times m$. Its columns form an orthonormal basis for the column space of $A$.
• R: The upper triangular matrix. When mode='economic', if $m \ge n$, $R$ is $n \times n$. If $m < n$, $R$ is $m \times n$.

QR decomposition is fundamental in algorithms for solving least-squares problems, eigenvalue computations (QR algorithm), and orthogonalization processes.

These SciPy functions provide tools for applying matrix decompositions. They handle the numerical complexities, allowing you to focus on how these decompositions can be used to solve problems in machine learning, such as dimensionality reduction (SVD), solving linear systems (LU), or finding least-squares solutions (QR). The next section provides a hands-on example focusing on SVD.