from types import FunctionType from sympy.core.numbers import Float, Integer from sympy.core.singleton import S from sympy.core.symbol import uniquely_named_symbol from sympy.core.mul import Mul from sympy.polys import PurePoly, cancel from sympy.functions.combinatorial.numbers import nC from sympy.polys.matrices.domainmatrix import DomainMatrix from .common import NonSquareMatrixError from .utilities import ( _get_intermediate_simp, _get_intermediate_simp_bool, _iszero, _is_zero_after_expand_mul, _dotprodsimp, _simplify) def _find_reasonable_pivot(col, iszerofunc=_iszero, simpfunc=_simplify): """ Find the lowest index of an item in ``col`` that is suitable for a pivot. If ``col`` consists only of Floats, the pivot with the largest norm is returned. Otherwise, the first element where ``iszerofunc`` returns False is used. If ``iszerofunc`` does not return false, items are simplified and retested until a suitable pivot is found. Returns a 4-tuple (pivot_offset, pivot_val, assumed_nonzero, newly_determined) where pivot_offset is the index of the pivot, pivot_val is the (possibly simplified) value of the pivot, assumed_nonzero is True if an assumption that the pivot was non-zero was made without being proved, and newly_determined are elements that were simplified during the process of pivot finding.""" newly_determined = [] col = list(col) # a column that contains a mix of floats and integers # but at least one float is considered a numerical # column, and so we do partial pivoting if all(isinstance(x, (Float, Integer)) for x in col) and any( isinstance(x, Float) for x in col): col_abs = [abs(x) for x in col] max_value = max(col_abs) if iszerofunc(max_value): # just because iszerofunc returned True, doesn't # mean the value is numerically zero. Make sure # to replace all entries with numerical zeros if max_value != 0: newly_determined = [(i, 0) for i, x in enumerate(col) if x != 0] return (None, None, False, newly_determined) index = col_abs.index(max_value) return (index, col[index], False, newly_determined) # PASS 1 (iszerofunc directly) possible_zeros = [] for i, x in enumerate(col): is_zero = iszerofunc(x) # is someone wrote a custom iszerofunc, it may return # BooleanFalse or BooleanTrue instead of True or False, # so use == for comparison instead of `is` if is_zero == False: # we found something that is definitely not zero return (i, x, False, newly_determined) possible_zeros.append(is_zero) # by this point, we've found no certain non-zeros if all(possible_zeros): # if everything is definitely zero, we have # no pivot return (None, None, False, newly_determined) # PASS 2 (iszerofunc after simplify) # we haven't found any for-sure non-zeros, so # go through the elements iszerofunc couldn't # make a determination about and opportunistically # simplify to see if we find something for i, x in enumerate(col): if possible_zeros[i] is not None: continue simped = simpfunc(x) is_zero = iszerofunc(simped) if is_zero in (True, False): newly_determined.append((i, simped)) if is_zero == False: return (i, simped, False, newly_determined) possible_zeros[i] = is_zero # after simplifying, some things that were recognized # as zeros might be zeros if all(possible_zeros): # if everything is definitely zero, we have # no pivot return (None, None, False, newly_determined) # PASS 3 (.equals(0)) # some expressions fail to simplify to zero, but # ``.equals(0)`` evaluates to True. As a last-ditch # attempt, apply ``.equals`` to these expressions for i, x in enumerate(col): if possible_zeros[i] is not None: continue if x.equals(S.Zero): # ``.iszero`` may return False with # an implicit assumption (e.g., ``x.equals(0)`` # when ``x`` is a symbol), so only treat it # as proved when ``.equals(0)`` returns True possible_zeros[i] = True newly_determined.append((i, S.Zero)) if all(possible_zeros): return (None, None, False, newly_determined) # at this point there is nothing that could definitely # be a pivot. To maintain compatibility with existing # behavior, we'll assume that an illdetermined thing is # non-zero. We should probably raise a warning in this case i = possible_zeros.index(None) return (i, col[i], True, newly_determined) def _find_reasonable_pivot_naive(col, iszerofunc=_iszero, simpfunc=None): """ Helper that computes the pivot value and location from a sequence of contiguous matrix column elements. As a side effect of the pivot search, this function may simplify some of the elements of the input column. A list of these simplified entries and their indices are also returned. This function mimics the behavior of _find_reasonable_pivot(), but does less work trying to determine if an indeterminate candidate pivot simplifies to zero. This more naive approach can be much faster, with the trade-off that it may erroneously return a pivot that is zero. ``col`` is a sequence of contiguous column entries to be searched for a suitable pivot. ``iszerofunc`` is a callable that returns a Boolean that indicates if its input is zero, or None if no such determination can be made. ``simpfunc`` is a callable that simplifies its input. It must return its input if it does not simplify its input. Passing in ``simpfunc=None`` indicates that the pivot search should not attempt to simplify any candidate pivots. Returns a 4-tuple: (pivot_offset, pivot_val, assumed_nonzero, newly_determined) ``pivot_offset`` is the sequence index of the pivot. ``pivot_val`` is the value of the pivot. pivot_val and col[pivot_index] are equivalent, but will be different when col[pivot_index] was simplified during the pivot search. ``assumed_nonzero`` is a boolean indicating if the pivot cannot be guaranteed to be zero. If assumed_nonzero is true, then the pivot may or may not be non-zero. If assumed_nonzero is false, then the pivot is non-zero. ``newly_determined`` is a list of index-value pairs of pivot candidates that were simplified during the pivot search. """ # indeterminates holds the index-value pairs of each pivot candidate # that is neither zero or non-zero, as determined by iszerofunc(). # If iszerofunc() indicates that a candidate pivot is guaranteed # non-zero, or that every candidate pivot is zero then the contents # of indeterminates are unused. # Otherwise, the only viable candidate pivots are symbolic. # In this case, indeterminates will have at least one entry, # and all but the first entry are ignored when simpfunc is None. indeterminates = [] for i, col_val in enumerate(col): col_val_is_zero = iszerofunc(col_val) if col_val_is_zero == False: # This pivot candidate is non-zero. return i, col_val, False, [] elif col_val_is_zero is None: # The candidate pivot's comparison with zero # is indeterminate. indeterminates.append((i, col_val)) if len(indeterminates) == 0: # All candidate pivots are guaranteed to be zero, i.e. there is # no pivot. return None, None, False, [] if simpfunc is None: # Caller did not pass in a simplification function that might # determine if an indeterminate pivot candidate is guaranteed # to be nonzero, so assume the first indeterminate candidate # is non-zero. return indeterminates[0][0], indeterminates[0][1], True, [] # newly_determined holds index-value pairs of candidate pivots # that were simplified during the search for a non-zero pivot. newly_determined = [] for i, col_val in indeterminates: tmp_col_val = simpfunc(col_val) if id(col_val) != id(tmp_col_val): # simpfunc() simplified this candidate pivot. newly_determined.append((i, tmp_col_val)) if iszerofunc(tmp_col_val) == False: # Candidate pivot simplified to a guaranteed non-zero value. return i, tmp_col_val, False, newly_determined return indeterminates[0][0], indeterminates[0][1], True, newly_determined # This functions is a candidate for caching if it gets implemented for matrices. def _berkowitz_toeplitz_matrix(M): """Return (A,T) where T the Toeplitz matrix used in the Berkowitz algorithm corresponding to ``M`` and A is the first principal submatrix. """ # the 0 x 0 case is trivial if M.rows == 0 and M.cols == 0: return M._new(1,1, [M.one]) # # Partition M = [ a_11 R ] # [ C A ] # a, R = M[0,0], M[0, 1:] C, A = M[1:, 0], M[1:,1:] # # The Toeplitz matrix looks like # # [ 1 ] # [ -a 1 ] # [ -RC -a 1 ] # [ -RAC -RC -a 1 ] # [ -RA**2C -RAC -RC -a 1 ] # etc. # Compute the diagonal entries. # Because multiplying matrix times vector is so much # more efficient than matrix times matrix, recursively # compute -R * A**n * C. diags = [C] for i in range(M.rows - 2): diags.append(A.multiply(diags[i], dotprodsimp=None)) diags = [(-R).multiply(d, dotprodsimp=None)[0, 0] for d in diags] diags = [M.one, -a] + diags def entry(i,j): if j > i: return M.zero return diags[i - j] toeplitz = M._new(M.cols + 1, M.rows, entry) return (A, toeplitz) # This functions is a candidate for caching if it gets implemented for matrices. def _berkowitz_vector(M): """ Run the Berkowitz algorithm and return a vector whose entries are the coefficients of the characteristic polynomial of ``M``. Given N x N matrix, efficiently compute coefficients of characteristic polynomials of ``M`` without division in the ground domain. This method is particularly useful for computing determinant, principal minors and characteristic polynomial when ``M`` has complicated coefficients e.g. polynomials. Semi-direct usage of this algorithm is also important in computing efficiently sub-resultant PRS. Assuming that M is a square matrix of dimension N x N and I is N x N identity matrix, then the Berkowitz vector is an N x 1 vector whose entries are coefficients of the polynomial charpoly(M) = det(t*I - M) As a consequence, all polynomials generated by Berkowitz algorithm are monic. For more information on the implemented algorithm refer to: [1] S.J. Berkowitz, On computing the determinant in small parallel time using a small number of processors, ACM, Information Processing Letters 18, 1984, pp. 147-150 [2] M. Keber, Division-Free computation of sub-resultants using Bezout matrices, Tech. Report MPI-I-2006-1-006, Saarbrucken, 2006 """ # handle the trivial cases if M.rows == 0 and M.cols == 0: return M._new(1, 1, [M.one]) elif M.rows == 1 and M.cols == 1: return M._new(2, 1, [M.one, -M[0,0]]) submat, toeplitz = _berkowitz_toeplitz_matrix(M) return toeplitz.multiply(_berkowitz_vector(submat), dotprodsimp=None) def _adjugate(M, method="berkowitz"): """Returns the adjugate, or classical adjoint, of a matrix. That is, the transpose of the matrix of cofactors. https://en.wikipedia.org/wiki/Adjugate Parameters ========== method : string, optional Method to use to find the cofactors, can be "bareiss", "berkowitz" or "lu". Examples ======== >>> from sympy import Matrix >>> M = Matrix([[1, 2], [3, 4]]) >>> M.adjugate() Matrix([ [ 4, -2], [-3, 1]]) See Also ======== cofactor_matrix sympy.matrices.common.MatrixCommon.transpose """ return M.cofactor_matrix(method=method).transpose() # This functions is a candidate for caching if it gets implemented for matrices. def _charpoly(M, x='lambda', simplify=_simplify): """Computes characteristic polynomial det(x*I - M) where I is the identity matrix. A PurePoly is returned, so using different variables for ``x`` does not affect the comparison or the polynomials: Parameters ========== x : string, optional Name for the "lambda" variable, defaults to "lambda". simplify : function, optional Simplification function to use on the characteristic polynomial calculated. Defaults to ``simplify``. Examples ======== >>> from sympy import Matrix >>> from sympy.abc import x, y >>> M = Matrix([[1, 3], [2, 0]]) >>> M.charpoly() PurePoly(lambda**2 - lambda - 6, lambda, domain='ZZ') >>> M.charpoly(x) == M.charpoly(y) True >>> M.charpoly(x) == M.charpoly(y) True Specifying ``x`` is optional; a symbol named ``lambda`` is used by default (which looks good when pretty-printed in unicode): >>> M.charpoly().as_expr() lambda**2 - lambda - 6 And if ``x`` clashes with an existing symbol, underscores will be prepended to the name to make it unique: >>> M = Matrix([[1, 2], [x, 0]]) >>> M.charpoly(x).as_expr() _x**2 - _x - 2*x Whether you pass a symbol or not, the generator can be obtained with the gen attribute since it may not be the same as the symbol that was passed: >>> M.charpoly(x).gen _x >>> M.charpoly(x).gen == x False Notes ===== The Samuelson-Berkowitz algorithm is used to compute the characteristic polynomial efficiently and without any division operations. Thus the characteristic polynomial over any commutative ring without zero divisors can be computed. If the determinant det(x*I - M) can be found out easily as in the case of an upper or a lower triangular matrix, then instead of Samuelson-Berkowitz algorithm, eigenvalues are computed and the characteristic polynomial with their help. See Also ======== det """ if not M.is_square: raise NonSquareMatrixError() if M.is_lower or M.is_upper: diagonal_elements = M.diagonal() x = uniquely_named_symbol(x, diagonal_elements, modify=lambda s: '_' + s) m = 1 for i in diagonal_elements: m = m * (x - simplify(i)) return PurePoly(m, x) berk_vector = _berkowitz_vector(M) x = uniquely_named_symbol(x, berk_vector, modify=lambda s: '_' + s) return PurePoly([simplify(a) for a in berk_vector], x) def _cofactor(M, i, j, method="berkowitz"): """Calculate the cofactor of an element. Parameters ========== method : string, optional Method to use to find the cofactors, can be "bareiss", "berkowitz" or "lu". Examples ======== >>> from sympy import Matrix >>> M = Matrix([[1, 2], [3, 4]]) >>> M.cofactor(0, 1) -3 See Also ======== cofactor_matrix minor minor_submatrix """ if not M.is_square or M.rows < 1: raise NonSquareMatrixError() return S.NegativeOne**((i + j) % 2) * M.minor(i, j, method) def _cofactor_matrix(M, method="berkowitz"): """Return a matrix containing the cofactor of each element. Parameters ========== method : string, optional Method to use to find the cofactors, can be "bareiss", "berkowitz" or "lu". Examples ======== >>> from sympy import Matrix >>> M = Matrix([[1, 2], [3, 4]]) >>> M.cofactor_matrix() Matrix([ [ 4, -3], [-2, 1]]) See Also ======== cofactor minor minor_submatrix """ if not M.is_square or M.rows < 1: raise NonSquareMatrixError() return M._new(M.rows, M.cols, lambda i, j: M.cofactor(i, j, method)) def _per(M): """Returns the permanent of a matrix. Unlike determinant, permanent is defined for both square and non-square matrices. For an m x n matrix, with m less than or equal to n, it is given as the sum over the permutations s of size less than or equal to m on [1, 2, . . . n] of the product from i = 1 to m of M[i, s[i]]. Taking the transpose will not affect the value of the permanent. In the case of a square matrix, this is the same as the permutation definition of the determinant, but it does not take the sign of the permutation into account. Computing the permanent with this definition is quite inefficient, so here the Ryser formula is used. Examples ======== >>> from sympy import Matrix >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> M.per() 450 >>> M = Matrix([1, 5, 7]) >>> M.per() 13 References ========== .. [1] Prof. Frank Ben's notes: https://math.berkeley.edu/~bernd/ban275.pdf .. [2] Wikipedia article on Permanent: https://en.wikipedia.org/wiki/Permanent_%28mathematics%29 .. [3] https://reference.wolfram.com/language/ref/Permanent.html .. [4] Permanent of a rectangular matrix : https://arxiv.org/pdf/0904.3251.pdf """ import itertools m, n = M.shape if m > n: M = M.T m, n = n, m s = list(range(n)) subsets = [] for i in range(1, m + 1): subsets += list(map(list, itertools.combinations(s, i))) perm = 0 for subset in subsets: prod = 1 sub_len = len(subset) for i in range(m): prod *= sum([M[i, j] for j in subset]) perm += prod * S.NegativeOne**sub_len * nC(n - sub_len, m - sub_len) perm *= S.NegativeOne**m return perm.simplify() def _det_DOM(M): DOM = DomainMatrix.from_Matrix(M, field=True, extension=True) K = DOM.domain return K.to_sympy(DOM.det()) # This functions is a candidate for caching if it gets implemented for matrices. def _det(M, method="bareiss", iszerofunc=None): """Computes the determinant of a matrix if ``M`` is a concrete matrix object otherwise return an expressions ``Determinant(M)`` if ``M`` is a ``MatrixSymbol`` or other expression. Parameters ========== method : string, optional Specifies the algorithm used for computing the matrix determinant. If the matrix is at most 3x3, a hard-coded formula is used and the specified method is ignored. Otherwise, it defaults to ``'bareiss'``. Also, if the matrix is an upper or a lower triangular matrix, determinant is computed by simple multiplication of diagonal elements, and the specified method is ignored. If it is set to ``'domain-ge'``, then Gaussian elimination method will be used via using DomainMatrix. If it is set to ``'bareiss'``, Bareiss' fraction-free algorithm will be used. If it is set to ``'berkowitz'``, Berkowitz' algorithm will be used. Otherwise, if it is set to ``'lu'``, LU decomposition will be used. .. note:: For backward compatibility, legacy keys like "bareis" and "det_lu" can still be used to indicate the corresponding methods. And the keys are also case-insensitive for now. However, it is suggested to use the precise keys for specifying the method. iszerofunc : FunctionType or None, optional If it is set to ``None``, it will be defaulted to ``_iszero`` if the method is set to ``'bareiss'``, and ``_is_zero_after_expand_mul`` if the method is set to ``'lu'``. It can also accept any user-specified zero testing function, if it is formatted as a function which accepts a single symbolic argument and returns ``True`` if it is tested as zero and ``False`` if it tested as non-zero, and also ``None`` if it is undecidable. Returns ======= det : Basic Result of determinant. Raises ====== ValueError If unrecognized keys are given for ``method`` or ``iszerofunc``. NonSquareMatrixError If attempted to calculate determinant from a non-square matrix. Examples ======== >>> from sympy import Matrix, eye, det >>> I3 = eye(3) >>> det(I3) 1 >>> M = Matrix([[1, 2], [3, 4]]) >>> det(M) -2 >>> det(M) == M.det() True >>> M.det(method="domain-ge") -2 """ # sanitize `method` method = method.lower() if method == "bareis": method = "bareiss" elif method == "det_lu": method = "lu" if method not in ("bareiss", "berkowitz", "lu", "domain-ge"): raise ValueError("Determinant method '%s' unrecognized" % method) if iszerofunc is None: if method == "bareiss": iszerofunc = _is_zero_after_expand_mul elif method == "lu": iszerofunc = _iszero elif not isinstance(iszerofunc, FunctionType): raise ValueError("Zero testing method '%s' unrecognized" % iszerofunc) n = M.rows if n == M.cols: # square check is done in individual method functions if n == 0: return M.one elif n == 1: return M[0, 0] elif n == 2: m = M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0] return _get_intermediate_simp(_dotprodsimp)(m) elif n == 3: m = (M[0, 0] * M[1, 1] * M[2, 2] + M[0, 1] * M[1, 2] * M[2, 0] + M[0, 2] * M[1, 0] * M[2, 1] - M[0, 2] * M[1, 1] * M[2, 0] - M[0, 0] * M[1, 2] * M[2, 1] - M[0, 1] * M[1, 0] * M[2, 2]) return _get_intermediate_simp(_dotprodsimp)(m) dets = [] for b in M.strongly_connected_components(): if method == "domain-ge": # uses DomainMatrix to evaluate determinant det = _det_DOM(M[b, b]) elif method == "bareiss": det = M[b, b]._eval_det_bareiss(iszerofunc=iszerofunc) elif method == "berkowitz": det = M[b, b]._eval_det_berkowitz() elif method == "lu": det = M[b, b]._eval_det_lu(iszerofunc=iszerofunc) dets.append(det) return Mul(*dets) # This functions is a candidate for caching if it gets implemented for matrices. def _det_bareiss(M, iszerofunc=_is_zero_after_expand_mul): """Compute matrix determinant using Bareiss' fraction-free algorithm which is an extension of the well known Gaussian elimination method. This approach is best suited for dense symbolic matrices and will result in a determinant with minimal number of fractions. It means that less term rewriting is needed on resulting formulae. Parameters ========== iszerofunc : function, optional The function to use to determine zeros when doing an LU decomposition. Defaults to ``lambda x: x.is_zero``. TODO: Implement algorithm for sparse matrices (SFF), http://www.eecis.udel.edu/~saunders/papers/sffge/it5.ps. """ # Recursively implemented Bareiss' algorithm as per Deanna Richelle Leggett's # thesis http://www.math.usm.edu/perry/Research/Thesis_DRL.pdf def bareiss(mat, cumm=1): if mat.rows == 0: return mat.one elif mat.rows == 1: return mat[0, 0] # find a pivot and extract the remaining matrix # With the default iszerofunc, _find_reasonable_pivot slows down # the computation by the factor of 2.5 in one test. # Relevant issues: #10279 and #13877. pivot_pos, pivot_val, _, _ = _find_reasonable_pivot(mat[:, 0], iszerofunc=iszerofunc) if pivot_pos is None: return mat.zero # if we have a valid pivot, we'll do a "row swap", so keep the # sign of the det sign = (-1) ** (pivot_pos % 2) # we want every row but the pivot row and every column rows = [i for i in range(mat.rows) if i != pivot_pos] cols = list(range(mat.cols)) tmp_mat = mat.extract(rows, cols) def entry(i, j): ret = (pivot_val*tmp_mat[i, j + 1] - mat[pivot_pos, j + 1]*tmp_mat[i, 0]) / cumm if _get_intermediate_simp_bool(True): return _dotprodsimp(ret) elif not ret.is_Atom: return cancel(ret) return ret return sign*bareiss(M._new(mat.rows - 1, mat.cols - 1, entry), pivot_val) if not M.is_square: raise NonSquareMatrixError() if M.rows == 0: return M.one # sympy/matrices/tests/test_matrices.py contains a test that # suggests that the determinant of a 0 x 0 matrix is one, by # convention. return bareiss(M) def _det_berkowitz(M): """ Use the Berkowitz algorithm to compute the determinant.""" if not M.is_square: raise NonSquareMatrixError() if M.rows == 0: return M.one # sympy/matrices/tests/test_matrices.py contains a test that # suggests that the determinant of a 0 x 0 matrix is one, by # convention. berk_vector = _berkowitz_vector(M) return (-1)**(len(berk_vector) - 1) * berk_vector[-1] # This functions is a candidate for caching if it gets implemented for matrices. def _det_LU(M, iszerofunc=_iszero, simpfunc=None): """ Computes the determinant of a matrix from its LU decomposition. This function uses the LU decomposition computed by LUDecomposition_Simple(). The keyword arguments iszerofunc and simpfunc are passed to LUDecomposition_Simple(). iszerofunc is a callable that returns a boolean indicating if its input is zero, or None if it cannot make the determination. simpfunc is a callable that simplifies its input. The default is simpfunc=None, which indicate that the pivot search algorithm should not attempt to simplify any candidate pivots. If simpfunc fails to simplify its input, then it must return its input instead of a copy. Parameters ========== iszerofunc : function, optional The function to use to determine zeros when doing an LU decomposition. Defaults to ``lambda x: x.is_zero``. simpfunc : function, optional The simplification function to use when looking for zeros for pivots. """ if not M.is_square: raise NonSquareMatrixError() if M.rows == 0: return M.one # sympy/matrices/tests/test_matrices.py contains a test that # suggests that the determinant of a 0 x 0 matrix is one, by # convention. lu, row_swaps = M.LUdecomposition_Simple(iszerofunc=iszerofunc, simpfunc=simpfunc) # P*A = L*U => det(A) = det(L)*det(U)/det(P) = det(P)*det(U). # Lower triangular factor L encoded in lu has unit diagonal => det(L) = 1. # P is a permutation matrix => det(P) in {-1, 1} => 1/det(P) = det(P). # LUdecomposition_Simple() returns a list of row exchange index pairs, rather # than a permutation matrix, but det(P) = (-1)**len(row_swaps). # Avoid forming the potentially time consuming product of U's diagonal entries # if the product is zero. # Bottom right entry of U is 0 => det(A) = 0. # It may be impossible to determine if this entry of U is zero when it is symbolic. if iszerofunc(lu[lu.rows-1, lu.rows-1]): return M.zero # Compute det(P) det = -M.one if len(row_swaps)%2 else M.one # Compute det(U) by calculating the product of U's diagonal entries. # The upper triangular portion of lu is the upper triangular portion of the # U factor in the LU decomposition. for k in range(lu.rows): det *= lu[k, k] # return det(P)*det(U) return det def _minor(M, i, j, method="berkowitz"): """Return the (i,j) minor of ``M``. That is, return the determinant of the matrix obtained by deleting the `i`th row and `j`th column from ``M``. Parameters ========== i, j : int The row and column to exclude to obtain the submatrix. method : string, optional Method to use to find the determinant of the submatrix, can be "bareiss", "berkowitz" or "lu". Examples ======== >>> from sympy import Matrix >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> M.minor(1, 1) -12 See Also ======== minor_submatrix cofactor det """ if not M.is_square: raise NonSquareMatrixError() return M.minor_submatrix(i, j).det(method=method) def _minor_submatrix(M, i, j): """Return the submatrix obtained by removing the `i`th row and `j`th column from ``M`` (works with Pythonic negative indices). Parameters ========== i, j : int The row and column to exclude to obtain the submatrix. Examples ======== >>> from sympy import Matrix >>> M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) >>> M.minor_submatrix(1, 1) Matrix([ [1, 3], [7, 9]]) See Also ======== minor cofactor """ if i < 0: i += M.rows if j < 0: j += M.cols if not 0 <= i < M.rows or not 0 <= j < M.cols: raise ValueError("`i` and `j` must satisfy 0 <= i < ``M.rows`` " "(%d)" % M.rows + "and 0 <= j < ``M.cols`` (%d)." % M.cols) rows = [a for a in range(M.rows) if a != i] cols = [a for a in range(M.cols) if a != j] return M.extract(rows, cols)