Coverage for /usr/lib/python3/dist-packages/sympy/matrices/eigen.py: 11%
475 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1from types import FunctionType
2from collections import Counter
4from mpmath import mp, workprec
5from mpmath.libmp.libmpf import prec_to_dps
7from sympy.core.sorting import default_sort_key
8from sympy.core.evalf import DEFAULT_MAXPREC, PrecisionExhausted
9from sympy.core.logic import fuzzy_and, fuzzy_or
10from sympy.core.numbers import Float
11from sympy.core.sympify import _sympify
12from sympy.functions.elementary.miscellaneous import sqrt
13from sympy.polys import roots, CRootOf, ZZ, QQ, EX
14from sympy.polys.matrices import DomainMatrix
15from sympy.polys.matrices.eigen import dom_eigenvects, dom_eigenvects_to_sympy
16from sympy.polys.polytools import gcd
18from .common import MatrixError, NonSquareMatrixError
19from .determinant import _find_reasonable_pivot
21from .utilities import _iszero, _simplify
24def _eigenvals_eigenvects_mpmath(M):
25 norm2 = lambda v: mp.sqrt(sum(i**2 for i in v))
27 v1 = None
28 prec = max([x._prec for x in M.atoms(Float)])
29 eps = 2**-prec
31 while prec < DEFAULT_MAXPREC:
32 with workprec(prec):
33 A = mp.matrix(M.evalf(n=prec_to_dps(prec)))
34 E, ER = mp.eig(A)
35 v2 = norm2([i for e in E for i in (mp.re(e), mp.im(e))])
36 if v1 is not None and mp.fabs(v1 - v2) < eps:
37 return E, ER
38 v1 = v2
39 prec *= 2
41 # we get here because the next step would have taken us
42 # past MAXPREC or because we never took a step; in case
43 # of the latter, we refuse to send back a solution since
44 # it would not have been verified; we also resist taking
45 # a small step to arrive exactly at MAXPREC since then
46 # the two calculations might be artificially close.
47 raise PrecisionExhausted
50def _eigenvals_mpmath(M, multiple=False):
51 """Compute eigenvalues using mpmath"""
52 E, _ = _eigenvals_eigenvects_mpmath(M)
53 result = [_sympify(x) for x in E]
54 if multiple:
55 return result
56 return dict(Counter(result))
59def _eigenvects_mpmath(M):
60 E, ER = _eigenvals_eigenvects_mpmath(M)
61 result = []
62 for i in range(M.rows):
63 eigenval = _sympify(E[i])
64 eigenvect = _sympify(ER[:, i])
65 result.append((eigenval, 1, [eigenvect]))
67 return result
70# This function is a candidate for caching if it gets implemented for matrices.
71def _eigenvals(
72 M, error_when_incomplete=True, *, simplify=False, multiple=False,
73 rational=False, **flags):
74 r"""Compute eigenvalues of the matrix.
76 Parameters
77 ==========
79 error_when_incomplete : bool, optional
80 If it is set to ``True``, it will raise an error if not all
81 eigenvalues are computed. This is caused by ``roots`` not returning
82 a full list of eigenvalues.
84 simplify : bool or function, optional
85 If it is set to ``True``, it attempts to return the most
86 simplified form of expressions returned by applying default
87 simplification method in every routine.
89 If it is set to ``False``, it will skip simplification in this
90 particular routine to save computation resources.
92 If a function is passed to, it will attempt to apply
93 the particular function as simplification method.
95 rational : bool, optional
96 If it is set to ``True``, every floating point numbers would be
97 replaced with rationals before computation. It can solve some
98 issues of ``roots`` routine not working well with floats.
100 multiple : bool, optional
101 If it is set to ``True``, the result will be in the form of a
102 list.
104 If it is set to ``False``, the result will be in the form of a
105 dictionary.
107 Returns
108 =======
110 eigs : list or dict
111 Eigenvalues of a matrix. The return format would be specified by
112 the key ``multiple``.
114 Raises
115 ======
117 MatrixError
118 If not enough roots had got computed.
120 NonSquareMatrixError
121 If attempted to compute eigenvalues from a non-square matrix.
123 Examples
124 ========
126 >>> from sympy import Matrix
127 >>> M = Matrix(3, 3, [0, 1, 1, 1, 0, 0, 1, 1, 1])
128 >>> M.eigenvals()
129 {-1: 1, 0: 1, 2: 1}
131 See Also
132 ========
134 MatrixDeterminant.charpoly
135 eigenvects
137 Notes
138 =====
140 Eigenvalues of a matrix $A$ can be computed by solving a matrix
141 equation $\det(A - \lambda I) = 0$
143 It's not always possible to return radical solutions for
144 eigenvalues for matrices larger than $4, 4$ shape due to
145 Abel-Ruffini theorem.
147 If there is no radical solution is found for the eigenvalue,
148 it may return eigenvalues in the form of
149 :class:`sympy.polys.rootoftools.ComplexRootOf`.
150 """
151 if not M:
152 if multiple:
153 return []
154 return {}
156 if not M.is_square:
157 raise NonSquareMatrixError("{} must be a square matrix.".format(M))
159 if M._rep.domain not in (ZZ, QQ):
160 # Skip this check for ZZ/QQ because it can be slow
161 if all(x.is_number for x in M) and M.has(Float):
162 return _eigenvals_mpmath(M, multiple=multiple)
164 if rational:
165 from sympy.simplify import nsimplify
166 M = M.applyfunc(
167 lambda x: nsimplify(x, rational=True) if x.has(Float) else x)
169 if multiple:
170 return _eigenvals_list(
171 M, error_when_incomplete=error_when_incomplete, simplify=simplify,
172 **flags)
173 return _eigenvals_dict(
174 M, error_when_incomplete=error_when_incomplete, simplify=simplify,
175 **flags)
178eigenvals_error_message = \
179"It is not always possible to express the eigenvalues of a matrix " + \
180"of size 5x5 or higher in radicals. " + \
181"We have CRootOf, but domains other than the rationals are not " + \
182"currently supported. " + \
183"If there are no symbols in the matrix, " + \
184"it should still be possible to compute numeric approximations " + \
185"of the eigenvalues using " + \
186"M.evalf().eigenvals() or M.charpoly().nroots()."
189def _eigenvals_list(
190 M, error_when_incomplete=True, simplify=False, **flags):
191 iblocks = M.strongly_connected_components()
192 all_eigs = []
193 is_dom = M._rep.domain in (ZZ, QQ)
194 for b in iblocks:
196 # Fast path for a 1x1 block:
197 if is_dom and len(b) == 1:
198 index = b[0]
199 val = M[index, index]
200 all_eigs.append(val)
201 continue
203 block = M[b, b]
205 if isinstance(simplify, FunctionType):
206 charpoly = block.charpoly(simplify=simplify)
207 else:
208 charpoly = block.charpoly()
210 eigs = roots(charpoly, multiple=True, **flags)
212 if len(eigs) != block.rows:
213 degree = int(charpoly.degree())
214 f = charpoly.as_expr()
215 x = charpoly.gen
216 try:
217 eigs = [CRootOf(f, x, idx) for idx in range(degree)]
218 except NotImplementedError:
219 if error_when_incomplete:
220 raise MatrixError(eigenvals_error_message)
221 else:
222 eigs = []
224 all_eigs += eigs
226 if not simplify:
227 return all_eigs
228 if not isinstance(simplify, FunctionType):
229 simplify = _simplify
230 return [simplify(value) for value in all_eigs]
233def _eigenvals_dict(
234 M, error_when_incomplete=True, simplify=False, **flags):
235 iblocks = M.strongly_connected_components()
236 all_eigs = {}
237 is_dom = M._rep.domain in (ZZ, QQ)
238 for b in iblocks:
240 # Fast path for a 1x1 block:
241 if is_dom and len(b) == 1:
242 index = b[0]
243 val = M[index, index]
244 all_eigs[val] = all_eigs.get(val, 0) + 1
245 continue
247 block = M[b, b]
249 if isinstance(simplify, FunctionType):
250 charpoly = block.charpoly(simplify=simplify)
251 else:
252 charpoly = block.charpoly()
254 eigs = roots(charpoly, multiple=False, **flags)
256 if sum(eigs.values()) != block.rows:
257 degree = int(charpoly.degree())
258 f = charpoly.as_expr()
259 x = charpoly.gen
260 try:
261 eigs = {CRootOf(f, x, idx): 1 for idx in range(degree)}
262 except NotImplementedError:
263 if error_when_incomplete:
264 raise MatrixError(eigenvals_error_message)
265 else:
266 eigs = {}
268 for k, v in eigs.items():
269 if k in all_eigs:
270 all_eigs[k] += v
271 else:
272 all_eigs[k] = v
274 if not simplify:
275 return all_eigs
276 if not isinstance(simplify, FunctionType):
277 simplify = _simplify
278 return {simplify(key): value for key, value in all_eigs.items()}
281def _eigenspace(M, eigenval, iszerofunc=_iszero, simplify=False):
282 """Get a basis for the eigenspace for a particular eigenvalue"""
283 m = M - M.eye(M.rows) * eigenval
284 ret = m.nullspace(iszerofunc=iszerofunc)
286 # The nullspace for a real eigenvalue should be non-trivial.
287 # If we didn't find an eigenvector, try once more a little harder
288 if len(ret) == 0 and simplify:
289 ret = m.nullspace(iszerofunc=iszerofunc, simplify=True)
290 if len(ret) == 0:
291 raise NotImplementedError(
292 "Can't evaluate eigenvector for eigenvalue {}".format(eigenval))
293 return ret
296def _eigenvects_DOM(M, **kwargs):
297 DOM = DomainMatrix.from_Matrix(M, field=True, extension=True)
298 DOM = DOM.to_dense()
300 if DOM.domain != EX:
301 rational, algebraic = dom_eigenvects(DOM)
302 eigenvects = dom_eigenvects_to_sympy(
303 rational, algebraic, M.__class__, **kwargs)
304 eigenvects = sorted(eigenvects, key=lambda x: default_sort_key(x[0]))
306 return eigenvects
307 return None
310def _eigenvects_sympy(M, iszerofunc, simplify=True, **flags):
311 eigenvals = M.eigenvals(rational=False, **flags)
313 # Make sure that we have all roots in radical form
314 for x in eigenvals:
315 if x.has(CRootOf):
316 raise MatrixError(
317 "Eigenvector computation is not implemented if the matrix have "
318 "eigenvalues in CRootOf form")
320 eigenvals = sorted(eigenvals.items(), key=default_sort_key)
321 ret = []
322 for val, mult in eigenvals:
323 vects = _eigenspace(M, val, iszerofunc=iszerofunc, simplify=simplify)
324 ret.append((val, mult, vects))
325 return ret
328# This functions is a candidate for caching if it gets implemented for matrices.
329def _eigenvects(M, error_when_incomplete=True, iszerofunc=_iszero, *, chop=False, **flags):
330 """Compute eigenvectors of the matrix.
332 Parameters
333 ==========
335 error_when_incomplete : bool, optional
336 Raise an error when not all eigenvalues are computed. This is
337 caused by ``roots`` not returning a full list of eigenvalues.
339 iszerofunc : function, optional
340 Specifies a zero testing function to be used in ``rref``.
342 Default value is ``_iszero``, which uses SymPy's naive and fast
343 default assumption handler.
345 It can also accept any user-specified zero testing function, if it
346 is formatted as a function which accepts a single symbolic argument
347 and returns ``True`` if it is tested as zero and ``False`` if it
348 is tested as non-zero, and ``None`` if it is undecidable.
350 simplify : bool or function, optional
351 If ``True``, ``as_content_primitive()`` will be used to tidy up
352 normalization artifacts.
354 It will also be used by the ``nullspace`` routine.
356 chop : bool or positive number, optional
357 If the matrix contains any Floats, they will be changed to Rationals
358 for computation purposes, but the answers will be returned after
359 being evaluated with evalf. The ``chop`` flag is passed to ``evalf``.
360 When ``chop=True`` a default precision will be used; a number will
361 be interpreted as the desired level of precision.
363 Returns
364 =======
366 ret : [(eigenval, multiplicity, eigenspace), ...]
367 A ragged list containing tuples of data obtained by ``eigenvals``
368 and ``nullspace``.
370 ``eigenspace`` is a list containing the ``eigenvector`` for each
371 eigenvalue.
373 ``eigenvector`` is a vector in the form of a ``Matrix``. e.g.
374 a vector of length 3 is returned as ``Matrix([a_1, a_2, a_3])``.
376 Raises
377 ======
379 NotImplementedError
380 If failed to compute nullspace.
382 Examples
383 ========
385 >>> from sympy import Matrix
386 >>> M = Matrix(3, 3, [0, 1, 1, 1, 0, 0, 1, 1, 1])
387 >>> M.eigenvects()
388 [(-1, 1, [Matrix([
389 [-1],
390 [ 1],
391 [ 0]])]), (0, 1, [Matrix([
392 [ 0],
393 [-1],
394 [ 1]])]), (2, 1, [Matrix([
395 [2/3],
396 [1/3],
397 [ 1]])])]
399 See Also
400 ========
402 eigenvals
403 MatrixSubspaces.nullspace
404 """
405 simplify = flags.get('simplify', True)
406 primitive = flags.get('simplify', False)
407 flags.pop('simplify', None) # remove this if it's there
408 flags.pop('multiple', None) # remove this if it's there
410 if not isinstance(simplify, FunctionType):
411 simpfunc = _simplify if simplify else lambda x: x
413 has_floats = M.has(Float)
414 if has_floats:
415 if all(x.is_number for x in M):
416 return _eigenvects_mpmath(M)
417 from sympy.simplify import nsimplify
418 M = M.applyfunc(lambda x: nsimplify(x, rational=True))
420 ret = _eigenvects_DOM(M)
421 if ret is None:
422 ret = _eigenvects_sympy(M, iszerofunc, simplify=simplify, **flags)
424 if primitive:
425 # if the primitive flag is set, get rid of any common
426 # integer denominators
427 def denom_clean(l):
428 return [(v / gcd(list(v))).applyfunc(simpfunc) for v in l]
430 ret = [(val, mult, denom_clean(es)) for val, mult, es in ret]
432 if has_floats:
433 # if we had floats to start with, turn the eigenvectors to floats
434 ret = [(val.evalf(chop=chop), mult, [v.evalf(chop=chop) for v in es])
435 for val, mult, es in ret]
437 return ret
440def _is_diagonalizable_with_eigen(M, reals_only=False):
441 """See _is_diagonalizable. This function returns the bool along with the
442 eigenvectors to avoid calculating them again in functions like
443 ``diagonalize``."""
445 if not M.is_square:
446 return False, []
448 eigenvecs = M.eigenvects(simplify=True)
450 for val, mult, basis in eigenvecs:
451 if reals_only and not val.is_real: # if we have a complex eigenvalue
452 return False, eigenvecs
454 if mult != len(basis): # if the geometric multiplicity doesn't equal the algebraic
455 return False, eigenvecs
457 return True, eigenvecs
459def _is_diagonalizable(M, reals_only=False, **kwargs):
460 """Returns ``True`` if a matrix is diagonalizable.
462 Parameters
463 ==========
465 reals_only : bool, optional
466 If ``True``, it tests whether the matrix can be diagonalized
467 to contain only real numbers on the diagonal.
470 If ``False``, it tests whether the matrix can be diagonalized
471 at all, even with numbers that may not be real.
473 Examples
474 ========
476 Example of a diagonalizable matrix:
478 >>> from sympy import Matrix
479 >>> M = Matrix([[1, 2, 0], [0, 3, 0], [2, -4, 2]])
480 >>> M.is_diagonalizable()
481 True
483 Example of a non-diagonalizable matrix:
485 >>> M = Matrix([[0, 1], [0, 0]])
486 >>> M.is_diagonalizable()
487 False
489 Example of a matrix that is diagonalized in terms of non-real entries:
491 >>> M = Matrix([[0, 1], [-1, 0]])
492 >>> M.is_diagonalizable(reals_only=False)
493 True
494 >>> M.is_diagonalizable(reals_only=True)
495 False
497 See Also
498 ========
500 is_diagonal
501 diagonalize
502 """
503 if not M.is_square:
504 return False
506 if all(e.is_real for e in M) and M.is_symmetric():
507 return True
509 if all(e.is_complex for e in M) and M.is_hermitian:
510 return True
512 return _is_diagonalizable_with_eigen(M, reals_only=reals_only)[0]
515#G&VL, Matrix Computations, Algo 5.4.2
516def _householder_vector(x):
517 if not x.cols == 1:
518 raise ValueError("Input must be a column matrix")
519 v = x.copy()
520 v_plus = x.copy()
521 v_minus = x.copy()
522 q = x[0, 0] / abs(x[0, 0])
523 norm_x = x.norm()
524 v_plus[0, 0] = x[0, 0] + q * norm_x
525 v_minus[0, 0] = x[0, 0] - q * norm_x
526 if x[1:, 0].norm() == 0:
527 bet = 0
528 v[0, 0] = 1
529 else:
530 if v_plus.norm() <= v_minus.norm():
531 v = v_plus
532 else:
533 v = v_minus
534 v = v / v[0]
535 bet = 2 / (v.norm() ** 2)
536 return v, bet
539def _bidiagonal_decmp_hholder(M):
540 m = M.rows
541 n = M.cols
542 A = M.as_mutable()
543 U, V = A.eye(m), A.eye(n)
544 for i in range(min(m, n)):
545 v, bet = _householder_vector(A[i:, i])
546 hh_mat = A.eye(m - i) - bet * v * v.H
547 A[i:, i:] = hh_mat * A[i:, i:]
548 temp = A.eye(m)
549 temp[i:, i:] = hh_mat
550 U = U * temp
551 if i + 1 <= n - 2:
552 v, bet = _householder_vector(A[i, i+1:].T)
553 hh_mat = A.eye(n - i - 1) - bet * v * v.H
554 A[i:, i+1:] = A[i:, i+1:] * hh_mat
555 temp = A.eye(n)
556 temp[i+1:, i+1:] = hh_mat
557 V = temp * V
558 return U, A, V
561def _eval_bidiag_hholder(M):
562 m = M.rows
563 n = M.cols
564 A = M.as_mutable()
565 for i in range(min(m, n)):
566 v, bet = _householder_vector(A[i:, i])
567 hh_mat = A.eye(m-i) - bet * v * v.H
568 A[i:, i:] = hh_mat * A[i:, i:]
569 if i + 1 <= n - 2:
570 v, bet = _householder_vector(A[i, i+1:].T)
571 hh_mat = A.eye(n - i - 1) - bet * v * v.H
572 A[i:, i+1:] = A[i:, i+1:] * hh_mat
573 return A
576def _bidiagonal_decomposition(M, upper=True):
577 """
578 Returns $(U,B,V.H)$ for
580 $$A = UBV^{H}$$
582 where $A$ is the input matrix, and $B$ is its Bidiagonalized form
584 Note: Bidiagonal Computation can hang for symbolic matrices.
586 Parameters
587 ==========
589 upper : bool. Whether to do upper bidiagnalization or lower.
590 True for upper and False for lower.
592 References
593 ==========
595 .. [1] Algorithm 5.4.2, Matrix computations by Golub and Van Loan, 4th edition
596 .. [2] Complex Matrix Bidiagonalization, https://github.com/vslobody/Householder-Bidiagonalization
598 """
600 if not isinstance(upper, bool):
601 raise ValueError("upper must be a boolean")
603 if upper:
604 return _bidiagonal_decmp_hholder(M)
606 X = _bidiagonal_decmp_hholder(M.H)
607 return X[2].H, X[1].H, X[0].H
610def _bidiagonalize(M, upper=True):
611 """
612 Returns $B$, the Bidiagonalized form of the input matrix.
614 Note: Bidiagonal Computation can hang for symbolic matrices.
616 Parameters
617 ==========
619 upper : bool. Whether to do upper bidiagnalization or lower.
620 True for upper and False for lower.
622 References
623 ==========
625 .. [1] Algorithm 5.4.2, Matrix computations by Golub and Van Loan, 4th edition
626 .. [2] Complex Matrix Bidiagonalization : https://github.com/vslobody/Householder-Bidiagonalization
628 """
630 if not isinstance(upper, bool):
631 raise ValueError("upper must be a boolean")
633 if upper:
634 return _eval_bidiag_hholder(M)
635 return _eval_bidiag_hholder(M.H).H
638def _diagonalize(M, reals_only=False, sort=False, normalize=False):
639 """
640 Return (P, D), where D is diagonal and
642 D = P^-1 * M * P
644 where M is current matrix.
646 Parameters
647 ==========
649 reals_only : bool. Whether to throw an error if complex numbers are need
650 to diagonalize. (Default: False)
652 sort : bool. Sort the eigenvalues along the diagonal. (Default: False)
654 normalize : bool. If True, normalize the columns of P. (Default: False)
656 Examples
657 ========
659 >>> from sympy import Matrix
660 >>> M = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
661 >>> M
662 Matrix([
663 [1, 2, 0],
664 [0, 3, 0],
665 [2, -4, 2]])
666 >>> (P, D) = M.diagonalize()
667 >>> D
668 Matrix([
669 [1, 0, 0],
670 [0, 2, 0],
671 [0, 0, 3]])
672 >>> P
673 Matrix([
674 [-1, 0, -1],
675 [ 0, 0, -1],
676 [ 2, 1, 2]])
677 >>> P.inv() * M * P
678 Matrix([
679 [1, 0, 0],
680 [0, 2, 0],
681 [0, 0, 3]])
683 See Also
684 ========
686 is_diagonal
687 is_diagonalizable
688 """
690 if not M.is_square:
691 raise NonSquareMatrixError()
693 is_diagonalizable, eigenvecs = _is_diagonalizable_with_eigen(M,
694 reals_only=reals_only)
696 if not is_diagonalizable:
697 raise MatrixError("Matrix is not diagonalizable")
699 if sort:
700 eigenvecs = sorted(eigenvecs, key=default_sort_key)
702 p_cols, diag = [], []
704 for val, mult, basis in eigenvecs:
705 diag += [val] * mult
706 p_cols += basis
708 if normalize:
709 p_cols = [v / v.norm() for v in p_cols]
711 return M.hstack(*p_cols), M.diag(*diag)
714def _fuzzy_positive_definite(M):
715 positive_diagonals = M._has_positive_diagonals()
716 if positive_diagonals is False:
717 return False
719 if positive_diagonals and M.is_strongly_diagonally_dominant:
720 return True
722 return None
725def _fuzzy_positive_semidefinite(M):
726 nonnegative_diagonals = M._has_nonnegative_diagonals()
727 if nonnegative_diagonals is False:
728 return False
730 if nonnegative_diagonals and M.is_weakly_diagonally_dominant:
731 return True
733 return None
736def _is_positive_definite(M):
737 if not M.is_hermitian:
738 if not M.is_square:
739 return False
740 M = M + M.H
742 fuzzy = _fuzzy_positive_definite(M)
743 if fuzzy is not None:
744 return fuzzy
746 return _is_positive_definite_GE(M)
749def _is_positive_semidefinite(M):
750 if not M.is_hermitian:
751 if not M.is_square:
752 return False
753 M = M + M.H
755 fuzzy = _fuzzy_positive_semidefinite(M)
756 if fuzzy is not None:
757 return fuzzy
759 return _is_positive_semidefinite_cholesky(M)
762def _is_negative_definite(M):
763 return _is_positive_definite(-M)
766def _is_negative_semidefinite(M):
767 return _is_positive_semidefinite(-M)
770def _is_indefinite(M):
771 if M.is_hermitian:
772 eigen = M.eigenvals()
773 args1 = [x.is_positive for x in eigen.keys()]
774 any_positive = fuzzy_or(args1)
775 args2 = [x.is_negative for x in eigen.keys()]
776 any_negative = fuzzy_or(args2)
778 return fuzzy_and([any_positive, any_negative])
780 elif M.is_square:
781 return (M + M.H).is_indefinite
783 return False
786def _is_positive_definite_GE(M):
787 """A division-free gaussian elimination method for testing
788 positive-definiteness."""
789 M = M.as_mutable()
790 size = M.rows
792 for i in range(size):
793 is_positive = M[i, i].is_positive
794 if is_positive is not True:
795 return is_positive
796 for j in range(i+1, size):
797 M[j, i+1:] = M[i, i] * M[j, i+1:] - M[j, i] * M[i, i+1:]
798 return True
801def _is_positive_semidefinite_cholesky(M):
802 """Uses Cholesky factorization with complete pivoting
804 References
805 ==========
807 .. [1] http://eprints.ma.man.ac.uk/1199/1/covered/MIMS_ep2008_116.pdf
809 .. [2] https://www.value-at-risk.net/cholesky-factorization/
810 """
811 M = M.as_mutable()
812 for k in range(M.rows):
813 diags = [M[i, i] for i in range(k, M.rows)]
814 pivot, pivot_val, nonzero, _ = _find_reasonable_pivot(diags)
816 if nonzero:
817 return None
819 if pivot is None:
820 for i in range(k+1, M.rows):
821 for j in range(k, M.cols):
822 iszero = M[i, j].is_zero
823 if iszero is None:
824 return None
825 elif iszero is False:
826 return False
827 return True
829 if M[k, k].is_negative or pivot_val.is_negative:
830 return False
831 elif not (M[k, k].is_nonnegative and pivot_val.is_nonnegative):
832 return None
834 if pivot > 0:
835 M.col_swap(k, k+pivot)
836 M.row_swap(k, k+pivot)
838 M[k, k] = sqrt(M[k, k])
839 M[k, k+1:] /= M[k, k]
840 M[k+1:, k+1:] -= M[k, k+1:].H * M[k, k+1:]
842 return M[-1, -1].is_nonnegative
845_doc_positive_definite = \
846 r"""Finds out the definiteness of a matrix.
848 Explanation
849 ===========
851 A square real matrix $A$ is:
853 - A positive definite matrix if $x^T A x > 0$
854 for all non-zero real vectors $x$.
855 - A positive semidefinite matrix if $x^T A x \geq 0$
856 for all non-zero real vectors $x$.
857 - A negative definite matrix if $x^T A x < 0$
858 for all non-zero real vectors $x$.
859 - A negative semidefinite matrix if $x^T A x \leq 0$
860 for all non-zero real vectors $x$.
861 - An indefinite matrix if there exists non-zero real vectors
862 $x, y$ with $x^T A x > 0 > y^T A y$.
864 A square complex matrix $A$ is:
866 - A positive definite matrix if $\text{re}(x^H A x) > 0$
867 for all non-zero complex vectors $x$.
868 - A positive semidefinite matrix if $\text{re}(x^H A x) \geq 0$
869 for all non-zero complex vectors $x$.
870 - A negative definite matrix if $\text{re}(x^H A x) < 0$
871 for all non-zero complex vectors $x$.
872 - A negative semidefinite matrix if $\text{re}(x^H A x) \leq 0$
873 for all non-zero complex vectors $x$.
874 - An indefinite matrix if there exists non-zero complex vectors
875 $x, y$ with $\text{re}(x^H A x) > 0 > \text{re}(y^H A y)$.
877 A matrix need not be symmetric or hermitian to be positive definite.
879 - A real non-symmetric matrix is positive definite if and only if
880 $\frac{A + A^T}{2}$ is positive definite.
881 - A complex non-hermitian matrix is positive definite if and only if
882 $\frac{A + A^H}{2}$ is positive definite.
884 And this extension can apply for all the definitions above.
886 However, for complex cases, you can restrict the definition of
887 $\text{re}(x^H A x) > 0$ to $x^H A x > 0$ and require the matrix
888 to be hermitian.
889 But we do not present this restriction for computation because you
890 can check ``M.is_hermitian`` independently with this and use
891 the same procedure.
893 Examples
894 ========
896 An example of symmetric positive definite matrix:
898 .. plot::
899 :context: reset
900 :format: doctest
901 :include-source: True
903 >>> from sympy import Matrix, symbols
904 >>> from sympy.plotting import plot3d
905 >>> a, b = symbols('a b')
906 >>> x = Matrix([a, b])
908 >>> A = Matrix([[1, 0], [0, 1]])
909 >>> A.is_positive_definite
910 True
911 >>> A.is_positive_semidefinite
912 True
914 >>> p = plot3d((x.T*A*x)[0, 0], (a, -1, 1), (b, -1, 1))
916 An example of symmetric positive semidefinite matrix:
918 .. plot::
919 :context: close-figs
920 :format: doctest
921 :include-source: True
923 >>> A = Matrix([[1, -1], [-1, 1]])
924 >>> A.is_positive_definite
925 False
926 >>> A.is_positive_semidefinite
927 True
929 >>> p = plot3d((x.T*A*x)[0, 0], (a, -1, 1), (b, -1, 1))
931 An example of symmetric negative definite matrix:
933 .. plot::
934 :context: close-figs
935 :format: doctest
936 :include-source: True
938 >>> A = Matrix([[-1, 0], [0, -1]])
939 >>> A.is_negative_definite
940 True
941 >>> A.is_negative_semidefinite
942 True
943 >>> A.is_indefinite
944 False
946 >>> p = plot3d((x.T*A*x)[0, 0], (a, -1, 1), (b, -1, 1))
948 An example of symmetric indefinite matrix:
950 .. plot::
951 :context: close-figs
952 :format: doctest
953 :include-source: True
955 >>> A = Matrix([[1, 2], [2, -1]])
956 >>> A.is_indefinite
957 True
959 >>> p = plot3d((x.T*A*x)[0, 0], (a, -1, 1), (b, -1, 1))
961 An example of non-symmetric positive definite matrix.
963 .. plot::
964 :context: close-figs
965 :format: doctest
966 :include-source: True
968 >>> A = Matrix([[1, 2], [-2, 1]])
969 >>> A.is_positive_definite
970 True
971 >>> A.is_positive_semidefinite
972 True
974 >>> p = plot3d((x.T*A*x)[0, 0], (a, -1, 1), (b, -1, 1))
976 Notes
977 =====
979 Although some people trivialize the definition of positive definite
980 matrices only for symmetric or hermitian matrices, this restriction
981 is not correct because it does not classify all instances of
982 positive definite matrices from the definition $x^T A x > 0$ or
983 $\text{re}(x^H A x) > 0$.
985 For instance, ``Matrix([[1, 2], [-2, 1]])`` presented in
986 the example above is an example of real positive definite matrix
987 that is not symmetric.
989 However, since the following formula holds true;
991 .. math::
992 \text{re}(x^H A x) > 0 \iff
993 \text{re}(x^H \frac{A + A^H}{2} x) > 0
995 We can classify all positive definite matrices that may or may not
996 be symmetric or hermitian by transforming the matrix to
997 $\frac{A + A^T}{2}$ or $\frac{A + A^H}{2}$
998 (which is guaranteed to be always real symmetric or complex
999 hermitian) and we can defer most of the studies to symmetric or
1000 hermitian positive definite matrices.
1002 But it is a different problem for the existance of Cholesky
1003 decomposition. Because even though a non symmetric or a non
1004 hermitian matrix can be positive definite, Cholesky or LDL
1005 decomposition does not exist because the decompositions require the
1006 matrix to be symmetric or hermitian.
1008 References
1009 ==========
1011 .. [1] https://en.wikipedia.org/wiki/Definiteness_of_a_matrix#Eigenvalues
1013 .. [2] https://mathworld.wolfram.com/PositiveDefiniteMatrix.html
1015 .. [3] Johnson, C. R. "Positive Definite Matrices." Amer.
1016 Math. Monthly 77, 259-264 1970.
1017 """
1019_is_positive_definite.__doc__ = _doc_positive_definite
1020_is_positive_semidefinite.__doc__ = _doc_positive_definite
1021_is_negative_definite.__doc__ = _doc_positive_definite
1022_is_negative_semidefinite.__doc__ = _doc_positive_definite
1023_is_indefinite.__doc__ = _doc_positive_definite
1026def _jordan_form(M, calc_transform=True, *, chop=False):
1027 """Return $(P, J)$ where $J$ is a Jordan block
1028 matrix and $P$ is a matrix such that $M = P J P^{-1}$
1030 Parameters
1031 ==========
1033 calc_transform : bool
1034 If ``False``, then only $J$ is returned.
1036 chop : bool
1037 All matrices are converted to exact types when computing
1038 eigenvalues and eigenvectors. As a result, there may be
1039 approximation errors. If ``chop==True``, these errors
1040 will be truncated.
1042 Examples
1043 ========
1045 >>> from sympy import Matrix
1046 >>> M = Matrix([[ 6, 5, -2, -3], [-3, -1, 3, 3], [ 2, 1, -2, -3], [-1, 1, 5, 5]])
1047 >>> P, J = M.jordan_form()
1048 >>> J
1049 Matrix([
1050 [2, 1, 0, 0],
1051 [0, 2, 0, 0],
1052 [0, 0, 2, 1],
1053 [0, 0, 0, 2]])
1055 See Also
1056 ========
1058 jordan_block
1059 """
1061 if not M.is_square:
1062 raise NonSquareMatrixError("Only square matrices have Jordan forms")
1064 mat = M
1065 has_floats = M.has(Float)
1067 if has_floats:
1068 try:
1069 max_prec = max(term._prec for term in M.values() if isinstance(term, Float))
1070 except ValueError:
1071 # if no term in the matrix is explicitly a Float calling max()
1072 # will throw a error so setting max_prec to default value of 53
1073 max_prec = 53
1075 # setting minimum max_dps to 15 to prevent loss of precision in
1076 # matrix containing non evaluated expressions
1077 max_dps = max(prec_to_dps(max_prec), 15)
1079 def restore_floats(*args):
1080 """If ``has_floats`` is `True`, cast all ``args`` as
1081 matrices of floats."""
1083 if has_floats:
1084 args = [m.evalf(n=max_dps, chop=chop) for m in args]
1085 if len(args) == 1:
1086 return args[0]
1088 return args
1090 # cache calculations for some speedup
1091 mat_cache = {}
1093 def eig_mat(val, pow):
1094 """Cache computations of ``(M - val*I)**pow`` for quick
1095 retrieval"""
1097 if (val, pow) in mat_cache:
1098 return mat_cache[(val, pow)]
1100 if (val, pow - 1) in mat_cache:
1101 mat_cache[(val, pow)] = mat_cache[(val, pow - 1)].multiply(
1102 mat_cache[(val, 1)], dotprodsimp=None)
1103 else:
1104 mat_cache[(val, pow)] = (mat - val*M.eye(M.rows)).pow(pow)
1106 return mat_cache[(val, pow)]
1108 # helper functions
1109 def nullity_chain(val, algebraic_multiplicity):
1110 """Calculate the sequence [0, nullity(E), nullity(E**2), ...]
1111 until it is constant where ``E = M - val*I``"""
1113 # mat.rank() is faster than computing the null space,
1114 # so use the rank-nullity theorem
1115 cols = M.cols
1116 ret = [0]
1117 nullity = cols - eig_mat(val, 1).rank()
1118 i = 2
1120 while nullity != ret[-1]:
1121 ret.append(nullity)
1123 if nullity == algebraic_multiplicity:
1124 break
1126 nullity = cols - eig_mat(val, i).rank()
1127 i += 1
1129 # Due to issues like #7146 and #15872, SymPy sometimes
1130 # gives the wrong rank. In this case, raise an error
1131 # instead of returning an incorrect matrix
1132 if nullity < ret[-1] or nullity > algebraic_multiplicity:
1133 raise MatrixError(
1134 "SymPy had encountered an inconsistent "
1135 "result while computing Jordan block: "
1136 "{}".format(M))
1138 return ret
1140 def blocks_from_nullity_chain(d):
1141 """Return a list of the size of each Jordan block.
1142 If d_n is the nullity of E**n, then the number
1143 of Jordan blocks of size n is
1145 2*d_n - d_(n-1) - d_(n+1)"""
1147 # d[0] is always the number of columns, so skip past it
1148 mid = [2*d[n] - d[n - 1] - d[n + 1] for n in range(1, len(d) - 1)]
1149 # d is assumed to plateau with "d[ len(d) ] == d[-1]", so
1150 # 2*d_n - d_(n-1) - d_(n+1) == d_n - d_(n-1)
1151 end = [d[-1] - d[-2]] if len(d) > 1 else [d[0]]
1153 return mid + end
1155 def pick_vec(small_basis, big_basis):
1156 """Picks a vector from big_basis that isn't in
1157 the subspace spanned by small_basis"""
1159 if len(small_basis) == 0:
1160 return big_basis[0]
1162 for v in big_basis:
1163 _, pivots = M.hstack(*(small_basis + [v])).echelon_form(
1164 with_pivots=True)
1166 if pivots[-1] == len(small_basis):
1167 return v
1169 # roots doesn't like Floats, so replace them with Rationals
1170 if has_floats:
1171 from sympy.simplify import nsimplify
1172 mat = mat.applyfunc(lambda x: nsimplify(x, rational=True))
1174 # first calculate the jordan block structure
1175 eigs = mat.eigenvals()
1177 # Make sure that we have all roots in radical form
1178 for x in eigs:
1179 if x.has(CRootOf):
1180 raise MatrixError(
1181 "Jordan normal form is not implemented if the matrix have "
1182 "eigenvalues in CRootOf form")
1184 # most matrices have distinct eigenvalues
1185 # and so are diagonalizable. In this case, don't
1186 # do extra work!
1187 if len(eigs.keys()) == mat.cols:
1188 blocks = sorted(eigs.keys(), key=default_sort_key)
1189 jordan_mat = mat.diag(*blocks)
1191 if not calc_transform:
1192 return restore_floats(jordan_mat)
1194 jordan_basis = [eig_mat(eig, 1).nullspace()[0]
1195 for eig in blocks]
1196 basis_mat = mat.hstack(*jordan_basis)
1198 return restore_floats(basis_mat, jordan_mat)
1200 block_structure = []
1202 for eig in sorted(eigs.keys(), key=default_sort_key):
1203 algebraic_multiplicity = eigs[eig]
1204 chain = nullity_chain(eig, algebraic_multiplicity)
1205 block_sizes = blocks_from_nullity_chain(chain)
1207 # if block_sizes = = [a, b, c, ...], then the number of
1208 # Jordan blocks of size 1 is a, of size 2 is b, etc.
1209 # create an array that has (eig, block_size) with one
1210 # entry for each block
1211 size_nums = [(i+1, num) for i, num in enumerate(block_sizes)]
1213 # we expect larger Jordan blocks to come earlier
1214 size_nums.reverse()
1216 block_structure.extend(
1217 [(eig, size) for size, num in size_nums for _ in range(num)])
1219 jordan_form_size = sum(size for eig, size in block_structure)
1221 if jordan_form_size != M.rows:
1222 raise MatrixError(
1223 "SymPy had encountered an inconsistent result while "
1224 "computing Jordan block. : {}".format(M))
1226 blocks = (mat.jordan_block(size=size, eigenvalue=eig) for eig, size in block_structure)
1227 jordan_mat = mat.diag(*blocks)
1229 if not calc_transform:
1230 return restore_floats(jordan_mat)
1232 # For each generalized eigenspace, calculate a basis.
1233 # We start by looking for a vector in null( (A - eig*I)**n )
1234 # which isn't in null( (A - eig*I)**(n-1) ) where n is
1235 # the size of the Jordan block
1236 #
1237 # Ideally we'd just loop through block_structure and
1238 # compute each generalized eigenspace. However, this
1239 # causes a lot of unneeded computation. Instead, we
1240 # go through the eigenvalues separately, since we know
1241 # their generalized eigenspaces must have bases that
1242 # are linearly independent.
1243 jordan_basis = []
1245 for eig in sorted(eigs.keys(), key=default_sort_key):
1246 eig_basis = []
1248 for block_eig, size in block_structure:
1249 if block_eig != eig:
1250 continue
1252 null_big = (eig_mat(eig, size)).nullspace()
1253 null_small = (eig_mat(eig, size - 1)).nullspace()
1255 # we want to pick something that is in the big basis
1256 # and not the small, but also something that is independent
1257 # of any other generalized eigenvectors from a different
1258 # generalized eigenspace sharing the same eigenvalue.
1259 vec = pick_vec(null_small + eig_basis, null_big)
1260 new_vecs = [eig_mat(eig, i).multiply(vec, dotprodsimp=None)
1261 for i in range(size)]
1263 eig_basis.extend(new_vecs)
1264 jordan_basis.extend(reversed(new_vecs))
1266 basis_mat = mat.hstack(*jordan_basis)
1268 return restore_floats(basis_mat, jordan_mat)
1271def _left_eigenvects(M, **flags):
1272 """Returns left eigenvectors and eigenvalues.
1274 This function returns the list of triples (eigenval, multiplicity,
1275 basis) for the left eigenvectors. Options are the same as for
1276 eigenvects(), i.e. the ``**flags`` arguments gets passed directly to
1277 eigenvects().
1279 Examples
1280 ========
1282 >>> from sympy import Matrix
1283 >>> M = Matrix([[0, 1, 1], [1, 0, 0], [1, 1, 1]])
1284 >>> M.eigenvects()
1285 [(-1, 1, [Matrix([
1286 [-1],
1287 [ 1],
1288 [ 0]])]), (0, 1, [Matrix([
1289 [ 0],
1290 [-1],
1291 [ 1]])]), (2, 1, [Matrix([
1292 [2/3],
1293 [1/3],
1294 [ 1]])])]
1295 >>> M.left_eigenvects()
1296 [(-1, 1, [Matrix([[-2, 1, 1]])]), (0, 1, [Matrix([[-1, -1, 1]])]), (2,
1297 1, [Matrix([[1, 1, 1]])])]
1299 """
1301 eigs = M.transpose().eigenvects(**flags)
1303 return [(val, mult, [l.transpose() for l in basis]) for val, mult, basis in eigs]
1306def _singular_values(M):
1307 """Compute the singular values of a Matrix
1309 Examples
1310 ========
1312 >>> from sympy import Matrix, Symbol
1313 >>> x = Symbol('x', real=True)
1314 >>> M = Matrix([[0, 1, 0], [0, x, 0], [-1, 0, 0]])
1315 >>> M.singular_values()
1316 [sqrt(x**2 + 1), 1, 0]
1318 See Also
1319 ========
1321 condition_number
1322 """
1324 if M.rows >= M.cols:
1325 valmultpairs = M.H.multiply(M).eigenvals()
1326 else:
1327 valmultpairs = M.multiply(M.H).eigenvals()
1329 # Expands result from eigenvals into a simple list
1330 vals = []
1332 for k, v in valmultpairs.items():
1333 vals += [sqrt(k)] * v # dangerous! same k in several spots!
1335 # Pad with zeros if singular values are computed in reverse way,
1336 # to give consistent format.
1337 if len(vals) < M.cols:
1338 vals += [M.zero] * (M.cols - len(vals))
1340 # sort them in descending order
1341 vals.sort(reverse=True, key=default_sort_key)
1343 return vals