Untitled
unknown
python
a year ago
12 kB
5
Indexable
from typing import List, Callable, TypeVar, Tuple, Optional, Union
from functools import reduce
from math import isinf, isnan
from common.utils import randomize_int, fast_power
from matrix.Vector import Vector
from matrix.Matrix import Matrix
T = TypeVar('T')
class SqMatrix(Matrix):
def __init__(self, content):
"""Set self.dim and self.content, assert that number of rows equal number of columns
Also initialize superclass properties on self."""
self.dim = len(content)
self.content = tuple(tuple(row) for row in content)
assert all(len(row) == self.dim for row in content)
super().__init__(content)
def __repr__(self):
"""Print a square matrix as 'SqMatrix(...)' with one row of content per line,
and indentation."""
lines = ",\n ".join(f"{row}" for row in self.content)
return f"SqMatrix(({lines}))"
@staticmethod
def random(dim: int, randomizer: Callable[[], int] = randomize_int()) -> 'SqMatrix':
"""Create and return a new square matrix of the given size, being filled
with randomly selected elements. Each such element is the return
value of the given randomizer function. randomizer will be called
once for each entry in the resulting matrix."""
content = [[randomizer() for _ in range(dim)] for _ in range(dim)]
return SqMatrix(content)
@staticmethod
def identity(dim: int) -> 'SqMatrix':
"""Return a square matrix of the given size with 1 (integer)
along the main diagonal, and 0 (integer) off the main diagonal."""
content = [[1 if i == j else 0 for j in range(dim)] for i in range(dim)]
return SqMatrix(content)
@staticmethod
def zero(dim: int) -> 'SqMatrix':
"""Return a square matrix of the given size consisting entirely
of 0 (integer)."""
content = [[0 for _ in range(dim)] for _ in range(dim)]
return SqMatrix(content)
@staticmethod
def tabulate(dim: int, f: Callable[[int, int], float]) -> 'SqMatrix':
"""Create and return a new square matrix of the given size.
The SqMatrix is filled with the return value of the given function, f.
The function, f, is called once per entry in the resulting matrix,
and the return value of f(i,j) is the value of m[i][j]"""
content = [[f(i, j) for j in range(dim)] for i in range(dim)]
return SqMatrix(content)
@staticmethod
def diagonal(entries) -> 'SqMatrix':
"""Given a sequence of values, return a square matrix with those
values along the main diagonal. If the sequence has size n,
then the resulting matrix has dimension n as well. The diagonal
matrix has 0 (integer) off the main diagonal."""
dim = len(entries)
content = [[entries[i] if i == j else 0 for j in range(dim)] for i in range(dim)]
return SqMatrix(content)
def gaussian_elimination_back_substitution(self, v: Vector) -> Optional[Vector]:
"""Adjoin the given column vector to the matrix, self.
Then perform elementary row operations on the given matrix to reduce
it to row echelon form. Then perform back substitution to return
the vector which solves the system."""
assert isinstance(v, Vector)
assert v.dim == self.dim
# compute an echelon form of the matrix, self.
# this form has zeros below the main diagonal.
# and non-zeros on the main diagonal unless the matrix is
# singular in which case there are zeros on the diagonal.
ech = self.adjoin_col(v).make_row_echelon()
if ech is None:
return None
assert ech.rows == self.rows
assert ech.cols == self.cols + 1
return ech.back_substitution()
def gauss_jordan_elimination(self, v: Vector) -> Optional[Vector]:
"""Adjoin the given column vector to the matrix, self.
Then use Gauss Jordan elimination (via elementary row operations)
to reduce the matrix to diagonal. This may result in a matrix with
a zero on the diagonal, in which case None is returned, otherwise
normalize the diagonal and return the vector which solves the system.
"""
assert isinstance(v, Vector)
assert v.dim == self.dim
diag, _ = self.adjoin_col(v).make_unit_diagonal()
if diag is None:
return None
assert diag.rows == self.rows
assert diag.cols == self.cols + 1
return diag.col_vec(v.dim)
def gauss_jordan_inverse(self) -> Tuple[Optional['SqMatrix'], T]:
"""returns 2-tuple of two values: (inverse,determinant),
If the determinant is zero, then inverse=None
"""
diag, det = self.adjoin_cols(SqMatrix.identity(self.dim)).make_unit_diagonal()
if diag is None:
return None, 0
assert diag.rows == self.rows
assert diag.cols == 2 * self.cols
return diag.extract_cols(range(self.dim, 2 * self.dim)), det
def laplacian_expansion(self, zero: T = 0) -> T:
"""Compute the determinant of the given square matrix using
Laplacian expansion."""
if self.dim == 1:
return self[0][0]
else:
det = 0
for j in range(self.dim):
det += ((-1) ** j) * self[0][j] * self.minor(0, j).laplacian_expansion()
return det
def cramers_rule(self, b: Vector) -> Vector:
"""Return a new vector which solves the system of equations Ax=b,
where A is the matrix self. If the determinant of A is 0, then None
is returned. The solution is found using Cramer's Rule:
I.e., to compute the k'th component of the returned Vector,
we replaced the k'th column of a by the column vector b,
and calculate the determinant of that matrix, then divide by
the determinant of A.
"""
assert self.dim == b.dim
det_A = self.laplacian_expansion()
if det_A == 0:
return None
solutions = []
for i in range(self.dim):
mat = self.copy()
mat.set_col(i, b)
det_mat = mat.laplacian_expansion()
solutions.append(det_mat / det_A)
return Vector(*solutions)
def power(self, p: int) -> 'SqMatrix':
"""Raise the matrix to the p'th power."""
assert isinstance(p, int)
assert p >= 0
identity_matrix = SqMatrix.identity(self.dim)
return fast_power(self, lambda x, y: x * y, p, identity_matrix)
def characteristic_polynomial(self) -> 'Polynomial':
"""Compute the characteristic polynomial of the square matrix.
This is done by generating the matrix M - λI, which is a SqMatrix
having a Polynomial in each entry. Then use laplacian_expansion
to compute the determinant.
"""
from matrix.Polynomial import Polynomial
variable = Polynomial.variable()
identity_matrix = SqMatrix.identity(self.dim)
diff_matrix = self - (identity_matrix * variable)
return diff_matrix.laplacian_expansion().to_polynomial(variable)
def eigenvalues(self, epsilon: float) -> List[float]:
"""Compute the (some) eigenvalues of the matrix by calling
characteristic_polynomial to obtain a Polynomial, then using
the function poly_roots to compute its roots.
The eigenvalues method returns the roots/eigenvalues in increasing
order.
Just as poly_roots, this method is allowed to fail to compute
the eigenvalues in some cases.
"""
from algebra.roots import poly_roots
assert isinstance(epsilon, float), f"epsilon={epsilon} is not a float"
char_poly = self.characteristic_polynomial()
roots = poly_roots(char_poly, epsilon)
return sorted(roots)
def eigenvectors(self, spectrum: List[float], epsilon=0.0000001) -> List[Vector]:
"""Given a list of non-zero eigenvalues, `spectrum`,
return the corresponding set of eigenvectors.
It is assumed that spectrum is already sorted in increasing order.
If the given eigenvalues are not distinct, return None.
If any eigenvalue is zero (abs(e)<epsilon), return None.
If there are not exactly self.dim number of eigenvalues given,
then return None"""
assert spectrum == sorted(spectrum)
if len(spectrum) != self.dim:
return None
# check for zero eigenvalues
if any(abs(e) < epsilon for e in spectrum):
return None
# check for duplicate eigenvalues
for k in range(len(spectrum) - 1):
# we know spectrum[k+1] >= spectrum[k]
if spectrum[k + 1] - spectrum[k] < epsilon:
# we found two equal eigenvalues (within epsilon of each other)
return None
def condition(m):
"""replace numbers close to zero by zero"""
return SqMatrix.tabulate(self.dim,
lambda r, c: 0 if abs(m[r][c]) < epsilon else m[r][c])
id = SqMatrix.identity(self.dim)
amli = [condition(self - id.scale(spectrum[k]))
for k in range(self.dim)] # A - lambda[k] I
def prod(ms):
return reduce(lambda a, b: a * b,
ms,
id)
# pamli: product-of-A-minus-lambda-I
# here we form a list of matrix products.
# the k-th element in this list is the product of all
# elements of amli except the k'th one
pamli = [prod(ms)
for k in range(len(spectrum))
for ms in [[amli[j]
for j in range(len(spectrum))
if j != k]]]
# each element in pamli is a matrix having a column of zeros.
# we want to select a column vector from each, but we want to avoid
# the zero vectors. thus for each matrix in pamli, we find the column
# vector with the maximum norm. This works if all the eigenvalues
# or distinct; however, if there is a repeated eigenvalue, then this
# search will result in a repeated eigenvector. This is a bug which
# we need to fix.
vectors = []
for k in range(len(spectrum)):
# find all columns of pamli[k],
# excluding any zero vector in `vectors`
independent = [cv # column vector
for j in range(pamli[k].cols)
for cv in [pamli[k].col_vec(j)]
if cv.norm() > epsilon]
if independent:
vectors.append(min(independent, key=lambda v: v.norm()))
else:
vectors.append(None)
return vectors
def eigen(self, epsilon: float) -> Optional[List[Tuple[float, Vector]]]:
"""Compute (eigenvalue, eigenvector) pairs of the matrix,
if one (or more) of the eigenvalues is 0, return None.
By 'e is 0' we mean abs(e) < epsilon.
If we cannot find self.dim number of unique eigenvalues,
then return None.
Else return a list of tuples (e,v) where e is an eigenvalue
and v is its normalized eigenvector. Note that if v is a normalized
eigenvector, then so is -v, so use v or -v as you prefer.
The list is sorted e1 < e2 < ... < en.
"""
spectrum = self.eigenvalues(epsilon)
if len(set(spectrum)) != self.dim:
return None
if any(abs(e) < epsilon for e in spectrum):
return None
vectors = self.eigenvectors(spectrum)
return [(spectrum[k], vectors[k])
for k in range(len(vectors))]
def minor(self, i, j):
"""Return the minor of the matrix after removing the i-th row and j-th column."""
return SqMatrix([row[:j] + row[j+1:] for row in (self[:i]+self[i+1:])])
Editor is loading...
Leave a Comment