Array and Matrix Operations

In scientific computing, array and matrix operations are the most fundamental and important skills. SciPy provides richer and more efficient linear algebra functionality through the scipy.linalg module, building upon NumPy. This chapter will introduce the usage of these features in detail.

Overview of scipy.linalg Module

The scipy.linalg module includes all functions in numpy.linalg, and provides additional advanced features:

• More complete set of linear algebra functions
• Better performance optimization
• More matrix decomposition methods
• Handling of special matrices
• Matrix function computation

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt

# View main functions of scipy.linalg
print("Main functions of scipy.linalg:")
functions = [attr for attr in dir(linalg) if not attr.startswith('_')]
print(f"Total {len(functions)} functions and classes")
print("Main functions:", functions[:10])

Basic Matrix Operations

1. Matrix Creation and Basic Properties

# Create various types of matrices
n = 4

# Random matrix
A = np.random.random((n, n))
print(f"Random matrix A:\n{A}")

# Symmetric matrix
S = A + A.T
print(f"\nSymmetric matrix S:\n{S}")

# Positive definite matrix
P = A @ A.T
print(f"\nPositive definite matrix P:\n{P}")

# Identity matrix
I = np.eye(n)
print(f"\nIdentity matrix I:\n{I}")

# Zero matrix
Z = np.zeros((n, n))
print(f"\nZero matrix Z:\n{Z}")

# Diagonal matrix
D = np.diag([1, 2, 3, 4])
print(f"\nDiagonal matrix D:\n{D}")

2. Basic Matrix Operations

# Matrix addition and subtraction
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])

print(f"Matrix A:\n{A}")
print(f"Matrix B:\n{B}")
print(f"A + B:\n{A + B}")
print(f"A - B:\n{A - B}")

# Matrix multiplication
print(f"\nMatrix multiplication A @ B:\n{A @ B}")
print(f"Element-wise multiplication A * B:\n{A * B}")

# Matrix power
print(f"\nMatrix square A²:\n{linalg.matrix_power(A, 2)}")
print(f"Matrix cube A³:\n{linalg.matrix_power(A, 3)}")

# Matrix transpose
print(f"\nMatrix transpose A.T:\n{A.T}")

# Matrix trace (sum of diagonal elements)
print(f"\nTrace of matrix tr(A): {np.trace(A)}")

3. Matrix Norms

# Different types of matrix norms
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

print(f"Matrix A:\n{A}")
print(f"\nVarious norms:")
print(f"Frobenius norm: {linalg.norm(A, 'fro'):.4f}")
print(f"1-norm (max column sum): {linalg.norm(A, 1):.4f}")
print(f"∞-norm (max row sum): {linalg.norm(A, np.inf):.4f}")
print(f"2-norm (largest singular value): {linalg.norm(A, 2):.4f}")
print(f"Nuclear norm (sum of singular values): {linalg.norm(A, 'nuc'):.4f}")

# Vector norms
v = np.array([3, 4, 5])
print(f"\nVector v: {v}")
print(f"L1 norm: {linalg.norm(v, 1):.4f}")
print(f"L2 norm: {linalg.norm(v, 2):.4f}")
print(f"L∞ norm: {linalg.norm(v, np.inf):.4f}")

Solving Linear Systems

1. Basic Linear Systems

# Solve linear system Ax = b
A = np.array([[3, 2, -1], [2, -2, 4], [-1, 0.5, -1]])
b = np.array([1, -2, 0])

print(f"Coefficient matrix A:\n{A}")
print(f"Constant vector b: {b}")

# Using scipy.linalg.solve
x = linalg.solve(A, b)
print(f"\nSolution vector x: {x}")

# Verify correctness of solution
residual = A @ x - b
print(f"Residual ||Ax - b||: {linalg.norm(residual):.2e}")

# Check matrix condition number
cond_num = linalg.cond(A)
print(f"Condition number: {cond_num:.2f}")
if cond_num > 1e12:
    print("Warning: Matrix is nearly singular, solution may be unstable")

2. Special Types of Linear Systems

# Solving symmetric positive definite matrices (using Cholesky decomposition)
n = 5
A_spd = np.random.random((n, n))
A_spd = A_spd @ A_spd.T  # Ensure positive definiteness
b = np.random.random(n)

# General solver
start_time = time.time()
x1 = linalg.solve(A_spd, b)
general_time = time.time() - start_time

# Utilizing positive definiteness
start_time = time.time()
x2 = linalg.solve(A_spd, b, assume_a='pos')
cholesky_time = time.time() - start_time

print(f"General solve time: {general_time:.6f} seconds")
print(f"Cholesky solve time: {cholesky_time:.6f} seconds")
print(f"Performance improvement: {general_time / cholesky_time:.2f}x")
print(f"Solution difference: {linalg.norm(x1 - x2):.2e}")

# Tridiagonal matrix solver
from scipy.linalg import solve_banded

# Create tridiagonal matrix
n = 100
diag_main = 2 * np.ones(n)
diag_upper = -1 * np.ones(n-1)
diag_lower = -1 * np.ones(n-1)

# Construct banded matrix format
ab = np.zeros((3, n))
ab[0, 1:] = diag_upper  # Upper diagonal
ab[1, :] = diag_main    # Main diagonal
ab[2, :-1] = diag_lower # Lower diagonal

b_tri = np.random.random(n)

# Using banded matrix solver
x_banded = solve_banded((1, 1), ab, b_tri)
print(f"\nTridiagonal matrix solve complete, solution vector length: {len(x_banded)}")

3. Least Squares Problems

# Least squares solution for overdetermined linear systems
m, n = 100, 50  # m > n, overdetermined system
A_over = np.random.random((m, n))
b_over = np.random.random(m)

# Least squares solution
x_lstsq, residuals, rank, s = linalg.lstsq(A_over, b_over)

print(f"Overdetermined system: {m} equations, {n} unknowns")
print(f"Matrix rank: {rank}")
print(f"Sum of squared residuals: {residuals[0]:.4f}")
print(f"Solution vector norm: {linalg.norm(x_lstsq):.4f}")

# Regularized least squares (Ridge regression)
alpha = 0.1  # Regularization parameter
A_reg = np.vstack([A_over, alpha * np.eye(n)])
b_reg = np.hstack([b_over, np.zeros(n)])
x_ridge = linalg.lstsq(A_reg, b_reg)[0]

print(f"\nRegularized solution vector norm: {linalg.norm(x_ridge):.4f}")
print(f"Regularization effect: norm reduced by {(1 - linalg.norm(x_ridge)/linalg.norm(x_lstsq))*100:.1f}%")

Matrix Decomposition

1. LU Decomposition

# LU decomposition: A = PLU
A = np.array([[2, 1, 1], [1, 3, 2], [1, 0, 0]])

P, L, U = linalg.lu(A)

print(f"Original matrix A:\n{A}")
print(f"\nPermutation matrix P:\n{P}")
print(f"\nLower triangular matrix L:\n{L}")
print(f"\nUpper triangular matrix U:\n{U}")

# Verify decomposition
reconstructed = P @ L @ U
print(f"\nReconstructed matrix PLU:\n{reconstructed}")
print(f"Reconstruction error: {linalg.norm(A - reconstructed):.2e}")

# Using LU decomposition to solve linear system
lu, piv = linalg.lu_factor(A)
b = np.array([4, 5, 6])
x = linalg.lu_solve((lu, piv), b)
print(f"\nSolving Ax = b using LU decomposition:")
print(f"Solution: {x}")
print(f"Verification: Ax = {A @ x}")

2. QR Decomposition

# QR decomposition: A = QR
A = np.random.random((5, 3))

Q, R = linalg.qr(A)

print(f"Original matrix A shape: {A.shape}")
print(f"Orthogonal matrix Q shape: {Q.shape}")
print(f"Upper triangular matrix R shape: {R.shape}")

# Verify orthogonality
print(f"\nVerifying orthogonality of Q:")
print(f"Q.T @ Q:\n{Q.T @ Q}")
print(f"Close to identity: {np.allclose(Q.T @ Q, np.eye(Q.shape[1]))}")

# Verify decomposition
reconstructed = Q @ R
print(f"\nReconstruction error: {linalg.norm(A - reconstructed):.2e}")

# Application of QR decomposition: least squares problem
m, n = 10, 5
A_qr = np.random.random((m, n))
b_qr = np.random.random(m)

Q, R = linalg.qr(A_qr)
# Solve Rx = Q.T @ b
x_qr = linalg.solve_triangular(R, Q.T @ b_qr)

print(f"\nSolving least squares problem using QR decomposition:")
print(f"Solution vector: {x_qr}")
print(f"Residual: {linalg.norm(A_qr @ x_qr - b_qr):.4f}")

3. Singular Value Decomposition (SVD)

# SVD decomposition: A = UΣV^T
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9], [10, 11, 12]])

U, s, Vt = linalg.svd(A)

print(f"Original matrix A shape: {A.shape}")
print(f"U matrix shape: {U.shape}")
print(f"Singular value vector s: {s}")
print(f"V^T matrix shape: {Vt.shape}")

# Reconstruct matrix
S = np.zeros(A.shape)
S[:len(s), :len(s)] = np.diag(s)
reconstructed = U @ S @ Vt

print(f"\nReconstruction error: {linalg.norm(A - reconstructed):.2e}")

# SVD application: matrix rank and condition number
rank = np.sum(s > 1e-10)
cond_num = s[0] / s[-1] if s[-1] > 1e-15 else np.inf

print(f"\nNumerical rank of matrix: {rank}")
print(f"Condition number: {cond_num:.2f}")

# Low-rank approximation
for k in [1, 2, 3]:
    A_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    error = linalg.norm(A - A_k, 'fro')
    print(f"Rank-{k} approximation error: {error:.4f}")

4. Cholesky Decomposition

# Cholesky decomposition: A = L @ L.T (only for positive definite matrices)
n = 4
A_random = np.random.random((n, n))
A_pd = A_random @ A_random.T  # Construct positive definite matrix

try:
    L = linalg.cholesky(A_pd, lower=True)
    print(f"Positive definite matrix A:\n{A_pd}")
    print(f"\nCholesky factor L:\n{L}")

    # Verify decomposition
    reconstructed = L @ L.T
    print(f"\nReconstruction error: {linalg.norm(A_pd - reconstructed):.2e}")

    # Using Cholesky decomposition to solve linear system
    b = np.random.random(n)

    # Two-step solution: Ly = b, L^T x = y
    y = linalg.solve_triangular(L, b, lower=True)
    x = linalg.solve_triangular(L.T, y, lower=False)

    print(f"\nSolving using Cholesky decomposition:")
    print(f"Solution: {x}")
    print(f"Verification: {linalg.norm(A_pd @ x - b):.2e}")

except linalg.LinAlgError:
    print("Matrix is not positive definite, cannot perform Cholesky decomposition")

Eigenvalues and Eigenvectors

1. Eigenvalue Decomposition of General Matrices

# Calculate eigenvalues and eigenvectors
A = np.array([[1, 2], [2, 1]])

eigenvalues, eigenvectors = linalg.eig(A)

print(f"Matrix A:\n{A}")
print(f"\nEigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")

# Verify eigenvalue equation Av = λv
for i in range(len(eigenvalues)):
    lambda_i = eigenvalues[i]
    v_i = eigenvectors[:, i]

    Av = A @ v_i
    lambda_v = lambda_i * v_i

    print(f"\nEigenvalue {i+1}: λ = {lambda_i:.4f}")
    print(f"Av = {Av}")
    print(f"λv = {lambda_v}")
    print(f"Error: {linalg.norm(Av - lambda_v):.2e}")

# Matrix diagonalization
D = np.diag(eigenvalues)
P = eigenvectors
A_reconstructed = P @ D @ linalg.inv(P)

print(f"\nDiagonalization verification:")
print(f"Reconstructed matrix:\n{A_reconstructed}")
print(f"Reconstruction error: {linalg.norm(A - A_reconstructed):.2e}")

2. Eigenvalue Decomposition of Symmetric Matrices

# Symmetric matrices have real eigenvalues and orthogonal eigenvectors
A_sym = np.array([[4, 2, 1], [2, 3, 0], [1, 0, 2]])

# Using specialized algorithm for symmetric matrices
eigenvalues_sym, eigenvectors_sym = linalg.eigh(A_sym)

print(f"Symmetric matrix A:\n{A_sym}")
print(f"\nEigenvalues (in ascending order): {eigenvalues_sym}")
print(f"Eigenvector matrix:\n{eigenvectors_sym}")

# Verify orthogonality
print(f"\nOrthogonality of eigenvectors:")
print(f"V^T @ V:\n{eigenvectors_sym.T @ eigenvectors_sym}")
print(f"Is orthogonal: {np.allclose(eigenvectors_sym.T @ eigenvectors_sym, np.eye(3))}")

# Spectral decomposition
spectral_decomp = eigenvectors_sym @ np.diag(eigenvalues_sym) @ eigenvectors_sym.T
print(f"\nSpectral decomposition reconstruction error: {linalg.norm(A_sym - spectral_decomp):.2e}")

3. Generalized Eigenvalue Problem

# Generalized eigenvalue problem: Av = λBv
A = np.array([[1, 2], [3, 4]])
B = np.array([[2, 1], [1, 2]])

eigenvalues_gen, eigenvectors_gen = linalg.eig(A, B)

print(f"Matrix A:\n{A}")
print(f"Matrix B:\n{B}")
print(f"\nGeneralized eigenvalues: {eigenvalues_gen}")
print(f"Generalized eigenvectors:\n{eigenvectors_gen}")

# Verify generalized eigenvalue equation
for i in range(len(eigenvalues_gen)):
    lambda_i = eigenvalues_gen[i]
    v_i = eigenvectors_gen[:, i]

    Av = A @ v_i
    lambda_Bv = lambda_i * (B @ v_i)

    print(f"\nGeneralized eigenvalue {i+1}: λ = {lambda_i:.4f}")
    print(f"Error: {linalg.norm(Av - lambda_Bv):.2e}")

Matrix Functions

1. Matrix Exponential

# Matrix exponential exp(A)
A = np.array([[1, 2], [0, 1]])

# Calculate matrix exponential
exp_A = linalg.expm(A)

print(f"Matrix A:\n{A}")
print(f"\nMatrix exponential exp(A):\n{exp_A}")

# Verify property: exp(A + B) = exp(A) * exp(B) (when AB = BA)
B = np.array([[0, 1], [0, 0]])
if np.allclose(A @ B, B @ A):
    exp_A_plus_B = linalg.expm(A + B)
    exp_A_exp_B = linalg.expm(A) @ linalg.expm(B)

    print(f"\nexp(A + B):\n{exp_A_plus_B}")
    print(f"exp(A) * exp(B):\n{exp_A_exp_B}")
    print(f"Difference: {linalg.norm(exp_A_plus_B - exp_A_exp_B):.2e}")

# Application of matrix exponential: solving linear differential equation systems
# dx/dt = Ax, x(0) = x0
# Solution is x(t) = exp(At) * x0
t = 1.0
x0 = np.array([1, 0])
exp_At = linalg.expm(A * t)
x_t = exp_At @ x0

print(f"\nSolution of linear differential equation system:")
print(f"Solution at t = {t}: {x_t}")

2. Matrix Logarithm and Power Functions

# Matrix logarithm
A_pos = np.array([[2, 1], [0, 3]])  # Ensure positive eigenvalues

log_A = linalg.logm(A_pos)
print(f"Matrix A:\n{A_pos}")
print(f"\nMatrix logarithm log(A):\n{log_A}")

# Verify: exp(log(A)) = A
exp_log_A = linalg.expm(log_A)
print(f"\nexp(log(A)):\n{exp_log_A}")
print(f"Reconstruction error: {linalg.norm(A_pos - exp_log_A):.2e}")

# Matrix power function
for p in [0.5, 2, -1]:
    A_power = linalg.fractional_matrix_power(A_pos, p)
    print(f"\nA^{p}:\n{A_power}")

    if p == -1:
        # Verify inverse matrix
        A_inv_check = linalg.inv(A_pos)
        print(f"Inverse matrix calculated using inv():\n{A_inv_check}")
        print(f"Difference: {linalg.norm(A_power - A_inv_check):.2e}")

3. Matrix Square Root

# Matrix square root
A_sqrt = linalg.sqrtm(A_pos)
print(f"Matrix square root sqrt(A):\n{A_sqrt}")

# Verify: sqrt(A) * sqrt(A) = A
A_reconstructed = A_sqrt @ A_sqrt
print(f"\nsqrt(A) * sqrt(A):\n{A_reconstructed}")
print(f"Reconstruction error: {linalg.norm(A_pos - A_reconstructed):.2e}")

# For symmetric positive definite matrices, can use eigenvalue decomposition to calculate square root
A_spd = np.array([[4, 2], [2, 3]])
evals, evecs = linalg.eigh(A_spd)
A_sqrt_eig = evecs @ np.diag(np.sqrt(evals)) @ evecs.T

print(f"\nSquare root calculated using eigenvalue decomposition:\n{A_sqrt_eig}")
print(f"Verification: {linalg.norm(A_sqrt_eig @ A_sqrt_eig - A_spd):.2e}")

Sparse Matrix Operations

1. Creation and Basic Operations of Sparse Matrices

from scipy import sparse

# Create sparse matrix
n = 1000
# Create tridiagonal matrix
diagonals = [np.ones(n-1), 2*np.ones(n), np.ones(n-1)]
offsets = [-1, 0, 1]
A_sparse = sparse.diags(diagonals, offsets, format='csr')

print(f"Sparse matrix shape: {A_sparse.shape}")
print(f"Number of non-zero elements: {A_sparse.nnz}")
print(f"Sparsity: {1 - A_sparse.nnz / (n * n):.4f}")

# Solving linear systems with sparse matrices
from scipy.sparse.linalg import spsolve

b_sparse = np.random.random(n)
x_sparse = spsolve(A_sparse, b_sparse)

print(f"\nSparse linear system solved")
print(f"Solution vector norm: {linalg.norm(x_sparse):.4f}")
print(f"Residual: {linalg.norm(A_sparse @ x_sparse - b_sparse):.2e}")

# Eigenvalue calculation for sparse matrices
from scipy.sparse.linalg import eigsh

# Calculate largest few eigenvalues
k = 5  # Calculate first 5 eigenvalues
eigenvalues_sparse, eigenvectors_sparse = eigsh(A_sparse, k=k, which='LA')

print(f"\nLargest {k} eigenvalues: {eigenvalues_sparse}")
print(f"Corresponding eigenvectors shape: {eigenvectors_sparse.shape}")

Practical Application Examples

Example 1: Image Compression (SVD Application)

(Content omitted for brevity - includes SVD image compression code)

Example 2: Principal Component Analysis (PCA)

(Content omitted for brevity - includes PCA implementation code)

Example 3: Linear Regression in Matrix Form

(Content omitted for brevity - includes linear regression code)

Performance Optimization Tips

1. Choosing Appropriate Algorithms

(Content omitted for brevity - includes performance comparison code)

2. Memory Optimization

(Content omitted for brevity - includes memory optimization code)

Summary

In this chapter, we thoroughly learned about SciPy's array and matrix operations:

1. Basic Operations: Matrix creation, operations, norm calculations
2. Linear Systems: Various solution methods and special case handling
3. Matrix Decomposition: LU, QR, SVD, Cholesky decomposition and applications
4. Eigenvalue Problems: Eigenvalue decomposition for general and symmetric matrices
5. Matrix Functions: Exponential, logarithm, power functions, etc.
6. Sparse Matrices: Efficient handling of large-scale sparse problems
7. Practical Applications: Image compression, PCA, linear regression, etc.
8. Performance Optimization: Algorithm selection and memory management

These skills are fundamental to scientific computing and provide a solid foundation for learning more advanced SciPy features. Next, we will learn about SciPy's powerful statistical features in Statistical Analysis.

Exercises

1. Matrix Decomposition Application: Implement a function that uses different matrix decomposition methods to solve linear systems and compare their performance
2. Eigenvalue Application: Write code to analyze network connectivity (using eigenvalues of the graph's adjacency matrix)
3. SVD Application: Implement a simple recommendation system using SVD for matrix factorization
4. Performance Optimization: Compare performance differences between dense and sparse matrix operations

# Reference answers for exercises

# 1. Matrix decomposition application
def solve_comparison(A, b):
    """
    Solve linear systems using different methods and compare performance
    """
    methods = {}

    # Direct solve
    start = time.time()
    x_direct = linalg.solve(A, b)
    methods['Direct solve'] = (time.time() - start, x_direct)

    # LU decomposition
    start = time.time()
    lu, piv = linalg.lu_factor(A)
    x_lu = linalg.lu_solve((lu, piv), b)
    methods['LU decomposition'] = (time.time() - start, x_lu)

    # QR decomposition
    start = time.time()
    Q, R = linalg.qr(A)
    x_qr = linalg.solve_triangular(R, Q.T @ b)
    methods['QR decomposition'] = (time.time() - start, x_qr)

    return methods

# Test
n = 500
A_test = np.random.random((n, n))
b_test = np.random.random(n)
results = solve_comparison(A_test, b_test)

for method, (time_taken, solution) in results.items():
    error = linalg.norm(A_test @ solution - b_test)
    print(f"{method}: {time_taken:.4f}s, error: {error:.2e}")