You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

327 lines
10 KiB

5 months ago
from types import FunctionType
from .utilities import _get_intermediate_simp, _iszero, _dotprodsimp, _simplify
from .determinant import _find_reasonable_pivot
def _row_reduce_list(mat, rows, cols, one, iszerofunc, simpfunc,
normalize_last=True, normalize=True, zero_above=True):
"""Row reduce a flat list representation of a matrix and return a tuple
(rref_matrix, pivot_cols, swaps) where ``rref_matrix`` is a flat list,
``pivot_cols`` are the pivot columns and ``swaps`` are any row swaps that
were used in the process of row reduction.
Parameters
==========
mat : list
list of matrix elements, must be ``rows`` * ``cols`` in length
rows, cols : integer
number of rows and columns in flat list representation
one : SymPy object
represents the value one, from ``Matrix.one``
iszerofunc : determines if an entry can be used as a pivot
simpfunc : used to simplify elements and test if they are
zero if ``iszerofunc`` returns `None`
normalize_last : indicates where all row reduction should
happen in a fraction-free manner and then the rows are
normalized (so that the pivots are 1), or whether
rows should be normalized along the way (like the naive
row reduction algorithm)
normalize : whether pivot rows should be normalized so that
the pivot value is 1
zero_above : whether entries above the pivot should be zeroed.
If ``zero_above=False``, an echelon matrix will be returned.
"""
def get_col(i):
return mat[i::cols]
def row_swap(i, j):
mat[i*cols:(i + 1)*cols], mat[j*cols:(j + 1)*cols] = \
mat[j*cols:(j + 1)*cols], mat[i*cols:(i + 1)*cols]
def cross_cancel(a, i, b, j):
"""Does the row op row[i] = a*row[i] - b*row[j]"""
q = (j - i)*cols
for p in range(i*cols, (i + 1)*cols):
mat[p] = isimp(a*mat[p] - b*mat[p + q])
isimp = _get_intermediate_simp(_dotprodsimp)
piv_row, piv_col = 0, 0
pivot_cols = []
swaps = []
# use a fraction free method to zero above and below each pivot
while piv_col < cols and piv_row < rows:
pivot_offset, pivot_val, \
assumed_nonzero, newly_determined = _find_reasonable_pivot(
get_col(piv_col)[piv_row:], iszerofunc, simpfunc)
# _find_reasonable_pivot may have simplified some things
# in the process. Let's not let them go to waste
for (offset, val) in newly_determined:
offset += piv_row
mat[offset*cols + piv_col] = val
if pivot_offset is None:
piv_col += 1
continue
pivot_cols.append(piv_col)
if pivot_offset != 0:
row_swap(piv_row, pivot_offset + piv_row)
swaps.append((piv_row, pivot_offset + piv_row))
# if we aren't normalizing last, we normalize
# before we zero the other rows
if normalize_last is False:
i, j = piv_row, piv_col
mat[i*cols + j] = one
for p in range(i*cols + j + 1, (i + 1)*cols):
mat[p] = isimp(mat[p] / pivot_val)
# after normalizing, the pivot value is 1
pivot_val = one
# zero above and below the pivot
for row in range(rows):
# don't zero our current row
if row == piv_row:
continue
# don't zero above the pivot unless we're told.
if zero_above is False and row < piv_row:
continue
# if we're already a zero, don't do anything
val = mat[row*cols + piv_col]
if iszerofunc(val):
continue
cross_cancel(pivot_val, row, val, piv_row)
piv_row += 1
# normalize each row
if normalize_last is True and normalize is True:
for piv_i, piv_j in enumerate(pivot_cols):
pivot_val = mat[piv_i*cols + piv_j]
mat[piv_i*cols + piv_j] = one
for p in range(piv_i*cols + piv_j + 1, (piv_i + 1)*cols):
mat[p] = isimp(mat[p] / pivot_val)
return mat, tuple(pivot_cols), tuple(swaps)
# This functions is a candidate for caching if it gets implemented for matrices.
def _row_reduce(M, iszerofunc, simpfunc, normalize_last=True,
normalize=True, zero_above=True):
mat, pivot_cols, swaps = _row_reduce_list(list(M), M.rows, M.cols, M.one,
iszerofunc, simpfunc, normalize_last=normalize_last,
normalize=normalize, zero_above=zero_above)
return M._new(M.rows, M.cols, mat), pivot_cols, swaps
def _is_echelon(M, iszerofunc=_iszero):
"""Returns `True` if the matrix is in echelon form. That is, all rows of
zeros are at the bottom, and below each leading non-zero in a row are
exclusively zeros."""
if M.rows <= 0 or M.cols <= 0:
return True
zeros_below = all(iszerofunc(t) for t in M[1:, 0])
if iszerofunc(M[0, 0]):
return zeros_below and _is_echelon(M[:, 1:], iszerofunc)
return zeros_below and _is_echelon(M[1:, 1:], iszerofunc)
def _echelon_form(M, iszerofunc=_iszero, simplify=False, with_pivots=False):
"""Returns a matrix row-equivalent to ``M`` that is in echelon form. Note
that echelon form of a matrix is *not* unique, however, properties like the
row space and the null space are preserved.
Examples
========
>>> from sympy import Matrix
>>> M = Matrix([[1, 2], [3, 4]])
>>> M.echelon_form()
Matrix([
[1, 2],
[0, -2]])
"""
simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
mat, pivots, _ = _row_reduce(M, iszerofunc, simpfunc,
normalize_last=True, normalize=False, zero_above=False)
if with_pivots:
return mat, pivots
return mat
# This functions is a candidate for caching if it gets implemented for matrices.
def _rank(M, iszerofunc=_iszero, simplify=False):
"""Returns the rank of a matrix.
Examples
========
>>> from sympy import Matrix
>>> from sympy.abc import x
>>> m = Matrix([[1, 2], [x, 1 - 1/x]])
>>> m.rank()
2
>>> n = Matrix(3, 3, range(1, 10))
>>> n.rank()
2
"""
def _permute_complexity_right(M, iszerofunc):
"""Permute columns with complicated elements as
far right as they can go. Since the ``sympy`` row reduction
algorithms start on the left, having complexity right-shifted
speeds things up.
Returns a tuple (mat, perm) where perm is a permutation
of the columns to perform to shift the complex columns right, and mat
is the permuted matrix."""
def complexity(i):
# the complexity of a column will be judged by how many
# element's zero-ness cannot be determined
return sum(1 if iszerofunc(e) is None else 0 for e in M[:, i])
complex = [(complexity(i), i) for i in range(M.cols)]
perm = [j for (i, j) in sorted(complex)]
return (M.permute(perm, orientation='cols'), perm)
simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
# for small matrices, we compute the rank explicitly
# if is_zero on elements doesn't answer the question
# for small matrices, we fall back to the full routine.
if M.rows <= 0 or M.cols <= 0:
return 0
if M.rows <= 1 or M.cols <= 1:
zeros = [iszerofunc(x) for x in M]
if False in zeros:
return 1
if M.rows == 2 and M.cols == 2:
zeros = [iszerofunc(x) for x in M]
if False not in zeros and None not in zeros:
return 0
d = M.det()
if iszerofunc(d) and False in zeros:
return 1
if iszerofunc(d) is False:
return 2
mat, _ = _permute_complexity_right(M, iszerofunc=iszerofunc)
_, pivots, _ = _row_reduce(mat, iszerofunc, simpfunc, normalize_last=True,
normalize=False, zero_above=False)
return len(pivots)
def _rref(M, iszerofunc=_iszero, simplify=False, pivots=True,
normalize_last=True):
"""Return reduced row-echelon form of matrix and indices of pivot vars.
Parameters
==========
iszerofunc : Function
A function used for detecting whether an element can
act as a pivot. ``lambda x: x.is_zero`` is used by default.
simplify : Function
A function used to simplify elements when looking for a pivot.
By default SymPy's ``simplify`` is used.
pivots : True or False
If ``True``, a tuple containing the row-reduced matrix and a tuple
of pivot columns is returned. If ``False`` just the row-reduced
matrix is returned.
normalize_last : True or False
If ``True``, no pivots are normalized to `1` until after all
entries above and below each pivot are zeroed. This means the row
reduction algorithm is fraction free until the very last step.
If ``False``, the naive row reduction procedure is used where
each pivot is normalized to be `1` before row operations are
used to zero above and below the pivot.
Examples
========
>>> from sympy import Matrix
>>> from sympy.abc import x
>>> m = Matrix([[1, 2], [x, 1 - 1/x]])
>>> m.rref()
(Matrix([
[1, 0],
[0, 1]]), (0, 1))
>>> rref_matrix, rref_pivots = m.rref()
>>> rref_matrix
Matrix([
[1, 0],
[0, 1]])
>>> rref_pivots
(0, 1)
``iszerofunc`` can correct rounding errors in matrices with float
values. In the following example, calling ``rref()`` leads to
floating point errors, incorrectly row reducing the matrix.
``iszerofunc= lambda x: abs(x)<1e-9`` sets sufficiently small numbers
to zero, avoiding this error.
>>> m = Matrix([[0.9, -0.1, -0.2, 0], [-0.8, 0.9, -0.4, 0], [-0.1, -0.8, 0.6, 0]])
>>> m.rref()
(Matrix([
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0]]), (0, 1, 2))
>>> m.rref(iszerofunc=lambda x:abs(x)<1e-9)
(Matrix([
[1, 0, -0.301369863013699, 0],
[0, 1, -0.712328767123288, 0],
[0, 0, 0, 0]]), (0, 1))
Notes
=====
The default value of ``normalize_last=True`` can provide significant
speedup to row reduction, especially on matrices with symbols. However,
if you depend on the form row reduction algorithm leaves entries
of the matrix, set ``noramlize_last=False``
"""
simpfunc = simplify if isinstance(simplify, FunctionType) else _simplify
mat, pivot_cols, _ = _row_reduce(M, iszerofunc, simpfunc,
normalize_last, normalize=True, zero_above=True)
if pivots:
mat = (mat, pivot_cols)
return mat