Coverage for /usr/lib/python3/dist-packages/sympy/matrices/common.py: 26%
1012 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2Basic methods common to all matrices to be used
3when creating more advanced matrices (e.g., matrices over rings,
4etc.).
5"""
7from collections import defaultdict
8from collections.abc import Iterable
9from inspect import isfunction
10from functools import reduce
12from sympy.assumptions.refine import refine
13from sympy.core import SympifyError, Add
14from sympy.core.basic import Atom
15from sympy.core.decorators import call_highest_priority
16from sympy.core.kind import Kind, NumberKind
17from sympy.core.logic import fuzzy_and, FuzzyBool
18from sympy.core.mod import Mod
19from sympy.core.singleton import S
20from sympy.core.symbol import Symbol
21from sympy.core.sympify import sympify
22from sympy.functions.elementary.complexes import Abs, re, im
23from .utilities import _dotprodsimp, _simplify
24from sympy.polys.polytools import Poly
25from sympy.utilities.iterables import flatten, is_sequence
26from sympy.utilities.misc import as_int, filldedent
27from sympy.tensor.array import NDimArray
29from .utilities import _get_intermediate_simp_bool
32class MatrixError(Exception):
33 pass
36class ShapeError(ValueError, MatrixError):
37 """Wrong matrix shape"""
38 pass
41class NonSquareMatrixError(ShapeError):
42 pass
45class NonInvertibleMatrixError(ValueError, MatrixError):
46 """The matrix in not invertible (division by multidimensional zero error)."""
47 pass
50class NonPositiveDefiniteMatrixError(ValueError, MatrixError):
51 """The matrix is not a positive-definite matrix."""
52 pass
55class MatrixRequired:
56 """All subclasses of matrix objects must implement the
57 required matrix properties listed here."""
58 rows = None # type: int
59 cols = None # type: int
60 _simplify = None
62 @classmethod
63 def _new(cls, *args, **kwargs):
64 """`_new` must, at minimum, be callable as
65 `_new(rows, cols, mat) where mat is a flat list of the
66 elements of the matrix."""
67 raise NotImplementedError("Subclasses must implement this.")
69 def __eq__(self, other):
70 raise NotImplementedError("Subclasses must implement this.")
72 def __getitem__(self, key):
73 """Implementations of __getitem__ should accept ints, in which
74 case the matrix is indexed as a flat list, tuples (i,j) in which
75 case the (i,j) entry is returned, slices, or mixed tuples (a,b)
76 where a and b are any combination of slices and integers."""
77 raise NotImplementedError("Subclasses must implement this.")
79 def __len__(self):
80 """The total number of entries in the matrix."""
81 raise NotImplementedError("Subclasses must implement this.")
83 @property
84 def shape(self):
85 raise NotImplementedError("Subclasses must implement this.")
88class MatrixShaping(MatrixRequired):
89 """Provides basic matrix shaping and extracting of submatrices"""
91 def _eval_col_del(self, col):
92 def entry(i, j):
93 return self[i, j] if j < col else self[i, j + 1]
94 return self._new(self.rows, self.cols - 1, entry)
96 def _eval_col_insert(self, pos, other):
98 def entry(i, j):
99 if j < pos:
100 return self[i, j]
101 elif pos <= j < pos + other.cols:
102 return other[i, j - pos]
103 return self[i, j - other.cols]
105 return self._new(self.rows, self.cols + other.cols, entry)
107 def _eval_col_join(self, other):
108 rows = self.rows
110 def entry(i, j):
111 if i < rows:
112 return self[i, j]
113 return other[i - rows, j]
115 return classof(self, other)._new(self.rows + other.rows, self.cols,
116 entry)
118 def _eval_extract(self, rowsList, colsList):
119 mat = list(self)
120 cols = self.cols
121 indices = (i * cols + j for i in rowsList for j in colsList)
122 return self._new(len(rowsList), len(colsList),
123 [mat[i] for i in indices])
125 def _eval_get_diag_blocks(self):
126 sub_blocks = []
128 def recurse_sub_blocks(M):
129 i = 1
130 while i <= M.shape[0]:
131 if i == 1:
132 to_the_right = M[0, i:]
133 to_the_bottom = M[i:, 0]
134 else:
135 to_the_right = M[:i, i:]
136 to_the_bottom = M[i:, :i]
137 if any(to_the_right) or any(to_the_bottom):
138 i += 1
139 continue
140 else:
141 sub_blocks.append(M[:i, :i])
142 if M.shape == M[:i, :i].shape:
143 return
144 else:
145 recurse_sub_blocks(M[i:, i:])
146 return
148 recurse_sub_blocks(self)
149 return sub_blocks
151 def _eval_row_del(self, row):
152 def entry(i, j):
153 return self[i, j] if i < row else self[i + 1, j]
154 return self._new(self.rows - 1, self.cols, entry)
156 def _eval_row_insert(self, pos, other):
157 entries = list(self)
158 insert_pos = pos * self.cols
159 entries[insert_pos:insert_pos] = list(other)
160 return self._new(self.rows + other.rows, self.cols, entries)
162 def _eval_row_join(self, other):
163 cols = self.cols
165 def entry(i, j):
166 if j < cols:
167 return self[i, j]
168 return other[i, j - cols]
170 return classof(self, other)._new(self.rows, self.cols + other.cols,
171 entry)
173 def _eval_tolist(self):
174 return [list(self[i,:]) for i in range(self.rows)]
176 def _eval_todok(self):
177 dok = {}
178 rows, cols = self.shape
179 for i in range(rows):
180 for j in range(cols):
181 val = self[i, j]
182 if val != self.zero:
183 dok[i, j] = val
184 return dok
186 def _eval_vec(self):
187 rows = self.rows
189 def entry(n, _):
190 # we want to read off the columns first
191 j = n // rows
192 i = n - j * rows
193 return self[i, j]
195 return self._new(len(self), 1, entry)
197 def _eval_vech(self, diagonal):
198 c = self.cols
199 v = []
200 if diagonal:
201 for j in range(c):
202 for i in range(j, c):
203 v.append(self[i, j])
204 else:
205 for j in range(c):
206 for i in range(j + 1, c):
207 v.append(self[i, j])
208 return self._new(len(v), 1, v)
210 def col_del(self, col):
211 """Delete the specified column."""
212 if col < 0:
213 col += self.cols
214 if not 0 <= col < self.cols:
215 raise IndexError("Column {} is out of range.".format(col))
216 return self._eval_col_del(col)
218 def col_insert(self, pos, other):
219 """Insert one or more columns at the given column position.
221 Examples
222 ========
224 >>> from sympy import zeros, ones
225 >>> M = zeros(3)
226 >>> V = ones(3, 1)
227 >>> M.col_insert(1, V)
228 Matrix([
229 [0, 1, 0, 0],
230 [0, 1, 0, 0],
231 [0, 1, 0, 0]])
233 See Also
234 ========
236 col
237 row_insert
238 """
239 # Allows you to build a matrix even if it is null matrix
240 if not self:
241 return type(self)(other)
243 pos = as_int(pos)
245 if pos < 0:
246 pos = self.cols + pos
247 if pos < 0:
248 pos = 0
249 elif pos > self.cols:
250 pos = self.cols
252 if self.rows != other.rows:
253 raise ShapeError(
254 "The matrices have incompatible number of rows ({} and {})"
255 .format(self.rows, other.rows))
257 return self._eval_col_insert(pos, other)
259 def col_join(self, other):
260 """Concatenates two matrices along self's last and other's first row.
262 Examples
263 ========
265 >>> from sympy import zeros, ones
266 >>> M = zeros(3)
267 >>> V = ones(1, 3)
268 >>> M.col_join(V)
269 Matrix([
270 [0, 0, 0],
271 [0, 0, 0],
272 [0, 0, 0],
273 [1, 1, 1]])
275 See Also
276 ========
278 col
279 row_join
280 """
281 # A null matrix can always be stacked (see #10770)
282 if self.rows == 0 and self.cols != other.cols:
283 return self._new(0, other.cols, []).col_join(other)
285 if self.cols != other.cols:
286 raise ShapeError(
287 "The matrices have incompatible number of columns ({} and {})"
288 .format(self.cols, other.cols))
289 return self._eval_col_join(other)
291 def col(self, j):
292 """Elementary column selector.
294 Examples
295 ========
297 >>> from sympy import eye
298 >>> eye(2).col(0)
299 Matrix([
300 [1],
301 [0]])
303 See Also
304 ========
306 row
307 col_del
308 col_join
309 col_insert
310 """
311 return self[:, j]
313 def extract(self, rowsList, colsList):
314 r"""Return a submatrix by specifying a list of rows and columns.
315 Negative indices can be given. All indices must be in the range
316 $-n \le i < n$ where $n$ is the number of rows or columns.
318 Examples
319 ========
321 >>> from sympy import Matrix
322 >>> m = Matrix(4, 3, range(12))
323 >>> m
324 Matrix([
325 [0, 1, 2],
326 [3, 4, 5],
327 [6, 7, 8],
328 [9, 10, 11]])
329 >>> m.extract([0, 1, 3], [0, 1])
330 Matrix([
331 [0, 1],
332 [3, 4],
333 [9, 10]])
335 Rows or columns can be repeated:
337 >>> m.extract([0, 0, 1], [-1])
338 Matrix([
339 [2],
340 [2],
341 [5]])
343 Every other row can be taken by using range to provide the indices:
345 >>> m.extract(range(0, m.rows, 2), [-1])
346 Matrix([
347 [2],
348 [8]])
350 RowsList or colsList can also be a list of booleans, in which case
351 the rows or columns corresponding to the True values will be selected:
353 >>> m.extract([0, 1, 2, 3], [True, False, True])
354 Matrix([
355 [0, 2],
356 [3, 5],
357 [6, 8],
358 [9, 11]])
359 """
361 if not is_sequence(rowsList) or not is_sequence(colsList):
362 raise TypeError("rowsList and colsList must be iterable")
363 # ensure rowsList and colsList are lists of integers
364 if rowsList and all(isinstance(i, bool) for i in rowsList):
365 rowsList = [index for index, item in enumerate(rowsList) if item]
366 if colsList and all(isinstance(i, bool) for i in colsList):
367 colsList = [index for index, item in enumerate(colsList) if item]
369 # ensure everything is in range
370 rowsList = [a2idx(k, self.rows) for k in rowsList]
371 colsList = [a2idx(k, self.cols) for k in colsList]
373 return self._eval_extract(rowsList, colsList)
375 def get_diag_blocks(self):
376 """Obtains the square sub-matrices on the main diagonal of a square matrix.
378 Useful for inverting symbolic matrices or solving systems of
379 linear equations which may be decoupled by having a block diagonal
380 structure.
382 Examples
383 ========
385 >>> from sympy import Matrix
386 >>> from sympy.abc import x, y, z
387 >>> A = Matrix([[1, 3, 0, 0], [y, z*z, 0, 0], [0, 0, x, 0], [0, 0, 0, 0]])
388 >>> a1, a2, a3 = A.get_diag_blocks()
389 >>> a1
390 Matrix([
391 [1, 3],
392 [y, z**2]])
393 >>> a2
394 Matrix([[x]])
395 >>> a3
396 Matrix([[0]])
398 """
399 return self._eval_get_diag_blocks()
401 @classmethod
402 def hstack(cls, *args):
403 """Return a matrix formed by joining args horizontally (i.e.
404 by repeated application of row_join).
406 Examples
407 ========
409 >>> from sympy import Matrix, eye
410 >>> Matrix.hstack(eye(2), 2*eye(2))
411 Matrix([
412 [1, 0, 2, 0],
413 [0, 1, 0, 2]])
414 """
415 if len(args) == 0:
416 return cls._new()
418 kls = type(args[0])
419 return reduce(kls.row_join, args)
421 def reshape(self, rows, cols):
422 """Reshape the matrix. Total number of elements must remain the same.
424 Examples
425 ========
427 >>> from sympy import Matrix
428 >>> m = Matrix(2, 3, lambda i, j: 1)
429 >>> m
430 Matrix([
431 [1, 1, 1],
432 [1, 1, 1]])
433 >>> m.reshape(1, 6)
434 Matrix([[1, 1, 1, 1, 1, 1]])
435 >>> m.reshape(3, 2)
436 Matrix([
437 [1, 1],
438 [1, 1],
439 [1, 1]])
441 """
442 if self.rows * self.cols != rows * cols:
443 raise ValueError("Invalid reshape parameters %d %d" % (rows, cols))
444 return self._new(rows, cols, lambda i, j: self[i * cols + j])
446 def row_del(self, row):
447 """Delete the specified row."""
448 if row < 0:
449 row += self.rows
450 if not 0 <= row < self.rows:
451 raise IndexError("Row {} is out of range.".format(row))
453 return self._eval_row_del(row)
455 def row_insert(self, pos, other):
456 """Insert one or more rows at the given row position.
458 Examples
459 ========
461 >>> from sympy import zeros, ones
462 >>> M = zeros(3)
463 >>> V = ones(1, 3)
464 >>> M.row_insert(1, V)
465 Matrix([
466 [0, 0, 0],
467 [1, 1, 1],
468 [0, 0, 0],
469 [0, 0, 0]])
471 See Also
472 ========
474 row
475 col_insert
476 """
477 # Allows you to build a matrix even if it is null matrix
478 if not self:
479 return self._new(other)
481 pos = as_int(pos)
483 if pos < 0:
484 pos = self.rows + pos
485 if pos < 0:
486 pos = 0
487 elif pos > self.rows:
488 pos = self.rows
490 if self.cols != other.cols:
491 raise ShapeError(
492 "The matrices have incompatible number of columns ({} and {})"
493 .format(self.cols, other.cols))
495 return self._eval_row_insert(pos, other)
497 def row_join(self, other):
498 """Concatenates two matrices along self's last and rhs's first column
500 Examples
501 ========
503 >>> from sympy import zeros, ones
504 >>> M = zeros(3)
505 >>> V = ones(3, 1)
506 >>> M.row_join(V)
507 Matrix([
508 [0, 0, 0, 1],
509 [0, 0, 0, 1],
510 [0, 0, 0, 1]])
512 See Also
513 ========
515 row
516 col_join
517 """
518 # A null matrix can always be stacked (see #10770)
519 if self.cols == 0 and self.rows != other.rows:
520 return self._new(other.rows, 0, []).row_join(other)
522 if self.rows != other.rows:
523 raise ShapeError(
524 "The matrices have incompatible number of rows ({} and {})"
525 .format(self.rows, other.rows))
526 return self._eval_row_join(other)
528 def diagonal(self, k=0):
529 """Returns the kth diagonal of self. The main diagonal
530 corresponds to `k=0`; diagonals above and below correspond to
531 `k > 0` and `k < 0`, respectively. The values of `self[i, j]`
532 for which `j - i = k`, are returned in order of increasing
533 `i + j`, starting with `i + j = |k|`.
535 Examples
536 ========
538 >>> from sympy import Matrix
539 >>> m = Matrix(3, 3, lambda i, j: j - i); m
540 Matrix([
541 [ 0, 1, 2],
542 [-1, 0, 1],
543 [-2, -1, 0]])
544 >>> _.diagonal()
545 Matrix([[0, 0, 0]])
546 >>> m.diagonal(1)
547 Matrix([[1, 1]])
548 >>> m.diagonal(-2)
549 Matrix([[-2]])
551 Even though the diagonal is returned as a Matrix, the element
552 retrieval can be done with a single index:
554 >>> Matrix.diag(1, 2, 3).diagonal()[1] # instead of [0, 1]
555 2
557 See Also
558 ========
560 diag
561 """
562 rv = []
563 k = as_int(k)
564 r = 0 if k > 0 else -k
565 c = 0 if r else k
566 while True:
567 if r == self.rows or c == self.cols:
568 break
569 rv.append(self[r, c])
570 r += 1
571 c += 1
572 if not rv:
573 raise ValueError(filldedent('''
574 The %s diagonal is out of range [%s, %s]''' % (
575 k, 1 - self.rows, self.cols - 1)))
576 return self._new(1, len(rv), rv)
578 def row(self, i):
579 """Elementary row selector.
581 Examples
582 ========
584 >>> from sympy import eye
585 >>> eye(2).row(0)
586 Matrix([[1, 0]])
588 See Also
589 ========
591 col
592 row_del
593 row_join
594 row_insert
595 """
596 return self[i, :]
598 @property
599 def shape(self):
600 """The shape (dimensions) of the matrix as the 2-tuple (rows, cols).
602 Examples
603 ========
605 >>> from sympy import zeros
606 >>> M = zeros(2, 3)
607 >>> M.shape
608 (2, 3)
609 >>> M.rows
610 2
611 >>> M.cols
612 3
613 """
614 return (self.rows, self.cols)
616 def todok(self):
617 """Return the matrix as dictionary of keys.
619 Examples
620 ========
622 >>> from sympy import Matrix
623 >>> M = Matrix.eye(3)
624 >>> M.todok()
625 {(0, 0): 1, (1, 1): 1, (2, 2): 1}
626 """
627 return self._eval_todok()
629 def tolist(self):
630 """Return the Matrix as a nested Python list.
632 Examples
633 ========
635 >>> from sympy import Matrix, ones
636 >>> m = Matrix(3, 3, range(9))
637 >>> m
638 Matrix([
639 [0, 1, 2],
640 [3, 4, 5],
641 [6, 7, 8]])
642 >>> m.tolist()
643 [[0, 1, 2], [3, 4, 5], [6, 7, 8]]
644 >>> ones(3, 0).tolist()
645 [[], [], []]
647 When there are no rows then it will not be possible to tell how
648 many columns were in the original matrix:
650 >>> ones(0, 3).tolist()
651 []
653 """
654 if not self.rows:
655 return []
656 if not self.cols:
657 return [[] for i in range(self.rows)]
658 return self._eval_tolist()
660 def todod(M):
661 """Returns matrix as dict of dicts containing non-zero elements of the Matrix
663 Examples
664 ========
666 >>> from sympy import Matrix
667 >>> A = Matrix([[0, 1],[0, 3]])
668 >>> A
669 Matrix([
670 [0, 1],
671 [0, 3]])
672 >>> A.todod()
673 {0: {1: 1}, 1: {1: 3}}
676 """
677 rowsdict = {}
678 Mlol = M.tolist()
679 for i, Mi in enumerate(Mlol):
680 row = {j: Mij for j, Mij in enumerate(Mi) if Mij}
681 if row:
682 rowsdict[i] = row
683 return rowsdict
685 def vec(self):
686 """Return the Matrix converted into a one column matrix by stacking columns
688 Examples
689 ========
691 >>> from sympy import Matrix
692 >>> m=Matrix([[1, 3], [2, 4]])
693 >>> m
694 Matrix([
695 [1, 3],
696 [2, 4]])
697 >>> m.vec()
698 Matrix([
699 [1],
700 [2],
701 [3],
702 [4]])
704 See Also
705 ========
707 vech
708 """
709 return self._eval_vec()
711 def vech(self, diagonal=True, check_symmetry=True):
712 """Reshapes the matrix into a column vector by stacking the
713 elements in the lower triangle.
715 Parameters
716 ==========
718 diagonal : bool, optional
719 If ``True``, it includes the diagonal elements.
721 check_symmetry : bool, optional
722 If ``True``, it checks whether the matrix is symmetric.
724 Examples
725 ========
727 >>> from sympy import Matrix
728 >>> m=Matrix([[1, 2], [2, 3]])
729 >>> m
730 Matrix([
731 [1, 2],
732 [2, 3]])
733 >>> m.vech()
734 Matrix([
735 [1],
736 [2],
737 [3]])
738 >>> m.vech(diagonal=False)
739 Matrix([[2]])
741 Notes
742 =====
744 This should work for symmetric matrices and ``vech`` can
745 represent symmetric matrices in vector form with less size than
746 ``vec``.
748 See Also
749 ========
751 vec
752 """
753 if not self.is_square:
754 raise NonSquareMatrixError
756 if check_symmetry and not self.is_symmetric():
757 raise ValueError("The matrix is not symmetric.")
759 return self._eval_vech(diagonal)
761 @classmethod
762 def vstack(cls, *args):
763 """Return a matrix formed by joining args vertically (i.e.
764 by repeated application of col_join).
766 Examples
767 ========
769 >>> from sympy import Matrix, eye
770 >>> Matrix.vstack(eye(2), 2*eye(2))
771 Matrix([
772 [1, 0],
773 [0, 1],
774 [2, 0],
775 [0, 2]])
776 """
777 if len(args) == 0:
778 return cls._new()
780 kls = type(args[0])
781 return reduce(kls.col_join, args)
784class MatrixSpecial(MatrixRequired):
785 """Construction of special matrices"""
787 @classmethod
788 def _eval_diag(cls, rows, cols, diag_dict):
789 """diag_dict is a defaultdict containing
790 all the entries of the diagonal matrix."""
791 def entry(i, j):
792 return diag_dict[(i, j)]
793 return cls._new(rows, cols, entry)
795 @classmethod
796 def _eval_eye(cls, rows, cols):
797 vals = [cls.zero]*(rows*cols)
798 vals[::cols+1] = [cls.one]*min(rows, cols)
799 return cls._new(rows, cols, vals, copy=False)
801 @classmethod
802 def _eval_jordan_block(cls, size: int, eigenvalue, band='upper'):
803 if band == 'lower':
804 def entry(i, j):
805 if i == j:
806 return eigenvalue
807 elif j + 1 == i:
808 return cls.one
809 return cls.zero
810 else:
811 def entry(i, j):
812 if i == j:
813 return eigenvalue
814 elif i + 1 == j:
815 return cls.one
816 return cls.zero
817 return cls._new(size, size, entry)
819 @classmethod
820 def _eval_ones(cls, rows, cols):
821 def entry(i, j):
822 return cls.one
823 return cls._new(rows, cols, entry)
825 @classmethod
826 def _eval_zeros(cls, rows, cols):
827 return cls._new(rows, cols, [cls.zero]*(rows*cols), copy=False)
829 @classmethod
830 def _eval_wilkinson(cls, n):
831 def entry(i, j):
832 return cls.one if i + 1 == j else cls.zero
834 D = cls._new(2*n + 1, 2*n + 1, entry)
836 wminus = cls.diag(list(range(-n, n + 1)), unpack=True) + D + D.T
837 wplus = abs(cls.diag(list(range(-n, n + 1)), unpack=True)) + D + D.T
839 return wminus, wplus
841 @classmethod
842 def diag(kls, *args, strict=False, unpack=True, rows=None, cols=None, **kwargs):
843 """Returns a matrix with the specified diagonal.
844 If matrices are passed, a block-diagonal matrix
845 is created (i.e. the "direct sum" of the matrices).
847 kwargs
848 ======
850 rows : rows of the resulting matrix; computed if
851 not given.
853 cols : columns of the resulting matrix; computed if
854 not given.
856 cls : class for the resulting matrix
858 unpack : bool which, when True (default), unpacks a single
859 sequence rather than interpreting it as a Matrix.
861 strict : bool which, when False (default), allows Matrices to
862 have variable-length rows.
864 Examples
865 ========
867 >>> from sympy import Matrix
868 >>> Matrix.diag(1, 2, 3)
869 Matrix([
870 [1, 0, 0],
871 [0, 2, 0],
872 [0, 0, 3]])
874 The current default is to unpack a single sequence. If this is
875 not desired, set `unpack=False` and it will be interpreted as
876 a matrix.
878 >>> Matrix.diag([1, 2, 3]) == Matrix.diag(1, 2, 3)
879 True
881 When more than one element is passed, each is interpreted as
882 something to put on the diagonal. Lists are converted to
883 matrices. Filling of the diagonal always continues from
884 the bottom right hand corner of the previous item: this
885 will create a block-diagonal matrix whether the matrices
886 are square or not.
888 >>> col = [1, 2, 3]
889 >>> row = [[4, 5]]
890 >>> Matrix.diag(col, row)
891 Matrix([
892 [1, 0, 0],
893 [2, 0, 0],
894 [3, 0, 0],
895 [0, 4, 5]])
897 When `unpack` is False, elements within a list need not all be
898 of the same length. Setting `strict` to True would raise a
899 ValueError for the following:
901 >>> Matrix.diag([[1, 2, 3], [4, 5], [6]], unpack=False)
902 Matrix([
903 [1, 2, 3],
904 [4, 5, 0],
905 [6, 0, 0]])
907 The type of the returned matrix can be set with the ``cls``
908 keyword.
910 >>> from sympy import ImmutableMatrix
911 >>> from sympy.utilities.misc import func_name
912 >>> func_name(Matrix.diag(1, cls=ImmutableMatrix))
913 'ImmutableDenseMatrix'
915 A zero dimension matrix can be used to position the start of
916 the filling at the start of an arbitrary row or column:
918 >>> from sympy import ones
919 >>> r2 = ones(0, 2)
920 >>> Matrix.diag(r2, 1, 2)
921 Matrix([
922 [0, 0, 1, 0],
923 [0, 0, 0, 2]])
925 See Also
926 ========
927 eye
928 diagonal
929 .dense.diag
930 .expressions.blockmatrix.BlockMatrix
931 .sparsetools.banded
932 """
933 from sympy.matrices.matrices import MatrixBase
934 from sympy.matrices.dense import Matrix
935 from sympy.matrices import SparseMatrix
936 klass = kwargs.get('cls', kls)
937 if unpack and len(args) == 1 and is_sequence(args[0]) and \
938 not isinstance(args[0], MatrixBase):
939 args = args[0]
941 # fill a default dict with the diagonal entries
942 diag_entries = defaultdict(int)
943 rmax = cmax = 0 # keep track of the biggest index seen
944 for m in args:
945 if isinstance(m, list):
946 if strict:
947 # if malformed, Matrix will raise an error
948 _ = Matrix(m)
949 r, c = _.shape
950 m = _.tolist()
951 else:
952 r, c, smat = SparseMatrix._handle_creation_inputs(m)
953 for (i, j), _ in smat.items():
954 diag_entries[(i + rmax, j + cmax)] = _
955 m = [] # to skip process below
956 elif hasattr(m, 'shape'): # a Matrix
957 # convert to list of lists
958 r, c = m.shape
959 m = m.tolist()
960 else: # in this case, we're a single value
961 diag_entries[(rmax, cmax)] = m
962 rmax += 1
963 cmax += 1
964 continue
965 # process list of lists
966 for i, mi in enumerate(m):
967 for j, _ in enumerate(mi):
968 diag_entries[(i + rmax, j + cmax)] = _
969 rmax += r
970 cmax += c
971 if rows is None:
972 rows, cols = cols, rows
973 if rows is None:
974 rows, cols = rmax, cmax
975 else:
976 cols = rows if cols is None else cols
977 if rows < rmax or cols < cmax:
978 raise ValueError(filldedent('''
979 The constructed matrix is {} x {} but a size of {} x {}
980 was specified.'''.format(rmax, cmax, rows, cols)))
981 return klass._eval_diag(rows, cols, diag_entries)
983 @classmethod
984 def eye(kls, rows, cols=None, **kwargs):
985 """Returns an identity matrix.
987 Parameters
988 ==========
990 rows : rows of the matrix
991 cols : cols of the matrix (if None, cols=rows)
993 kwargs
994 ======
995 cls : class of the returned matrix
996 """
997 if cols is None:
998 cols = rows
999 if rows < 0 or cols < 0:
1000 raise ValueError("Cannot create a {} x {} matrix. "
1001 "Both dimensions must be positive".format(rows, cols))
1002 klass = kwargs.get('cls', kls)
1003 rows, cols = as_int(rows), as_int(cols)
1005 return klass._eval_eye(rows, cols)
1007 @classmethod
1008 def jordan_block(kls, size=None, eigenvalue=None, *, band='upper', **kwargs):
1009 """Returns a Jordan block
1011 Parameters
1012 ==========
1014 size : Integer, optional
1015 Specifies the shape of the Jordan block matrix.
1017 eigenvalue : Number or Symbol
1018 Specifies the value for the main diagonal of the matrix.
1020 .. note::
1021 The keyword ``eigenval`` is also specified as an alias
1022 of this keyword, but it is not recommended to use.
1024 We may deprecate the alias in later release.
1026 band : 'upper' or 'lower', optional
1027 Specifies the position of the off-diagonal to put `1` s on.
1029 cls : Matrix, optional
1030 Specifies the matrix class of the output form.
1032 If it is not specified, the class type where the method is
1033 being executed on will be returned.
1035 Returns
1036 =======
1038 Matrix
1039 A Jordan block matrix.
1041 Raises
1042 ======
1044 ValueError
1045 If insufficient arguments are given for matrix size
1046 specification, or no eigenvalue is given.
1048 Examples
1049 ========
1051 Creating a default Jordan block:
1053 >>> from sympy import Matrix
1054 >>> from sympy.abc import x
1055 >>> Matrix.jordan_block(4, x)
1056 Matrix([
1057 [x, 1, 0, 0],
1058 [0, x, 1, 0],
1059 [0, 0, x, 1],
1060 [0, 0, 0, x]])
1062 Creating an alternative Jordan block matrix where `1` is on
1063 lower off-diagonal:
1065 >>> Matrix.jordan_block(4, x, band='lower')
1066 Matrix([
1067 [x, 0, 0, 0],
1068 [1, x, 0, 0],
1069 [0, 1, x, 0],
1070 [0, 0, 1, x]])
1072 Creating a Jordan block with keyword arguments
1074 >>> Matrix.jordan_block(size=4, eigenvalue=x)
1075 Matrix([
1076 [x, 1, 0, 0],
1077 [0, x, 1, 0],
1078 [0, 0, x, 1],
1079 [0, 0, 0, x]])
1081 References
1082 ==========
1084 .. [1] https://en.wikipedia.org/wiki/Jordan_matrix
1085 """
1086 klass = kwargs.pop('cls', kls)
1088 eigenval = kwargs.get('eigenval', None)
1089 if eigenvalue is None and eigenval is None:
1090 raise ValueError("Must supply an eigenvalue")
1091 elif eigenvalue != eigenval and None not in (eigenval, eigenvalue):
1092 raise ValueError(
1093 "Inconsistent values are given: 'eigenval'={}, "
1094 "'eigenvalue'={}".format(eigenval, eigenvalue))
1095 else:
1096 if eigenval is not None:
1097 eigenvalue = eigenval
1099 if size is None:
1100 raise ValueError("Must supply a matrix size")
1102 size = as_int(size)
1103 return klass._eval_jordan_block(size, eigenvalue, band)
1105 @classmethod
1106 def ones(kls, rows, cols=None, **kwargs):
1107 """Returns a matrix of ones.
1109 Parameters
1110 ==========
1112 rows : rows of the matrix
1113 cols : cols of the matrix (if None, cols=rows)
1115 kwargs
1116 ======
1117 cls : class of the returned matrix
1118 """
1119 if cols is None:
1120 cols = rows
1121 klass = kwargs.get('cls', kls)
1122 rows, cols = as_int(rows), as_int(cols)
1124 return klass._eval_ones(rows, cols)
1126 @classmethod
1127 def zeros(kls, rows, cols=None, **kwargs):
1128 """Returns a matrix of zeros.
1130 Parameters
1131 ==========
1133 rows : rows of the matrix
1134 cols : cols of the matrix (if None, cols=rows)
1136 kwargs
1137 ======
1138 cls : class of the returned matrix
1139 """
1140 if cols is None:
1141 cols = rows
1142 if rows < 0 or cols < 0:
1143 raise ValueError("Cannot create a {} x {} matrix. "
1144 "Both dimensions must be positive".format(rows, cols))
1145 klass = kwargs.get('cls', kls)
1146 rows, cols = as_int(rows), as_int(cols)
1148 return klass._eval_zeros(rows, cols)
1150 @classmethod
1151 def companion(kls, poly):
1152 """Returns a companion matrix of a polynomial.
1154 Examples
1155 ========
1157 >>> from sympy import Matrix, Poly, Symbol, symbols
1158 >>> x = Symbol('x')
1159 >>> c0, c1, c2, c3, c4 = symbols('c0:5')
1160 >>> p = Poly(c0 + c1*x + c2*x**2 + c3*x**3 + c4*x**4 + x**5, x)
1161 >>> Matrix.companion(p)
1162 Matrix([
1163 [0, 0, 0, 0, -c0],
1164 [1, 0, 0, 0, -c1],
1165 [0, 1, 0, 0, -c2],
1166 [0, 0, 1, 0, -c3],
1167 [0, 0, 0, 1, -c4]])
1168 """
1169 poly = kls._sympify(poly)
1170 if not isinstance(poly, Poly):
1171 raise ValueError("{} must be a Poly instance.".format(poly))
1172 if not poly.is_monic:
1173 raise ValueError("{} must be a monic polynomial.".format(poly))
1174 if not poly.is_univariate:
1175 raise ValueError(
1176 "{} must be a univariate polynomial.".format(poly))
1178 size = poly.degree()
1179 if not size >= 1:
1180 raise ValueError(
1181 "{} must have degree not less than 1.".format(poly))
1183 coeffs = poly.all_coeffs()
1184 def entry(i, j):
1185 if j == size - 1:
1186 return -coeffs[-1 - i]
1187 elif i == j + 1:
1188 return kls.one
1189 return kls.zero
1190 return kls._new(size, size, entry)
1193 @classmethod
1194 def wilkinson(kls, n, **kwargs):
1195 """Returns two square Wilkinson Matrix of size 2*n + 1
1196 $W_{2n + 1}^-, W_{2n + 1}^+ =$ Wilkinson(n)
1198 Examples
1199 ========
1201 >>> from sympy import Matrix
1202 >>> wminus, wplus = Matrix.wilkinson(3)
1203 >>> wminus
1204 Matrix([
1205 [-3, 1, 0, 0, 0, 0, 0],
1206 [ 1, -2, 1, 0, 0, 0, 0],
1207 [ 0, 1, -1, 1, 0, 0, 0],
1208 [ 0, 0, 1, 0, 1, 0, 0],
1209 [ 0, 0, 0, 1, 1, 1, 0],
1210 [ 0, 0, 0, 0, 1, 2, 1],
1211 [ 0, 0, 0, 0, 0, 1, 3]])
1212 >>> wplus
1213 Matrix([
1214 [3, 1, 0, 0, 0, 0, 0],
1215 [1, 2, 1, 0, 0, 0, 0],
1216 [0, 1, 1, 1, 0, 0, 0],
1217 [0, 0, 1, 0, 1, 0, 0],
1218 [0, 0, 0, 1, 1, 1, 0],
1219 [0, 0, 0, 0, 1, 2, 1],
1220 [0, 0, 0, 0, 0, 1, 3]])
1222 References
1223 ==========
1225 .. [1] https://blogs.mathworks.com/cleve/2013/04/15/wilkinsons-matrices-2/
1226 .. [2] J. H. Wilkinson, The Algebraic Eigenvalue Problem, Claredon Press, Oxford, 1965, 662 pp.
1228 """
1229 klass = kwargs.get('cls', kls)
1230 n = as_int(n)
1231 return klass._eval_wilkinson(n)
1233class MatrixProperties(MatrixRequired):
1234 """Provides basic properties of a matrix."""
1236 def _eval_atoms(self, *types):
1237 result = set()
1238 for i in self:
1239 result.update(i.atoms(*types))
1240 return result
1242 def _eval_free_symbols(self):
1243 return set().union(*(i.free_symbols for i in self if i))
1245 def _eval_has(self, *patterns):
1246 return any(a.has(*patterns) for a in self)
1248 def _eval_is_anti_symmetric(self, simpfunc):
1249 if not all(simpfunc(self[i, j] + self[j, i]).is_zero for i in range(self.rows) for j in range(self.cols)):
1250 return False
1251 return True
1253 def _eval_is_diagonal(self):
1254 for i in range(self.rows):
1255 for j in range(self.cols):
1256 if i != j and self[i, j]:
1257 return False
1258 return True
1260 # _eval_is_hermitian is called by some general SymPy
1261 # routines and has a different *args signature. Make
1262 # sure the names don't clash by adding `_matrix_` in name.
1263 def _eval_is_matrix_hermitian(self, simpfunc):
1264 mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i].conjugate()))
1265 return mat.is_zero_matrix
1267 def _eval_is_Identity(self) -> FuzzyBool:
1268 def dirac(i, j):
1269 if i == j:
1270 return 1
1271 return 0
1273 return all(self[i, j] == dirac(i, j)
1274 for i in range(self.rows)
1275 for j in range(self.cols))
1277 def _eval_is_lower_hessenberg(self):
1278 return all(self[i, j].is_zero
1279 for i in range(self.rows)
1280 for j in range(i + 2, self.cols))
1282 def _eval_is_lower(self):
1283 return all(self[i, j].is_zero
1284 for i in range(self.rows)
1285 for j in range(i + 1, self.cols))
1287 def _eval_is_symbolic(self):
1288 return self.has(Symbol)
1290 def _eval_is_symmetric(self, simpfunc):
1291 mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i]))
1292 return mat.is_zero_matrix
1294 def _eval_is_zero_matrix(self):
1295 if any(i.is_zero == False for i in self):
1296 return False
1297 if any(i.is_zero is None for i in self):
1298 return None
1299 return True
1301 def _eval_is_upper_hessenberg(self):
1302 return all(self[i, j].is_zero
1303 for i in range(2, self.rows)
1304 for j in range(min(self.cols, (i - 1))))
1306 def _eval_values(self):
1307 return [i for i in self if not i.is_zero]
1309 def _has_positive_diagonals(self):
1310 diagonal_entries = (self[i, i] for i in range(self.rows))
1311 return fuzzy_and(x.is_positive for x in diagonal_entries)
1313 def _has_nonnegative_diagonals(self):
1314 diagonal_entries = (self[i, i] for i in range(self.rows))
1315 return fuzzy_and(x.is_nonnegative for x in diagonal_entries)
1317 def atoms(self, *types):
1318 """Returns the atoms that form the current object.
1320 Examples
1321 ========
1323 >>> from sympy.abc import x, y
1324 >>> from sympy import Matrix
1325 >>> Matrix([[x]])
1326 Matrix([[x]])
1327 >>> _.atoms()
1328 {x}
1329 >>> Matrix([[x, y], [y, x]])
1330 Matrix([
1331 [x, y],
1332 [y, x]])
1333 >>> _.atoms()
1334 {x, y}
1335 """
1337 types = tuple(t if isinstance(t, type) else type(t) for t in types)
1338 if not types:
1339 types = (Atom,)
1340 return self._eval_atoms(*types)
1342 @property
1343 def free_symbols(self):
1344 """Returns the free symbols within the matrix.
1346 Examples
1347 ========
1349 >>> from sympy.abc import x
1350 >>> from sympy import Matrix
1351 >>> Matrix([[x], [1]]).free_symbols
1352 {x}
1353 """
1354 return self._eval_free_symbols()
1356 def has(self, *patterns):
1357 """Test whether any subexpression matches any of the patterns.
1359 Examples
1360 ========
1362 >>> from sympy import Matrix, SparseMatrix, Float
1363 >>> from sympy.abc import x, y
1364 >>> A = Matrix(((1, x), (0.2, 3)))
1365 >>> B = SparseMatrix(((1, x), (0.2, 3)))
1366 >>> A.has(x)
1367 True
1368 >>> A.has(y)
1369 False
1370 >>> A.has(Float)
1371 True
1372 >>> B.has(x)
1373 True
1374 >>> B.has(y)
1375 False
1376 >>> B.has(Float)
1377 True
1378 """
1379 return self._eval_has(*patterns)
1381 def is_anti_symmetric(self, simplify=True):
1382 """Check if matrix M is an antisymmetric matrix,
1383 that is, M is a square matrix with all M[i, j] == -M[j, i].
1385 When ``simplify=True`` (default), the sum M[i, j] + M[j, i] is
1386 simplified before testing to see if it is zero. By default,
1387 the SymPy simplify function is used. To use a custom function
1388 set simplify to a function that accepts a single argument which
1389 returns a simplified expression. To skip simplification, set
1390 simplify to False but note that although this will be faster,
1391 it may induce false negatives.
1393 Examples
1394 ========
1396 >>> from sympy import Matrix, symbols
1397 >>> m = Matrix(2, 2, [0, 1, -1, 0])
1398 >>> m
1399 Matrix([
1400 [ 0, 1],
1401 [-1, 0]])
1402 >>> m.is_anti_symmetric()
1403 True
1404 >>> x, y = symbols('x y')
1405 >>> m = Matrix(2, 3, [0, 0, x, -y, 0, 0])
1406 >>> m
1407 Matrix([
1408 [ 0, 0, x],
1409 [-y, 0, 0]])
1410 >>> m.is_anti_symmetric()
1411 False
1413 >>> from sympy.abc import x, y
1414 >>> m = Matrix(3, 3, [0, x**2 + 2*x + 1, y,
1415 ... -(x + 1)**2, 0, x*y,
1416 ... -y, -x*y, 0])
1418 Simplification of matrix elements is done by default so even
1419 though two elements which should be equal and opposite would not
1420 pass an equality test, the matrix is still reported as
1421 anti-symmetric:
1423 >>> m[0, 1] == -m[1, 0]
1424 False
1425 >>> m.is_anti_symmetric()
1426 True
1428 If ``simplify=False`` is used for the case when a Matrix is already
1429 simplified, this will speed things up. Here, we see that without
1430 simplification the matrix does not appear anti-symmetric:
1432 >>> m.is_anti_symmetric(simplify=False)
1433 False
1435 But if the matrix were already expanded, then it would appear
1436 anti-symmetric and simplification in the is_anti_symmetric routine
1437 is not needed:
1439 >>> m = m.expand()
1440 >>> m.is_anti_symmetric(simplify=False)
1441 True
1442 """
1443 # accept custom simplification
1444 simpfunc = simplify
1445 if not isfunction(simplify):
1446 simpfunc = _simplify if simplify else lambda x: x
1448 if not self.is_square:
1449 return False
1450 return self._eval_is_anti_symmetric(simpfunc)
1452 def is_diagonal(self):
1453 """Check if matrix is diagonal,
1454 that is matrix in which the entries outside the main diagonal are all zero.
1456 Examples
1457 ========
1459 >>> from sympy import Matrix, diag
1460 >>> m = Matrix(2, 2, [1, 0, 0, 2])
1461 >>> m
1462 Matrix([
1463 [1, 0],
1464 [0, 2]])
1465 >>> m.is_diagonal()
1466 True
1468 >>> m = Matrix(2, 2, [1, 1, 0, 2])
1469 >>> m
1470 Matrix([
1471 [1, 1],
1472 [0, 2]])
1473 >>> m.is_diagonal()
1474 False
1476 >>> m = diag(1, 2, 3)
1477 >>> m
1478 Matrix([
1479 [1, 0, 0],
1480 [0, 2, 0],
1481 [0, 0, 3]])
1482 >>> m.is_diagonal()
1483 True
1485 See Also
1486 ========
1488 is_lower
1489 is_upper
1490 sympy.matrices.matrices.MatrixEigen.is_diagonalizable
1491 diagonalize
1492 """
1493 return self._eval_is_diagonal()
1495 @property
1496 def is_weakly_diagonally_dominant(self):
1497 r"""Tests if the matrix is row weakly diagonally dominant.
1499 Explanation
1500 ===========
1502 A $n, n$ matrix $A$ is row weakly diagonally dominant if
1504 .. math::
1505 \left|A_{i, i}\right| \ge \sum_{j = 0, j \neq i}^{n-1}
1506 \left|A_{i, j}\right| \quad {\text{for all }}
1507 i \in \{ 0, ..., n-1 \}
1509 Examples
1510 ========
1512 >>> from sympy import Matrix
1513 >>> A = Matrix([[3, -2, 1], [1, -3, 2], [-1, 2, 4]])
1514 >>> A.is_weakly_diagonally_dominant
1515 True
1517 >>> A = Matrix([[-2, 2, 1], [1, 3, 2], [1, -2, 0]])
1518 >>> A.is_weakly_diagonally_dominant
1519 False
1521 >>> A = Matrix([[-4, 2, 1], [1, 6, 2], [1, -2, 5]])
1522 >>> A.is_weakly_diagonally_dominant
1523 True
1525 Notes
1526 =====
1528 If you want to test whether a matrix is column diagonally
1529 dominant, you can apply the test after transposing the matrix.
1530 """
1531 if not self.is_square:
1532 return False
1534 rows, cols = self.shape
1536 def test_row(i):
1537 summation = self.zero
1538 for j in range(cols):
1539 if i != j:
1540 summation += Abs(self[i, j])
1541 return (Abs(self[i, i]) - summation).is_nonnegative
1543 return fuzzy_and(test_row(i) for i in range(rows))
1545 @property
1546 def is_strongly_diagonally_dominant(self):
1547 r"""Tests if the matrix is row strongly diagonally dominant.
1549 Explanation
1550 ===========
1552 A $n, n$ matrix $A$ is row strongly diagonally dominant if
1554 .. math::
1555 \left|A_{i, i}\right| > \sum_{j = 0, j \neq i}^{n-1}
1556 \left|A_{i, j}\right| \quad {\text{for all }}
1557 i \in \{ 0, ..., n-1 \}
1559 Examples
1560 ========
1562 >>> from sympy import Matrix
1563 >>> A = Matrix([[3, -2, 1], [1, -3, 2], [-1, 2, 4]])
1564 >>> A.is_strongly_diagonally_dominant
1565 False
1567 >>> A = Matrix([[-2, 2, 1], [1, 3, 2], [1, -2, 0]])
1568 >>> A.is_strongly_diagonally_dominant
1569 False
1571 >>> A = Matrix([[-4, 2, 1], [1, 6, 2], [1, -2, 5]])
1572 >>> A.is_strongly_diagonally_dominant
1573 True
1575 Notes
1576 =====
1578 If you want to test whether a matrix is column diagonally
1579 dominant, you can apply the test after transposing the matrix.
1580 """
1581 if not self.is_square:
1582 return False
1584 rows, cols = self.shape
1586 def test_row(i):
1587 summation = self.zero
1588 for j in range(cols):
1589 if i != j:
1590 summation += Abs(self[i, j])
1591 return (Abs(self[i, i]) - summation).is_positive
1593 return fuzzy_and(test_row(i) for i in range(rows))
1595 @property
1596 def is_hermitian(self):
1597 """Checks if the matrix is Hermitian.
1599 In a Hermitian matrix element i,j is the complex conjugate of
1600 element j,i.
1602 Examples
1603 ========
1605 >>> from sympy import Matrix
1606 >>> from sympy import I
1607 >>> from sympy.abc import x
1608 >>> a = Matrix([[1, I], [-I, 1]])
1609 >>> a
1610 Matrix([
1611 [ 1, I],
1612 [-I, 1]])
1613 >>> a.is_hermitian
1614 True
1615 >>> a[0, 0] = 2*I
1616 >>> a.is_hermitian
1617 False
1618 >>> a[0, 0] = x
1619 >>> a.is_hermitian
1620 >>> a[0, 1] = a[1, 0]*I
1621 >>> a.is_hermitian
1622 False
1623 """
1624 if not self.is_square:
1625 return False
1627 return self._eval_is_matrix_hermitian(_simplify)
1629 @property
1630 def is_Identity(self) -> FuzzyBool:
1631 if not self.is_square:
1632 return False
1633 return self._eval_is_Identity()
1635 @property
1636 def is_lower_hessenberg(self):
1637 r"""Checks if the matrix is in the lower-Hessenberg form.
1639 The lower hessenberg matrix has zero entries
1640 above the first superdiagonal.
1642 Examples
1643 ========
1645 >>> from sympy import Matrix
1646 >>> a = Matrix([[1, 2, 0, 0], [5, 2, 3, 0], [3, 4, 3, 7], [5, 6, 1, 1]])
1647 >>> a
1648 Matrix([
1649 [1, 2, 0, 0],
1650 [5, 2, 3, 0],
1651 [3, 4, 3, 7],
1652 [5, 6, 1, 1]])
1653 >>> a.is_lower_hessenberg
1654 True
1656 See Also
1657 ========
1659 is_upper_hessenberg
1660 is_lower
1661 """
1662 return self._eval_is_lower_hessenberg()
1664 @property
1665 def is_lower(self):
1666 """Check if matrix is a lower triangular matrix. True can be returned
1667 even if the matrix is not square.
1669 Examples
1670 ========
1672 >>> from sympy import Matrix
1673 >>> m = Matrix(2, 2, [1, 0, 0, 1])
1674 >>> m
1675 Matrix([
1676 [1, 0],
1677 [0, 1]])
1678 >>> m.is_lower
1679 True
1681 >>> m = Matrix(4, 3, [0, 0, 0, 2, 0, 0, 1, 4, 0, 6, 6, 5])
1682 >>> m
1683 Matrix([
1684 [0, 0, 0],
1685 [2, 0, 0],
1686 [1, 4, 0],
1687 [6, 6, 5]])
1688 >>> m.is_lower
1689 True
1691 >>> from sympy.abc import x, y
1692 >>> m = Matrix(2, 2, [x**2 + y, y**2 + x, 0, x + y])
1693 >>> m
1694 Matrix([
1695 [x**2 + y, x + y**2],
1696 [ 0, x + y]])
1697 >>> m.is_lower
1698 False
1700 See Also
1701 ========
1703 is_upper
1704 is_diagonal
1705 is_lower_hessenberg
1706 """
1707 return self._eval_is_lower()
1709 @property
1710 def is_square(self):
1711 """Checks if a matrix is square.
1713 A matrix is square if the number of rows equals the number of columns.
1714 The empty matrix is square by definition, since the number of rows and
1715 the number of columns are both zero.
1717 Examples
1718 ========
1720 >>> from sympy import Matrix
1721 >>> a = Matrix([[1, 2, 3], [4, 5, 6]])
1722 >>> b = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
1723 >>> c = Matrix([])
1724 >>> a.is_square
1725 False
1726 >>> b.is_square
1727 True
1728 >>> c.is_square
1729 True
1730 """
1731 return self.rows == self.cols
1733 def is_symbolic(self):
1734 """Checks if any elements contain Symbols.
1736 Examples
1737 ========
1739 >>> from sympy import Matrix
1740 >>> from sympy.abc import x, y
1741 >>> M = Matrix([[x, y], [1, 0]])
1742 >>> M.is_symbolic()
1743 True
1745 """
1746 return self._eval_is_symbolic()
1748 def is_symmetric(self, simplify=True):
1749 """Check if matrix is symmetric matrix,
1750 that is square matrix and is equal to its transpose.
1752 By default, simplifications occur before testing symmetry.
1753 They can be skipped using 'simplify=False'; while speeding things a bit,
1754 this may however induce false negatives.
1756 Examples
1757 ========
1759 >>> from sympy import Matrix
1760 >>> m = Matrix(2, 2, [0, 1, 1, 2])
1761 >>> m
1762 Matrix([
1763 [0, 1],
1764 [1, 2]])
1765 >>> m.is_symmetric()
1766 True
1768 >>> m = Matrix(2, 2, [0, 1, 2, 0])
1769 >>> m
1770 Matrix([
1771 [0, 1],
1772 [2, 0]])
1773 >>> m.is_symmetric()
1774 False
1776 >>> m = Matrix(2, 3, [0, 0, 0, 0, 0, 0])
1777 >>> m
1778 Matrix([
1779 [0, 0, 0],
1780 [0, 0, 0]])
1781 >>> m.is_symmetric()
1782 False
1784 >>> from sympy.abc import x, y
1785 >>> m = Matrix(3, 3, [1, x**2 + 2*x + 1, y, (x + 1)**2, 2, 0, y, 0, 3])
1786 >>> m
1787 Matrix([
1788 [ 1, x**2 + 2*x + 1, y],
1789 [(x + 1)**2, 2, 0],
1790 [ y, 0, 3]])
1791 >>> m.is_symmetric()
1792 True
1794 If the matrix is already simplified, you may speed-up is_symmetric()
1795 test by using 'simplify=False'.
1797 >>> bool(m.is_symmetric(simplify=False))
1798 False
1799 >>> m1 = m.expand()
1800 >>> m1.is_symmetric(simplify=False)
1801 True
1802 """
1803 simpfunc = simplify
1804 if not isfunction(simplify):
1805 simpfunc = _simplify if simplify else lambda x: x
1807 if not self.is_square:
1808 return False
1810 return self._eval_is_symmetric(simpfunc)
1812 @property
1813 def is_upper_hessenberg(self):
1814 """Checks if the matrix is the upper-Hessenberg form.
1816 The upper hessenberg matrix has zero entries
1817 below the first subdiagonal.
1819 Examples
1820 ========
1822 >>> from sympy import Matrix
1823 >>> a = Matrix([[1, 4, 2, 3], [3, 4, 1, 7], [0, 2, 3, 4], [0, 0, 1, 3]])
1824 >>> a
1825 Matrix([
1826 [1, 4, 2, 3],
1827 [3, 4, 1, 7],
1828 [0, 2, 3, 4],
1829 [0, 0, 1, 3]])
1830 >>> a.is_upper_hessenberg
1831 True
1833 See Also
1834 ========
1836 is_lower_hessenberg
1837 is_upper
1838 """
1839 return self._eval_is_upper_hessenberg()
1841 @property
1842 def is_upper(self):
1843 """Check if matrix is an upper triangular matrix. True can be returned
1844 even if the matrix is not square.
1846 Examples
1847 ========
1849 >>> from sympy import Matrix
1850 >>> m = Matrix(2, 2, [1, 0, 0, 1])
1851 >>> m
1852 Matrix([
1853 [1, 0],
1854 [0, 1]])
1855 >>> m.is_upper
1856 True
1858 >>> m = Matrix(4, 3, [5, 1, 9, 0, 4, 6, 0, 0, 5, 0, 0, 0])
1859 >>> m
1860 Matrix([
1861 [5, 1, 9],
1862 [0, 4, 6],
1863 [0, 0, 5],
1864 [0, 0, 0]])
1865 >>> m.is_upper
1866 True
1868 >>> m = Matrix(2, 3, [4, 2, 5, 6, 1, 1])
1869 >>> m
1870 Matrix([
1871 [4, 2, 5],
1872 [6, 1, 1]])
1873 >>> m.is_upper
1874 False
1876 See Also
1877 ========
1879 is_lower
1880 is_diagonal
1881 is_upper_hessenberg
1882 """
1883 return all(self[i, j].is_zero
1884 for i in range(1, self.rows)
1885 for j in range(min(i, self.cols)))
1887 @property
1888 def is_zero_matrix(self):
1889 """Checks if a matrix is a zero matrix.
1891 A matrix is zero if every element is zero. A matrix need not be square
1892 to be considered zero. The empty matrix is zero by the principle of
1893 vacuous truth. For a matrix that may or may not be zero (e.g.
1894 contains a symbol), this will be None
1896 Examples
1897 ========
1899 >>> from sympy import Matrix, zeros
1900 >>> from sympy.abc import x
1901 >>> a = Matrix([[0, 0], [0, 0]])
1902 >>> b = zeros(3, 4)
1903 >>> c = Matrix([[0, 1], [0, 0]])
1904 >>> d = Matrix([])
1905 >>> e = Matrix([[x, 0], [0, 0]])
1906 >>> a.is_zero_matrix
1907 True
1908 >>> b.is_zero_matrix
1909 True
1910 >>> c.is_zero_matrix
1911 False
1912 >>> d.is_zero_matrix
1913 True
1914 >>> e.is_zero_matrix
1915 """
1916 return self._eval_is_zero_matrix()
1918 def values(self):
1919 """Return non-zero values of self."""
1920 return self._eval_values()
1923class MatrixOperations(MatrixRequired):
1924 """Provides basic matrix shape and elementwise
1925 operations. Should not be instantiated directly."""
1927 def _eval_adjoint(self):
1928 return self.transpose().conjugate()
1930 def _eval_applyfunc(self, f):
1931 out = self._new(self.rows, self.cols, [f(x) for x in self])
1932 return out
1934 def _eval_as_real_imag(self): # type: ignore
1935 return (self.applyfunc(re), self.applyfunc(im))
1937 def _eval_conjugate(self):
1938 return self.applyfunc(lambda x: x.conjugate())
1940 def _eval_permute_cols(self, perm):
1941 # apply the permutation to a list
1942 mapping = list(perm)
1944 def entry(i, j):
1945 return self[i, mapping[j]]
1947 return self._new(self.rows, self.cols, entry)
1949 def _eval_permute_rows(self, perm):
1950 # apply the permutation to a list
1951 mapping = list(perm)
1953 def entry(i, j):
1954 return self[mapping[i], j]
1956 return self._new(self.rows, self.cols, entry)
1958 def _eval_trace(self):
1959 return sum(self[i, i] for i in range(self.rows))
1961 def _eval_transpose(self):
1962 return self._new(self.cols, self.rows, lambda i, j: self[j, i])
1964 def adjoint(self):
1965 """Conjugate transpose or Hermitian conjugation."""
1966 return self._eval_adjoint()
1968 def applyfunc(self, f):
1969 """Apply a function to each element of the matrix.
1971 Examples
1972 ========
1974 >>> from sympy import Matrix
1975 >>> m = Matrix(2, 2, lambda i, j: i*2+j)
1976 >>> m
1977 Matrix([
1978 [0, 1],
1979 [2, 3]])
1980 >>> m.applyfunc(lambda i: 2*i)
1981 Matrix([
1982 [0, 2],
1983 [4, 6]])
1985 """
1986 if not callable(f):
1987 raise TypeError("`f` must be callable.")
1989 return self._eval_applyfunc(f)
1991 def as_real_imag(self, deep=True, **hints):
1992 """Returns a tuple containing the (real, imaginary) part of matrix."""
1993 # XXX: Ignoring deep and hints...
1994 return self._eval_as_real_imag()
1996 def conjugate(self):
1997 """Return the by-element conjugation.
1999 Examples
2000 ========
2002 >>> from sympy import SparseMatrix, I
2003 >>> a = SparseMatrix(((1, 2 + I), (3, 4), (I, -I)))
2004 >>> a
2005 Matrix([
2006 [1, 2 + I],
2007 [3, 4],
2008 [I, -I]])
2009 >>> a.C
2010 Matrix([
2011 [ 1, 2 - I],
2012 [ 3, 4],
2013 [-I, I]])
2015 See Also
2016 ========
2018 transpose: Matrix transposition
2019 H: Hermite conjugation
2020 sympy.matrices.matrices.MatrixBase.D: Dirac conjugation
2021 """
2022 return self._eval_conjugate()
2024 def doit(self, **hints):
2025 return self.applyfunc(lambda x: x.doit(**hints))
2027 def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None, verbose=False):
2028 """Apply evalf() to each element of self."""
2029 options = {'subs':subs, 'maxn':maxn, 'chop':chop, 'strict':strict,
2030 'quad':quad, 'verbose':verbose}
2031 return self.applyfunc(lambda i: i.evalf(n, **options))
2033 def expand(self, deep=True, modulus=None, power_base=True, power_exp=True,
2034 mul=True, log=True, multinomial=True, basic=True, **hints):
2035 """Apply core.function.expand to each entry of the matrix.
2037 Examples
2038 ========
2040 >>> from sympy.abc import x
2041 >>> from sympy import Matrix
2042 >>> Matrix(1, 1, [x*(x+1)])
2043 Matrix([[x*(x + 1)]])
2044 >>> _.expand()
2045 Matrix([[x**2 + x]])
2047 """
2048 return self.applyfunc(lambda x: x.expand(
2049 deep, modulus, power_base, power_exp, mul, log, multinomial, basic,
2050 **hints))
2052 @property
2053 def H(self):
2054 """Return Hermite conjugate.
2056 Examples
2057 ========
2059 >>> from sympy import Matrix, I
2060 >>> m = Matrix((0, 1 + I, 2, 3))
2061 >>> m
2062 Matrix([
2063 [ 0],
2064 [1 + I],
2065 [ 2],
2066 [ 3]])
2067 >>> m.H
2068 Matrix([[0, 1 - I, 2, 3]])
2070 See Also
2071 ========
2073 conjugate: By-element conjugation
2074 sympy.matrices.matrices.MatrixBase.D: Dirac conjugation
2075 """
2076 return self.T.C
2078 def permute(self, perm, orientation='rows', direction='forward'):
2079 r"""Permute the rows or columns of a matrix by the given list of
2080 swaps.
2082 Parameters
2083 ==========
2085 perm : Permutation, list, or list of lists
2086 A representation for the permutation.
2088 If it is ``Permutation``, it is used directly with some
2089 resizing with respect to the matrix size.
2091 If it is specified as list of lists,
2092 (e.g., ``[[0, 1], [0, 2]]``), then the permutation is formed
2093 from applying the product of cycles. The direction how the
2094 cyclic product is applied is described in below.
2096 If it is specified as a list, the list should represent
2097 an array form of a permutation. (e.g., ``[1, 2, 0]``) which
2098 would would form the swapping function
2099 `0 \mapsto 1, 1 \mapsto 2, 2\mapsto 0`.
2101 orientation : 'rows', 'cols'
2102 A flag to control whether to permute the rows or the columns
2104 direction : 'forward', 'backward'
2105 A flag to control whether to apply the permutations from
2106 the start of the list first, or from the back of the list
2107 first.
2109 For example, if the permutation specification is
2110 ``[[0, 1], [0, 2]]``,
2112 If the flag is set to ``'forward'``, the cycle would be
2113 formed as `0 \mapsto 2, 2 \mapsto 1, 1 \mapsto 0`.
2115 If the flag is set to ``'backward'``, the cycle would be
2116 formed as `0 \mapsto 1, 1 \mapsto 2, 2 \mapsto 0`.
2118 If the argument ``perm`` is not in a form of list of lists,
2119 this flag takes no effect.
2121 Examples
2122 ========
2124 >>> from sympy import eye
2125 >>> M = eye(3)
2126 >>> M.permute([[0, 1], [0, 2]], orientation='rows', direction='forward')
2127 Matrix([
2128 [0, 0, 1],
2129 [1, 0, 0],
2130 [0, 1, 0]])
2132 >>> from sympy import eye
2133 >>> M = eye(3)
2134 >>> M.permute([[0, 1], [0, 2]], orientation='rows', direction='backward')
2135 Matrix([
2136 [0, 1, 0],
2137 [0, 0, 1],
2138 [1, 0, 0]])
2140 Notes
2141 =====
2143 If a bijective function
2144 `\sigma : \mathbb{N}_0 \rightarrow \mathbb{N}_0` denotes the
2145 permutation.
2147 If the matrix `A` is the matrix to permute, represented as
2148 a horizontal or a vertical stack of vectors:
2150 .. math::
2151 A =
2152 \begin{bmatrix}
2153 a_0 \\ a_1 \\ \vdots \\ a_{n-1}
2154 \end{bmatrix} =
2155 \begin{bmatrix}
2156 \alpha_0 & \alpha_1 & \cdots & \alpha_{n-1}
2157 \end{bmatrix}
2159 If the matrix `B` is the result, the permutation of matrix rows
2160 is defined as:
2162 .. math::
2163 B := \begin{bmatrix}
2164 a_{\sigma(0)} \\ a_{\sigma(1)} \\ \vdots \\ a_{\sigma(n-1)}
2165 \end{bmatrix}
2167 And the permutation of matrix columns is defined as:
2169 .. math::
2170 B := \begin{bmatrix}
2171 \alpha_{\sigma(0)} & \alpha_{\sigma(1)} &
2172 \cdots & \alpha_{\sigma(n-1)}
2173 \end{bmatrix}
2174 """
2175 from sympy.combinatorics import Permutation
2177 # allow british variants and `columns`
2178 if direction == 'forwards':
2179 direction = 'forward'
2180 if direction == 'backwards':
2181 direction = 'backward'
2182 if orientation == 'columns':
2183 orientation = 'cols'
2185 if direction not in ('forward', 'backward'):
2186 raise TypeError("direction='{}' is an invalid kwarg. "
2187 "Try 'forward' or 'backward'".format(direction))
2188 if orientation not in ('rows', 'cols'):
2189 raise TypeError("orientation='{}' is an invalid kwarg. "
2190 "Try 'rows' or 'cols'".format(orientation))
2192 if not isinstance(perm, (Permutation, Iterable)):
2193 raise ValueError(
2194 "{} must be a list, a list of lists, "
2195 "or a SymPy permutation object.".format(perm))
2197 # ensure all swaps are in range
2198 max_index = self.rows if orientation == 'rows' else self.cols
2199 if not all(0 <= t <= max_index for t in flatten(list(perm))):
2200 raise IndexError("`swap` indices out of range.")
2202 if perm and not isinstance(perm, Permutation) and \
2203 isinstance(perm[0], Iterable):
2204 if direction == 'forward':
2205 perm = list(reversed(perm))
2206 perm = Permutation(perm, size=max_index+1)
2207 else:
2208 perm = Permutation(perm, size=max_index+1)
2210 if orientation == 'rows':
2211 return self._eval_permute_rows(perm)
2212 if orientation == 'cols':
2213 return self._eval_permute_cols(perm)
2215 def permute_cols(self, swaps, direction='forward'):
2216 """Alias for
2217 ``self.permute(swaps, orientation='cols', direction=direction)``
2219 See Also
2220 ========
2222 permute
2223 """
2224 return self.permute(swaps, orientation='cols', direction=direction)
2226 def permute_rows(self, swaps, direction='forward'):
2227 """Alias for
2228 ``self.permute(swaps, orientation='rows', direction=direction)``
2230 See Also
2231 ========
2233 permute
2234 """
2235 return self.permute(swaps, orientation='rows', direction=direction)
2237 def refine(self, assumptions=True):
2238 """Apply refine to each element of the matrix.
2240 Examples
2241 ========
2243 >>> from sympy import Symbol, Matrix, Abs, sqrt, Q
2244 >>> x = Symbol('x')
2245 >>> Matrix([[Abs(x)**2, sqrt(x**2)],[sqrt(x**2), Abs(x)**2]])
2246 Matrix([
2247 [ Abs(x)**2, sqrt(x**2)],
2248 [sqrt(x**2), Abs(x)**2]])
2249 >>> _.refine(Q.real(x))
2250 Matrix([
2251 [ x**2, Abs(x)],
2252 [Abs(x), x**2]])
2254 """
2255 return self.applyfunc(lambda x: refine(x, assumptions))
2257 def replace(self, F, G, map=False, simultaneous=True, exact=None):
2258 """Replaces Function F in Matrix entries with Function G.
2260 Examples
2261 ========
2263 >>> from sympy import symbols, Function, Matrix
2264 >>> F, G = symbols('F, G', cls=Function)
2265 >>> M = Matrix(2, 2, lambda i, j: F(i+j)) ; M
2266 Matrix([
2267 [F(0), F(1)],
2268 [F(1), F(2)]])
2269 >>> N = M.replace(F,G)
2270 >>> N
2271 Matrix([
2272 [G(0), G(1)],
2273 [G(1), G(2)]])
2274 """
2275 return self.applyfunc(
2276 lambda x: x.replace(F, G, map=map, simultaneous=simultaneous, exact=exact))
2278 def rot90(self, k=1):
2279 """Rotates Matrix by 90 degrees
2281 Parameters
2282 ==========
2284 k : int
2285 Specifies how many times the matrix is rotated by 90 degrees
2286 (clockwise when positive, counter-clockwise when negative).
2288 Examples
2289 ========
2291 >>> from sympy import Matrix, symbols
2292 >>> A = Matrix(2, 2, symbols('a:d'))
2293 >>> A
2294 Matrix([
2295 [a, b],
2296 [c, d]])
2298 Rotating the matrix clockwise one time:
2300 >>> A.rot90(1)
2301 Matrix([
2302 [c, a],
2303 [d, b]])
2305 Rotating the matrix anticlockwise two times:
2307 >>> A.rot90(-2)
2308 Matrix([
2309 [d, c],
2310 [b, a]])
2311 """
2313 mod = k%4
2314 if mod == 0:
2315 return self
2316 if mod == 1:
2317 return self[::-1, ::].T
2318 if mod == 2:
2319 return self[::-1, ::-1]
2320 if mod == 3:
2321 return self[::, ::-1].T
2323 def simplify(self, **kwargs):
2324 """Apply simplify to each element of the matrix.
2326 Examples
2327 ========
2329 >>> from sympy.abc import x, y
2330 >>> from sympy import SparseMatrix, sin, cos
2331 >>> SparseMatrix(1, 1, [x*sin(y)**2 + x*cos(y)**2])
2332 Matrix([[x*sin(y)**2 + x*cos(y)**2]])
2333 >>> _.simplify()
2334 Matrix([[x]])
2335 """
2336 return self.applyfunc(lambda x: x.simplify(**kwargs))
2338 def subs(self, *args, **kwargs): # should mirror core.basic.subs
2339 """Return a new matrix with subs applied to each entry.
2341 Examples
2342 ========
2344 >>> from sympy.abc import x, y
2345 >>> from sympy import SparseMatrix, Matrix
2346 >>> SparseMatrix(1, 1, [x])
2347 Matrix([[x]])
2348 >>> _.subs(x, y)
2349 Matrix([[y]])
2350 >>> Matrix(_).subs(y, x)
2351 Matrix([[x]])
2352 """
2354 if len(args) == 1 and not isinstance(args[0], (dict, set)) and iter(args[0]) and not is_sequence(args[0]):
2355 args = (list(args[0]),)
2357 return self.applyfunc(lambda x: x.subs(*args, **kwargs))
2359 def trace(self):
2360 """
2361 Returns the trace of a square matrix i.e. the sum of the
2362 diagonal elements.
2364 Examples
2365 ========
2367 >>> from sympy import Matrix
2368 >>> A = Matrix(2, 2, [1, 2, 3, 4])
2369 >>> A.trace()
2370 5
2372 """
2373 if self.rows != self.cols:
2374 raise NonSquareMatrixError()
2375 return self._eval_trace()
2377 def transpose(self):
2378 """
2379 Returns the transpose of the matrix.
2381 Examples
2382 ========
2384 >>> from sympy import Matrix
2385 >>> A = Matrix(2, 2, [1, 2, 3, 4])
2386 >>> A.transpose()
2387 Matrix([
2388 [1, 3],
2389 [2, 4]])
2391 >>> from sympy import Matrix, I
2392 >>> m=Matrix(((1, 2+I), (3, 4)))
2393 >>> m
2394 Matrix([
2395 [1, 2 + I],
2396 [3, 4]])
2397 >>> m.transpose()
2398 Matrix([
2399 [ 1, 3],
2400 [2 + I, 4]])
2401 >>> m.T == m.transpose()
2402 True
2404 See Also
2405 ========
2407 conjugate: By-element conjugation
2409 """
2410 return self._eval_transpose()
2412 @property
2413 def T(self):
2414 '''Matrix transposition'''
2415 return self.transpose()
2417 @property
2418 def C(self):
2419 '''By-element conjugation'''
2420 return self.conjugate()
2422 def n(self, *args, **kwargs):
2423 """Apply evalf() to each element of self."""
2424 return self.evalf(*args, **kwargs)
2426 def xreplace(self, rule): # should mirror core.basic.xreplace
2427 """Return a new matrix with xreplace applied to each entry.
2429 Examples
2430 ========
2432 >>> from sympy.abc import x, y
2433 >>> from sympy import SparseMatrix, Matrix
2434 >>> SparseMatrix(1, 1, [x])
2435 Matrix([[x]])
2436 >>> _.xreplace({x: y})
2437 Matrix([[y]])
2438 >>> Matrix(_).xreplace({y: x})
2439 Matrix([[x]])
2440 """
2441 return self.applyfunc(lambda x: x.xreplace(rule))
2443 def _eval_simplify(self, **kwargs):
2444 # XXX: We can't use self.simplify here as mutable subclasses will
2445 # override simplify and have it return None
2446 return MatrixOperations.simplify(self, **kwargs)
2448 def _eval_trigsimp(self, **opts):
2449 from sympy.simplify.trigsimp import trigsimp
2450 return self.applyfunc(lambda x: trigsimp(x, **opts))
2452 def upper_triangular(self, k=0):
2453 """Return the elements on and above the kth diagonal of a matrix.
2454 If k is not specified then simply returns upper-triangular portion
2455 of a matrix
2457 Examples
2458 ========
2460 >>> from sympy import ones
2461 >>> A = ones(4)
2462 >>> A.upper_triangular()
2463 Matrix([
2464 [1, 1, 1, 1],
2465 [0, 1, 1, 1],
2466 [0, 0, 1, 1],
2467 [0, 0, 0, 1]])
2469 >>> A.upper_triangular(2)
2470 Matrix([
2471 [0, 0, 1, 1],
2472 [0, 0, 0, 1],
2473 [0, 0, 0, 0],
2474 [0, 0, 0, 0]])
2476 >>> A.upper_triangular(-1)
2477 Matrix([
2478 [1, 1, 1, 1],
2479 [1, 1, 1, 1],
2480 [0, 1, 1, 1],
2481 [0, 0, 1, 1]])
2483 """
2485 def entry(i, j):
2486 return self[i, j] if i + k <= j else self.zero
2488 return self._new(self.rows, self.cols, entry)
2491 def lower_triangular(self, k=0):
2492 """Return the elements on and below the kth diagonal of a matrix.
2493 If k is not specified then simply returns lower-triangular portion
2494 of a matrix
2496 Examples
2497 ========
2499 >>> from sympy import ones
2500 >>> A = ones(4)
2501 >>> A.lower_triangular()
2502 Matrix([
2503 [1, 0, 0, 0],
2504 [1, 1, 0, 0],
2505 [1, 1, 1, 0],
2506 [1, 1, 1, 1]])
2508 >>> A.lower_triangular(-2)
2509 Matrix([
2510 [0, 0, 0, 0],
2511 [0, 0, 0, 0],
2512 [1, 0, 0, 0],
2513 [1, 1, 0, 0]])
2515 >>> A.lower_triangular(1)
2516 Matrix([
2517 [1, 1, 0, 0],
2518 [1, 1, 1, 0],
2519 [1, 1, 1, 1],
2520 [1, 1, 1, 1]])
2522 """
2524 def entry(i, j):
2525 return self[i, j] if i + k >= j else self.zero
2527 return self._new(self.rows, self.cols, entry)
2531class MatrixArithmetic(MatrixRequired):
2532 """Provides basic matrix arithmetic operations.
2533 Should not be instantiated directly."""
2535 _op_priority = 10.01
2537 def _eval_Abs(self):
2538 return self._new(self.rows, self.cols, lambda i, j: Abs(self[i, j]))
2540 def _eval_add(self, other):
2541 return self._new(self.rows, self.cols,
2542 lambda i, j: self[i, j] + other[i, j])
2544 def _eval_matrix_mul(self, other):
2545 def entry(i, j):
2546 vec = [self[i,k]*other[k,j] for k in range(self.cols)]
2547 try:
2548 return Add(*vec)
2549 except (TypeError, SympifyError):
2550 # Some matrices don't work with `sum` or `Add`
2551 # They don't work with `sum` because `sum` tries to add `0`
2552 # Fall back to a safe way to multiply if the `Add` fails.
2553 return reduce(lambda a, b: a + b, vec)
2555 return self._new(self.rows, other.cols, entry)
2557 def _eval_matrix_mul_elementwise(self, other):
2558 return self._new(self.rows, self.cols, lambda i, j: self[i,j]*other[i,j])
2560 def _eval_matrix_rmul(self, other):
2561 def entry(i, j):
2562 return sum(other[i,k]*self[k,j] for k in range(other.cols))
2563 return self._new(other.rows, self.cols, entry)
2565 def _eval_pow_by_recursion(self, num):
2566 if num == 1:
2567 return self
2569 if num % 2 == 1:
2570 a, b = self, self._eval_pow_by_recursion(num - 1)
2571 else:
2572 a = b = self._eval_pow_by_recursion(num // 2)
2574 return a.multiply(b)
2576 def _eval_pow_by_cayley(self, exp):
2577 from sympy.discrete.recurrences import linrec_coeffs
2578 row = self.shape[0]
2579 p = self.charpoly()
2581 coeffs = (-p).all_coeffs()[1:]
2582 coeffs = linrec_coeffs(coeffs, exp)
2583 new_mat = self.eye(row)
2584 ans = self.zeros(row)
2586 for i in range(row):
2587 ans += coeffs[i]*new_mat
2588 new_mat *= self
2590 return ans
2592 def _eval_pow_by_recursion_dotprodsimp(self, num, prevsimp=None):
2593 if prevsimp is None:
2594 prevsimp = [True]*len(self)
2596 if num == 1:
2597 return self
2599 if num % 2 == 1:
2600 a, b = self, self._eval_pow_by_recursion_dotprodsimp(num - 1,
2601 prevsimp=prevsimp)
2602 else:
2603 a = b = self._eval_pow_by_recursion_dotprodsimp(num // 2,
2604 prevsimp=prevsimp)
2606 m = a.multiply(b, dotprodsimp=False)
2607 lenm = len(m)
2608 elems = [None]*lenm
2610 for i in range(lenm):
2611 if prevsimp[i]:
2612 elems[i], prevsimp[i] = _dotprodsimp(m[i], withsimp=True)
2613 else:
2614 elems[i] = m[i]
2616 return m._new(m.rows, m.cols, elems)
2618 def _eval_scalar_mul(self, other):
2619 return self._new(self.rows, self.cols, lambda i, j: self[i,j]*other)
2621 def _eval_scalar_rmul(self, other):
2622 return self._new(self.rows, self.cols, lambda i, j: other*self[i,j])
2624 def _eval_Mod(self, other):
2625 return self._new(self.rows, self.cols, lambda i, j: Mod(self[i, j], other))
2627 # Python arithmetic functions
2628 def __abs__(self):
2629 """Returns a new matrix with entry-wise absolute values."""
2630 return self._eval_Abs()
2632 @call_highest_priority('__radd__')
2633 def __add__(self, other):
2634 """Return self + other, raising ShapeError if shapes do not match."""
2635 if isinstance(other, NDimArray): # Matrix and array addition is currently not implemented
2636 return NotImplemented
2637 other = _matrixify(other)
2638 # matrix-like objects can have shapes. This is
2639 # our first sanity check.
2640 if hasattr(other, 'shape'):
2641 if self.shape != other.shape:
2642 raise ShapeError("Matrix size mismatch: %s + %s" % (
2643 self.shape, other.shape))
2645 # honest SymPy matrices defer to their class's routine
2646 if getattr(other, 'is_Matrix', False):
2647 # call the highest-priority class's _eval_add
2648 a, b = self, other
2649 if a.__class__ != classof(a, b):
2650 b, a = a, b
2651 return a._eval_add(b)
2652 # Matrix-like objects can be passed to CommonMatrix routines directly.
2653 if getattr(other, 'is_MatrixLike', False):
2654 return MatrixArithmetic._eval_add(self, other)
2656 raise TypeError('cannot add %s and %s' % (type(self), type(other)))
2658 @call_highest_priority('__rtruediv__')
2659 def __truediv__(self, other):
2660 return self * (self.one / other)
2662 @call_highest_priority('__rmatmul__')
2663 def __matmul__(self, other):
2664 other = _matrixify(other)
2665 if not getattr(other, 'is_Matrix', False) and not getattr(other, 'is_MatrixLike', False):
2666 return NotImplemented
2668 return self.__mul__(other)
2670 def __mod__(self, other):
2671 return self.applyfunc(lambda x: x % other)
2673 @call_highest_priority('__rmul__')
2674 def __mul__(self, other):
2675 """Return self*other where other is either a scalar or a matrix
2676 of compatible dimensions.
2678 Examples
2679 ========
2681 >>> from sympy import Matrix
2682 >>> A = Matrix([[1, 2, 3], [4, 5, 6]])
2683 >>> 2*A == A*2 == Matrix([[2, 4, 6], [8, 10, 12]])
2684 True
2685 >>> B = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
2686 >>> A*B
2687 Matrix([
2688 [30, 36, 42],
2689 [66, 81, 96]])
2690 >>> B*A
2691 Traceback (most recent call last):
2692 ...
2693 ShapeError: Matrices size mismatch.
2694 >>>
2696 See Also
2697 ========
2699 matrix_multiply_elementwise
2700 """
2702 return self.multiply(other)
2704 def multiply(self, other, dotprodsimp=None):
2705 """Same as __mul__() but with optional simplification.
2707 Parameters
2708 ==========
2710 dotprodsimp : bool, optional
2711 Specifies whether intermediate term algebraic simplification is used
2712 during matrix multiplications to control expression blowup and thus
2713 speed up calculation. Default is off.
2714 """
2716 isimpbool = _get_intermediate_simp_bool(False, dotprodsimp)
2717 other = _matrixify(other)
2718 # matrix-like objects can have shapes. This is
2719 # our first sanity check. Double check other is not explicitly not a Matrix.
2720 if (hasattr(other, 'shape') and len(other.shape) == 2 and
2721 (getattr(other, 'is_Matrix', True) or
2722 getattr(other, 'is_MatrixLike', True))):
2723 if self.shape[1] != other.shape[0]:
2724 raise ShapeError("Matrix size mismatch: %s * %s." % (
2725 self.shape, other.shape))
2727 # honest SymPy matrices defer to their class's routine
2728 if getattr(other, 'is_Matrix', False):
2729 m = self._eval_matrix_mul(other)
2730 if isimpbool:
2731 return m._new(m.rows, m.cols, [_dotprodsimp(e) for e in m])
2732 return m
2734 # Matrix-like objects can be passed to CommonMatrix routines directly.
2735 if getattr(other, 'is_MatrixLike', False):
2736 return MatrixArithmetic._eval_matrix_mul(self, other)
2738 # if 'other' is not iterable then scalar multiplication.
2739 if not isinstance(other, Iterable):
2740 try:
2741 return self._eval_scalar_mul(other)
2742 except TypeError:
2743 pass
2745 return NotImplemented
2747 def multiply_elementwise(self, other):
2748 """Return the Hadamard product (elementwise product) of A and B
2750 Examples
2751 ========
2753 >>> from sympy import Matrix
2754 >>> A = Matrix([[0, 1, 2], [3, 4, 5]])
2755 >>> B = Matrix([[1, 10, 100], [100, 10, 1]])
2756 >>> A.multiply_elementwise(B)
2757 Matrix([
2758 [ 0, 10, 200],
2759 [300, 40, 5]])
2761 See Also
2762 ========
2764 sympy.matrices.matrices.MatrixBase.cross
2765 sympy.matrices.matrices.MatrixBase.dot
2766 multiply
2767 """
2768 if self.shape != other.shape:
2769 raise ShapeError("Matrix shapes must agree {} != {}".format(self.shape, other.shape))
2771 return self._eval_matrix_mul_elementwise(other)
2773 def __neg__(self):
2774 return self._eval_scalar_mul(-1)
2776 @call_highest_priority('__rpow__')
2777 def __pow__(self, exp):
2778 """Return self**exp a scalar or symbol."""
2780 return self.pow(exp)
2783 def pow(self, exp, method=None):
2784 r"""Return self**exp a scalar or symbol.
2786 Parameters
2787 ==========
2789 method : multiply, mulsimp, jordan, cayley
2790 If multiply then it returns exponentiation using recursion.
2791 If jordan then Jordan form exponentiation will be used.
2792 If cayley then the exponentiation is done using Cayley-Hamilton
2793 theorem.
2794 If mulsimp then the exponentiation is done using recursion
2795 with dotprodsimp. This specifies whether intermediate term
2796 algebraic simplification is used during naive matrix power to
2797 control expression blowup and thus speed up calculation.
2798 If None, then it heuristically decides which method to use.
2800 """
2802 if method is not None and method not in ['multiply', 'mulsimp', 'jordan', 'cayley']:
2803 raise TypeError('No such method')
2804 if self.rows != self.cols:
2805 raise NonSquareMatrixError()
2806 a = self
2807 jordan_pow = getattr(a, '_matrix_pow_by_jordan_blocks', None)
2808 exp = sympify(exp)
2810 if exp.is_zero:
2811 return a._new(a.rows, a.cols, lambda i, j: int(i == j))
2812 if exp == 1:
2813 return a
2815 diagonal = getattr(a, 'is_diagonal', None)
2816 if diagonal is not None and diagonal():
2817 return a._new(a.rows, a.cols, lambda i, j: a[i,j]**exp if i == j else 0)
2819 if exp.is_Number and exp % 1 == 0:
2820 if a.rows == 1:
2821 return a._new([[a[0]**exp]])
2822 if exp < 0:
2823 exp = -exp
2824 a = a.inv()
2825 # When certain conditions are met,
2826 # Jordan block algorithm is faster than
2827 # computation by recursion.
2828 if method == 'jordan':
2829 try:
2830 return jordan_pow(exp)
2831 except MatrixError:
2832 if method == 'jordan':
2833 raise
2835 elif method == 'cayley':
2836 if not exp.is_Number or exp % 1 != 0:
2837 raise ValueError("cayley method is only valid for integer powers")
2838 return a._eval_pow_by_cayley(exp)
2840 elif method == "mulsimp":
2841 if not exp.is_Number or exp % 1 != 0:
2842 raise ValueError("mulsimp method is only valid for integer powers")
2843 return a._eval_pow_by_recursion_dotprodsimp(exp)
2845 elif method == "multiply":
2846 if not exp.is_Number or exp % 1 != 0:
2847 raise ValueError("multiply method is only valid for integer powers")
2848 return a._eval_pow_by_recursion(exp)
2850 elif method is None and exp.is_Number and exp % 1 == 0:
2851 # Decide heuristically which method to apply
2852 if a.rows == 2 and exp > 100000:
2853 return jordan_pow(exp)
2854 elif _get_intermediate_simp_bool(True, None):
2855 return a._eval_pow_by_recursion_dotprodsimp(exp)
2856 elif exp > 10000:
2857 return a._eval_pow_by_cayley(exp)
2858 else:
2859 return a._eval_pow_by_recursion(exp)
2861 if jordan_pow:
2862 try:
2863 return jordan_pow(exp)
2864 except NonInvertibleMatrixError:
2865 # Raised by jordan_pow on zero determinant matrix unless exp is
2866 # definitely known to be a non-negative integer.
2867 # Here we raise if n is definitely not a non-negative integer
2868 # but otherwise we can leave this as an unevaluated MatPow.
2869 if exp.is_integer is False or exp.is_nonnegative is False:
2870 raise
2872 from sympy.matrices.expressions import MatPow
2873 return MatPow(a, exp)
2875 @call_highest_priority('__add__')
2876 def __radd__(self, other):
2877 return self + other
2879 @call_highest_priority('__matmul__')
2880 def __rmatmul__(self, other):
2881 other = _matrixify(other)
2882 if not getattr(other, 'is_Matrix', False) and not getattr(other, 'is_MatrixLike', False):
2883 return NotImplemented
2885 return self.__rmul__(other)
2887 @call_highest_priority('__mul__')
2888 def __rmul__(self, other):
2889 return self.rmultiply(other)
2891 def rmultiply(self, other, dotprodsimp=None):
2892 """Same as __rmul__() but with optional simplification.
2894 Parameters
2895 ==========
2897 dotprodsimp : bool, optional
2898 Specifies whether intermediate term algebraic simplification is used
2899 during matrix multiplications to control expression blowup and thus
2900 speed up calculation. Default is off.
2901 """
2902 isimpbool = _get_intermediate_simp_bool(False, dotprodsimp)
2903 other = _matrixify(other)
2904 # matrix-like objects can have shapes. This is
2905 # our first sanity check. Double check other is not explicitly not a Matrix.
2906 if (hasattr(other, 'shape') and len(other.shape) == 2 and
2907 (getattr(other, 'is_Matrix', True) or
2908 getattr(other, 'is_MatrixLike', True))):
2909 if self.shape[0] != other.shape[1]:
2910 raise ShapeError("Matrix size mismatch.")
2912 # honest SymPy matrices defer to their class's routine
2913 if getattr(other, 'is_Matrix', False):
2914 m = self._eval_matrix_rmul(other)
2915 if isimpbool:
2916 return m._new(m.rows, m.cols, [_dotprodsimp(e) for e in m])
2917 return m
2918 # Matrix-like objects can be passed to CommonMatrix routines directly.
2919 if getattr(other, 'is_MatrixLike', False):
2920 return MatrixArithmetic._eval_matrix_rmul(self, other)
2922 # if 'other' is not iterable then scalar multiplication.
2923 if not isinstance(other, Iterable):
2924 try:
2925 return self._eval_scalar_rmul(other)
2926 except TypeError:
2927 pass
2929 return NotImplemented
2931 @call_highest_priority('__sub__')
2932 def __rsub__(self, a):
2933 return (-self) + a
2935 @call_highest_priority('__rsub__')
2936 def __sub__(self, a):
2937 return self + (-a)
2939class MatrixCommon(MatrixArithmetic, MatrixOperations, MatrixProperties,
2940 MatrixSpecial, MatrixShaping):
2941 """All common matrix operations including basic arithmetic, shaping,
2942 and special matrices like `zeros`, and `eye`."""
2943 _diff_wrt = True # type: bool
2946class _MinimalMatrix:
2947 """Class providing the minimum functionality
2948 for a matrix-like object and implementing every method
2949 required for a `MatrixRequired`. This class does not have everything
2950 needed to become a full-fledged SymPy object, but it will satisfy the
2951 requirements of anything inheriting from `MatrixRequired`. If you wish
2952 to make a specialized matrix type, make sure to implement these
2953 methods and properties with the exception of `__init__` and `__repr__`
2954 which are included for convenience."""
2956 is_MatrixLike = True
2957 _sympify = staticmethod(sympify)
2958 _class_priority = 3
2959 zero = S.Zero
2960 one = S.One
2962 is_Matrix = True
2963 is_MatrixExpr = False
2965 @classmethod
2966 def _new(cls, *args, **kwargs):
2967 return cls(*args, **kwargs)
2969 def __init__(self, rows, cols=None, mat=None, copy=False):
2970 if isfunction(mat):
2971 # if we passed in a function, use that to populate the indices
2972 mat = [mat(i, j) for i in range(rows) for j in range(cols)]
2973 if cols is None and mat is None:
2974 mat = rows
2975 rows, cols = getattr(mat, 'shape', (rows, cols))
2976 try:
2977 # if we passed in a list of lists, flatten it and set the size
2978 if cols is None and mat is None:
2979 mat = rows
2980 cols = len(mat[0])
2981 rows = len(mat)
2982 mat = [x for l in mat for x in l]
2983 except (IndexError, TypeError):
2984 pass
2985 self.mat = tuple(self._sympify(x) for x in mat)
2986 self.rows, self.cols = rows, cols
2987 if self.rows is None or self.cols is None:
2988 raise NotImplementedError("Cannot initialize matrix with given parameters")
2990 def __getitem__(self, key):
2991 def _normalize_slices(row_slice, col_slice):
2992 """Ensure that row_slice and col_slice do not have
2993 `None` in their arguments. Any integers are converted
2994 to slices of length 1"""
2995 if not isinstance(row_slice, slice):
2996 row_slice = slice(row_slice, row_slice + 1, None)
2997 row_slice = slice(*row_slice.indices(self.rows))
2999 if not isinstance(col_slice, slice):
3000 col_slice = slice(col_slice, col_slice + 1, None)
3001 col_slice = slice(*col_slice.indices(self.cols))
3003 return (row_slice, col_slice)
3005 def _coord_to_index(i, j):
3006 """Return the index in _mat corresponding
3007 to the (i,j) position in the matrix. """
3008 return i * self.cols + j
3010 if isinstance(key, tuple):
3011 i, j = key
3012 if isinstance(i, slice) or isinstance(j, slice):
3013 # if the coordinates are not slices, make them so
3014 # and expand the slices so they don't contain `None`
3015 i, j = _normalize_slices(i, j)
3017 rowsList, colsList = list(range(self.rows))[i], \
3018 list(range(self.cols))[j]
3019 indices = (i * self.cols + j for i in rowsList for j in
3020 colsList)
3021 return self._new(len(rowsList), len(colsList),
3022 [self.mat[i] for i in indices])
3024 # if the key is a tuple of ints, change
3025 # it to an array index
3026 key = _coord_to_index(i, j)
3027 return self.mat[key]
3029 def __eq__(self, other):
3030 try:
3031 classof(self, other)
3032 except TypeError:
3033 return False
3034 return (
3035 self.shape == other.shape and list(self) == list(other))
3037 def __len__(self):
3038 return self.rows*self.cols
3040 def __repr__(self):
3041 return "_MinimalMatrix({}, {}, {})".format(self.rows, self.cols,
3042 self.mat)
3044 @property
3045 def shape(self):
3046 return (self.rows, self.cols)
3049class _CastableMatrix: # this is needed here ONLY FOR TESTS.
3050 def as_mutable(self):
3051 return self
3053 def as_immutable(self):
3054 return self
3057class _MatrixWrapper:
3058 """Wrapper class providing the minimum functionality for a matrix-like
3059 object: .rows, .cols, .shape, indexability, and iterability. CommonMatrix
3060 math operations should work on matrix-like objects. This one is intended for
3061 matrix-like objects which use the same indexing format as SymPy with respect
3062 to returning matrix elements instead of rows for non-tuple indexes.
3063 """
3065 is_Matrix = False # needs to be here because of __getattr__
3066 is_MatrixLike = True
3068 def __init__(self, mat, shape):
3069 self.mat = mat
3070 self.shape = shape
3071 self.rows, self.cols = shape
3073 def __getitem__(self, key):
3074 if isinstance(key, tuple):
3075 return sympify(self.mat.__getitem__(key))
3077 return sympify(self.mat.__getitem__((key // self.rows, key % self.cols)))
3079 def __iter__(self): # supports numpy.matrix and numpy.array
3080 mat = self.mat
3081 cols = self.cols
3083 return iter(sympify(mat[r, c]) for r in range(self.rows) for c in range(cols))
3086class MatrixKind(Kind):
3087 """
3088 Kind for all matrices in SymPy.
3090 Basic class for this kind is ``MatrixBase`` and ``MatrixExpr``,
3091 but any expression representing the matrix can have this.
3093 Parameters
3094 ==========
3096 element_kind : Kind
3097 Kind of the element. Default is
3098 :class:`sympy.core.kind.NumberKind`,
3099 which means that the matrix contains only numbers.
3101 Examples
3102 ========
3104 Any instance of matrix class has ``MatrixKind``:
3106 >>> from sympy import MatrixSymbol
3107 >>> A = MatrixSymbol('A', 2,2)
3108 >>> A.kind
3109 MatrixKind(NumberKind)
3111 Although expression representing a matrix may be not instance of
3112 matrix class, it will have ``MatrixKind`` as well:
3114 >>> from sympy import MatrixExpr, Integral
3115 >>> from sympy.abc import x
3116 >>> intM = Integral(A, x)
3117 >>> isinstance(intM, MatrixExpr)
3118 False
3119 >>> intM.kind
3120 MatrixKind(NumberKind)
3122 Use ``isinstance()`` to check for ``MatrixKind`` without specifying
3123 the element kind. Use ``is`` with specifying the element kind:
3125 >>> from sympy import Matrix
3126 >>> from sympy.core import NumberKind
3127 >>> from sympy.matrices import MatrixKind
3128 >>> M = Matrix([1, 2])
3129 >>> isinstance(M.kind, MatrixKind)
3130 True
3131 >>> M.kind is MatrixKind(NumberKind)
3132 True
3134 See Also
3135 ========
3137 sympy.core.kind.NumberKind
3138 sympy.core.kind.UndefinedKind
3139 sympy.core.containers.TupleKind
3140 sympy.sets.sets.SetKind
3142 """
3143 def __new__(cls, element_kind=NumberKind):
3144 obj = super().__new__(cls, element_kind)
3145 obj.element_kind = element_kind
3146 return obj
3148 def __repr__(self):
3149 return "MatrixKind(%s)" % self.element_kind
3152def _matrixify(mat):
3153 """If `mat` is a Matrix or is matrix-like,
3154 return a Matrix or MatrixWrapper object. Otherwise
3155 `mat` is passed through without modification."""
3157 if getattr(mat, 'is_Matrix', False) or getattr(mat, 'is_MatrixLike', False):
3158 return mat
3160 if not(getattr(mat, 'is_Matrix', True) or getattr(mat, 'is_MatrixLike', True)):
3161 return mat
3163 shape = None
3165 if hasattr(mat, 'shape'): # numpy, scipy.sparse
3166 if len(mat.shape) == 2:
3167 shape = mat.shape
3168 elif hasattr(mat, 'rows') and hasattr(mat, 'cols'): # mpmath
3169 shape = (mat.rows, mat.cols)
3171 if shape:
3172 return _MatrixWrapper(mat, shape)
3174 return mat
3177def a2idx(j, n=None):
3178 """Return integer after making positive and validating against n."""
3179 if not isinstance(j, int):
3180 jindex = getattr(j, '__index__', None)
3181 if jindex is not None:
3182 j = jindex()
3183 else:
3184 raise IndexError("Invalid index a[%r]" % (j,))
3185 if n is not None:
3186 if j < 0:
3187 j += n
3188 if not (j >= 0 and j < n):
3189 raise IndexError("Index out of range: a[%s]" % (j,))
3190 return int(j)
3193def classof(A, B):
3194 """
3195 Get the type of the result when combining matrices of different types.
3197 Currently the strategy is that immutability is contagious.
3199 Examples
3200 ========
3202 >>> from sympy import Matrix, ImmutableMatrix
3203 >>> from sympy.matrices.common import classof
3204 >>> M = Matrix([[1, 2], [3, 4]]) # a Mutable Matrix
3205 >>> IM = ImmutableMatrix([[1, 2], [3, 4]])
3206 >>> classof(M, IM)
3207 <class 'sympy.matrices.immutable.ImmutableDenseMatrix'>
3208 """
3209 priority_A = getattr(A, '_class_priority', None)
3210 priority_B = getattr(B, '_class_priority', None)
3211 if None not in (priority_A, priority_B):
3212 if A._class_priority > B._class_priority:
3213 return A.__class__
3214 else:
3215 return B.__class__
3217 try:
3218 import numpy
3219 except ImportError:
3220 pass
3221 else:
3222 if isinstance(A, numpy.ndarray):
3223 return B.__class__
3224 if isinstance(B, numpy.ndarray):
3225 return A.__class__
3227 raise TypeError("Incompatible classes %s, %s" % (A.__class__, B.__class__))