9.4 KiB
9.4 KiB
SymPy Matrices and Linear Algebra
This document covers SymPy's matrix operations, linear algebra capabilities, and solving systems of linear equations.
Matrix Creation
Basic Matrix Construction
from sympy import Matrix, eye, zeros, ones, diag
# From list of rows
M = Matrix([[1, 2], [3, 4]])
M = Matrix([
[1, 2, 3],
[4, 5, 6]
])
# Column vector
v = Matrix([1, 2, 3])
# Row vector
v = Matrix([[1, 2, 3]])
Special Matrices
# Identity matrix
I = eye(3) # 3x3 identity
# [[1, 0, 0],
# [0, 1, 0],
# [0, 0, 1]]
# Zero matrix
Z = zeros(2, 3) # 2 rows, 3 columns of zeros
# Ones matrix
O = ones(3, 2) # 3 rows, 2 columns of ones
# Diagonal matrix
D = diag(1, 2, 3)
# [[1, 0, 0],
# [0, 2, 0],
# [0, 0, 3]]
# Block diagonal
from sympy import BlockDiagMatrix
A = Matrix([[1, 2], [3, 4]])
B = Matrix([[5, 6], [7, 8]])
BD = BlockDiagMatrix(A, B)
Matrix Properties and Access
Shape and Dimensions
M = Matrix([[1, 2, 3], [4, 5, 6]])
M.shape # (2, 3) - returns tuple (rows, cols)
M.rows # 2
M.cols # 3
Accessing Elements
M = Matrix([[1, 2, 3], [4, 5, 6]])
# Single element
M[0, 0] # 1 (zero-indexed)
M[1, 2] # 6
# Row access
M[0, :] # Matrix([[1, 2, 3]])
M.row(0) # Same as above
# Column access
M[:, 1] # Matrix([[2], [5]])
M.col(1) # Same as above
# Slicing
M[0:2, 0:2] # Top-left 2x2 submatrix
Modification
M = Matrix([[1, 2], [3, 4]])
# Insert row
M = M.row_insert(1, Matrix([[5, 6]]))
# [[1, 2],
# [5, 6],
# [3, 4]]
# Insert column
M = M.col_insert(1, Matrix([7, 8]))
# Delete row
M = M.row_del(0)
# Delete column
M = M.col_del(1)
Basic Matrix Operations
Arithmetic Operations
from sympy import Matrix
A = Matrix([[1, 2], [3, 4]])
B = Matrix([[5, 6], [7, 8]])
# Addition
C = A + B
# Subtraction
C = A - B
# Scalar multiplication
C = 2 * A
# Matrix multiplication
C = A * B
# Element-wise multiplication (Hadamard product)
C = A.multiply_elementwise(B)
# Power
C = A**2 # Same as A * A
C = A**3 # A * A * A
Transpose and Conjugate
M = Matrix([[1, 2], [3, 4]])
# Transpose
M.T
# [[1, 3],
# [2, 4]]
# Conjugate (for complex matrices)
M.conjugate()
# Conjugate transpose (Hermitian transpose)
M.H # Same as M.conjugate().T
Inverse
M = Matrix([[1, 2], [3, 4]])
# Inverse
M_inv = M**-1
M_inv = M.inv()
# Verify
M * M_inv # Returns identity matrix
# Check if invertible
M.is_invertible() # True or False
Advanced Linear Algebra
Determinant
M = Matrix([[1, 2], [3, 4]])
M.det() # -2
# For symbolic matrices
from sympy import symbols
a, b, c, d = symbols('a b c d')
M = Matrix([[a, b], [c, d]])
M.det() # a*d - b*c
Trace
M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
M.trace() # 1 + 5 + 9 = 15
Row Echelon Form
M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
# Reduced Row Echelon Form
rref_M, pivot_cols = M.rref()
# rref_M is the RREF matrix
# pivot_cols is tuple of pivot column indices
Rank
M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
M.rank() # 2 (this matrix is rank-deficient)
Nullspace and Column Space
M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
# Nullspace (kernel)
null = M.nullspace()
# Returns list of basis vectors for nullspace
# Column space
col = M.columnspace()
# Returns list of basis vectors for column space
# Row space
row = M.rowspace()
# Returns list of basis vectors for row space
Orthogonalization
# Gram-Schmidt orthogonalization
vectors = [Matrix([1, 2, 3]), Matrix([4, 5, 6])]
ortho = Matrix.orthogonalize(*vectors)
# With normalization
ortho_norm = Matrix.orthogonalize(*vectors, normalize=True)
Eigenvalues and Eigenvectors
Computing Eigenvalues
M = Matrix([[1, 2], [2, 1]])
# Eigenvalues with multiplicities
eigenvals = M.eigenvals()
# Returns dict: {eigenvalue: multiplicity}
# Example: {3: 1, -1: 1}
# Just the eigenvalues as a list
eigs = list(M.eigenvals().keys())
Computing Eigenvectors
M = Matrix([[1, 2], [2, 1]])
# Eigenvectors with eigenvalues
eigenvects = M.eigenvects()
# Returns list of tuples: (eigenvalue, multiplicity, [eigenvectors])
# Example: [(3, 1, [Matrix([1, 1])]), (-1, 1, [Matrix([1, -1])])]
# Access individual eigenvectors
for eigenval, multiplicity, eigenvecs in M.eigenvects():
print(f"Eigenvalue: {eigenval}")
print(f"Eigenvectors: {eigenvecs}")
Diagonalization
M = Matrix([[1, 2], [2, 1]])
# Check if diagonalizable
M.is_diagonalizable() # True or False
# Diagonalize (M = P*D*P^-1)
P, D = M.diagonalize()
# P: matrix of eigenvectors
# D: diagonal matrix of eigenvalues
# Verify
P * D * P**-1 == M # True
Characteristic Polynomial
from sympy import symbols
lam = symbols('lambda')
M = Matrix([[1, 2], [2, 1]])
charpoly = M.charpoly(lam)
# Returns characteristic polynomial
Jordan Normal Form
M = Matrix([[2, 1, 0], [0, 2, 1], [0, 0, 2]])
# Jordan form (for non-diagonalizable matrices)
P, J = M.jordan_form()
# J is the Jordan normal form
# P is the transformation matrix
Matrix Decompositions
LU Decomposition
M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 10]])
# LU decomposition
L, U, perm = M.LUdecomposition()
# L: lower triangular
# U: upper triangular
# perm: permutation indices
QR Decomposition
M = Matrix([[1, 2], [3, 4], [5, 6]])
# QR decomposition
Q, R = M.QRdecomposition()
# Q: orthogonal matrix
# R: upper triangular matrix
Cholesky Decomposition
# For positive definite symmetric matrices
M = Matrix([[4, 2], [2, 3]])
L = M.cholesky()
# L is lower triangular such that M = L*L.T
Singular Value Decomposition (SVD)
M = Matrix([[1, 2], [3, 4], [5, 6]])
# SVD (note: may require numerical evaluation)
U, S, V = M.singular_value_decomposition()
# M = U * S * V
Solving Linear Systems
Using Matrix Equations
# Solve Ax = b
A = Matrix([[1, 2], [3, 4]])
b = Matrix([5, 6])
# Solution
x = A.solve(b) # or A**-1 * b
# Least squares (for overdetermined systems)
x = A.solve_least_squares(b)
Using linsolve
from sympy import linsolve, symbols
x, y = symbols('x y')
# Method 1: List of equations
eqs = [x + y - 5, 2*x - y - 1]
sol = linsolve(eqs, [x, y])
# {(2, 3)}
# Method 2: Augmented matrix
M = Matrix([[1, 1, 5], [2, -1, 1]])
sol = linsolve(M, [x, y])
# Method 3: A*x = b form
A = Matrix([[1, 1], [2, -1]])
b = Matrix([5, 1])
sol = linsolve((A, b), [x, y])
Underdetermined and Overdetermined Systems
# Underdetermined (infinite solutions)
A = Matrix([[1, 2, 3]])
b = Matrix([6])
sol = A.solve(b) # Returns parametric solution
# Overdetermined (least squares)
A = Matrix([[1, 2], [3, 4], [5, 6]])
b = Matrix([1, 2, 3])
sol = A.solve_least_squares(b)
Symbolic Matrices
Working with Symbolic Entries
from sympy import symbols, Matrix
a, b, c, d = symbols('a b c d')
M = Matrix([[a, b], [c, d]])
# All operations work symbolically
M.det() # a*d - b*c
M.inv() # Matrix([[d/(a*d - b*c), -b/(a*d - b*c)], ...])
M.eigenvals() # Symbolic eigenvalues
Matrix Functions
from sympy import exp, sin, cos, Matrix
M = Matrix([[0, 1], [-1, 0]])
# Matrix exponential
exp(M)
# Trigonometric functions
sin(M)
cos(M)
Mutable vs Immutable Matrices
from sympy import Matrix, ImmutableMatrix
# Mutable (default)
M = Matrix([[1, 2], [3, 4]])
M[0, 0] = 5 # Allowed
# Immutable (for use as dictionary keys, etc.)
I = ImmutableMatrix([[1, 2], [3, 4]])
# I[0, 0] = 5 # Error: ImmutableMatrix cannot be modified
Sparse Matrices
For large matrices with many zero entries:
from sympy import SparseMatrix
# Create sparse matrix
S = SparseMatrix(1000, 1000, {(0, 0): 1, (100, 100): 2})
# Only stores non-zero elements
# Convert dense to sparse
M = Matrix([[1, 0, 0], [0, 2, 0]])
S = SparseMatrix(M)
Common Linear Algebra Patterns
Pattern 1: Solving Ax = b for Multiple b Vectors
A = Matrix([[1, 2], [3, 4]])
A_inv = A.inv()
b1 = Matrix([5, 6])
b2 = Matrix([7, 8])
x1 = A_inv * b1
x2 = A_inv * b2
Pattern 2: Change of Basis
# Given vectors in old basis, convert to new basis
old_basis = [Matrix([1, 0]), Matrix([0, 1])]
new_basis = [Matrix([1, 1]), Matrix([1, -1])]
# Change of basis matrix
P = Matrix.hstack(*new_basis)
P_inv = P.inv()
# Convert vector v from old to new basis
v = Matrix([3, 4])
v_new = P_inv * v
Pattern 3: Matrix Condition Number
# Estimate condition number (ratio of largest to smallest singular value)
M = Matrix([[1, 2], [3, 4]])
eigenvals = M.eigenvals()
cond = max(eigenvals.keys()) / min(eigenvals.keys())
Pattern 4: Projection Matrices
# Project onto column space of A
A = Matrix([[1, 0], [0, 1], [1, 1]])
P = A * (A.T * A).inv() * A.T
# P is projection matrix onto column space of A
Important Notes
-
Zero-testing: SymPy's symbolic zero-testing can affect accuracy. For numerical work, consider using
evalf()or numerical libraries. -
Performance: For large numerical matrices, consider converting to NumPy using
lambdifyor using numerical linear algebra libraries directly. -
Symbolic computation: Matrix operations with symbolic entries can be computationally expensive for large matrices.
-
Assumptions: Use symbol assumptions (e.g.,
real=True,positive=True) to help SymPy simplify matrix expressions correctly.