test_eigen.py 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712
  1. from sympy.core.evalf import N
  2. from sympy.core.numbers import (Float, I, Rational)
  3. from sympy.core.symbol import (Symbol, symbols)
  4. from sympy.functions.elementary.complexes import Abs
  5. from sympy.functions.elementary.miscellaneous import sqrt
  6. from sympy.functions.elementary.trigonometric import (cos, sin)
  7. from sympy.matrices import eye, Matrix
  8. from sympy.core.singleton import S
  9. from sympy.testing.pytest import raises, XFAIL
  10. from sympy.matrices.exceptions import NonSquareMatrixError, MatrixError
  11. from sympy.matrices.expressions.fourier import DFT
  12. from sympy.simplify.simplify import simplify
  13. from sympy.matrices.immutable import ImmutableMatrix
  14. from sympy.testing.pytest import slow
  15. from sympy.testing.matrices import allclose
  16. def test_eigen():
  17. R = Rational
  18. M = Matrix.eye(3)
  19. assert M.eigenvals(multiple=False) == {S.One: 3}
  20. assert M.eigenvals(multiple=True) == [1, 1, 1]
  21. assert M.eigenvects() == (
  22. [(1, 3, [Matrix([1, 0, 0]),
  23. Matrix([0, 1, 0]),
  24. Matrix([0, 0, 1])])])
  25. assert M.left_eigenvects() == (
  26. [(1, 3, [Matrix([[1, 0, 0]]),
  27. Matrix([[0, 1, 0]]),
  28. Matrix([[0, 0, 1]])])])
  29. M = Matrix([[0, 1, 1],
  30. [1, 0, 0],
  31. [1, 1, 1]])
  32. assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}
  33. assert M.eigenvects() == (
  34. [
  35. (-1, 1, [Matrix([-1, 1, 0])]),
  36. ( 0, 1, [Matrix([0, -1, 1])]),
  37. ( 2, 1, [Matrix([R(2, 3), R(1, 3), 1])])
  38. ])
  39. assert M.left_eigenvects() == (
  40. [
  41. (-1, 1, [Matrix([[-2, 1, 1]])]),
  42. (0, 1, [Matrix([[-1, -1, 1]])]),
  43. (2, 1, [Matrix([[1, 1, 1]])])
  44. ])
  45. a = Symbol('a')
  46. M = Matrix([[a, 0],
  47. [0, 1]])
  48. assert M.eigenvals() == {a: 1, S.One: 1}
  49. M = Matrix([[1, -1],
  50. [1, 3]])
  51. assert M.eigenvects() == ([(2, 2, [Matrix(2, 1, [-1, 1])])])
  52. assert M.left_eigenvects() == ([(2, 2, [Matrix([[1, 1]])])])
  53. M = Matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
  54. a = R(15, 2)
  55. b = 3*33**R(1, 2)
  56. c = R(13, 2)
  57. d = (R(33, 8) + 3*b/8)
  58. e = (R(33, 8) - 3*b/8)
  59. def NS(e, n):
  60. return str(N(e, n))
  61. r = [
  62. (a - b/2, 1, [Matrix([(12 + 24/(c - b/2))/((c - b/2)*e) + 3/(c - b/2),
  63. (6 + 12/(c - b/2))/e, 1])]),
  64. ( 0, 1, [Matrix([1, -2, 1])]),
  65. (a + b/2, 1, [Matrix([(12 + 24/(c + b/2))/((c + b/2)*d) + 3/(c + b/2),
  66. (6 + 12/(c + b/2))/d, 1])]),
  67. ]
  68. r1 = [(NS(r[i][0], 2), NS(r[i][1], 2),
  69. [NS(j, 2) for j in r[i][2][0]]) for i in range(len(r))]
  70. r = M.eigenvects()
  71. r2 = [(NS(r[i][0], 2), NS(r[i][1], 2),
  72. [NS(j, 2) for j in r[i][2][0]]) for i in range(len(r))]
  73. assert sorted(r1) == sorted(r2)
  74. eps = Symbol('eps', real=True)
  75. M = Matrix([[abs(eps), I*eps ],
  76. [-I*eps, abs(eps) ]])
  77. assert M.eigenvects() == (
  78. [
  79. ( 0, 1, [Matrix([[-I*eps/abs(eps)], [1]])]),
  80. ( 2*abs(eps), 1, [ Matrix([[I*eps/abs(eps)], [1]]) ] ),
  81. ])
  82. assert M.left_eigenvects() == (
  83. [
  84. (0, 1, [Matrix([[I*eps/Abs(eps), 1]])]),
  85. (2*Abs(eps), 1, [Matrix([[-I*eps/Abs(eps), 1]])])
  86. ])
  87. M = Matrix(3, 3, [1, 2, 0, 0, 3, 0, 2, -4, 2])
  88. M._eigenvects = M.eigenvects(simplify=False)
  89. assert max(i.q for i in M._eigenvects[0][2][0]) > 1
  90. M._eigenvects = M.eigenvects(simplify=True)
  91. assert max(i.q for i in M._eigenvects[0][2][0]) == 1
  92. M = Matrix([[Rational(1, 4), 1], [1, 1]])
  93. assert M.eigenvects() == [
  94. (Rational(5, 8) - sqrt(73)/8, 1, [Matrix([[-sqrt(73)/8 - Rational(3, 8)], [1]])]),
  95. (Rational(5, 8) + sqrt(73)/8, 1, [Matrix([[Rational(-3, 8) + sqrt(73)/8], [1]])])]
  96. # issue 10719
  97. assert Matrix([]).eigenvals() == {}
  98. assert Matrix([]).eigenvals(multiple=True) == []
  99. assert Matrix([]).eigenvects() == []
  100. # issue 15119
  101. raises(NonSquareMatrixError,
  102. lambda: Matrix([[1, 2], [0, 4], [0, 0]]).eigenvals())
  103. raises(NonSquareMatrixError,
  104. lambda: Matrix([[1, 0], [3, 4], [5, 6]]).eigenvals())
  105. raises(NonSquareMatrixError,
  106. lambda: Matrix([[1, 2, 3], [0, 5, 6]]).eigenvals())
  107. raises(NonSquareMatrixError,
  108. lambda: Matrix([[1, 0, 0], [4, 5, 0]]).eigenvals())
  109. raises(NonSquareMatrixError,
  110. lambda: Matrix([[1, 2, 3], [0, 5, 6]]).eigenvals(
  111. error_when_incomplete = False))
  112. raises(NonSquareMatrixError,
  113. lambda: Matrix([[1, 0, 0], [4, 5, 0]]).eigenvals(
  114. error_when_incomplete = False))
  115. m = Matrix([[1, 2], [3, 4]])
  116. assert isinstance(m.eigenvals(simplify=True, multiple=False), dict)
  117. assert isinstance(m.eigenvals(simplify=True, multiple=True), list)
  118. assert isinstance(m.eigenvals(simplify=lambda x: x, multiple=False), dict)
  119. assert isinstance(m.eigenvals(simplify=lambda x: x, multiple=True), list)
  120. def test_float_eigenvals():
  121. m = Matrix([[1, .6, .6], [.6, .9, .9], [.9, .6, .6]])
  122. evals = [
  123. Rational(5, 4) - sqrt(385)/20,
  124. sqrt(385)/20 + Rational(5, 4),
  125. S.Zero]
  126. n_evals = m.eigenvals(rational=True, multiple=True)
  127. n_evals = sorted(n_evals)
  128. s_evals = [x.evalf() for x in evals]
  129. s_evals = sorted(s_evals)
  130. for x, y in zip(n_evals, s_evals):
  131. assert abs(x-y) < 10**-9
  132. @XFAIL
  133. def test_eigen_vects():
  134. m = Matrix(2, 2, [1, 0, 0, I])
  135. raises(NotImplementedError, lambda: m.is_diagonalizable(True))
  136. # !!! bug because of eigenvects() or roots(x**2 + (-1 - I)*x + I, x)
  137. # see issue 5292
  138. assert not m.is_diagonalizable(True)
  139. raises(MatrixError, lambda: m.diagonalize(True))
  140. (P, D) = m.diagonalize(True)
  141. def test_issue_8240():
  142. # Eigenvalues of large triangular matrices
  143. x, y = symbols('x y')
  144. n = 200
  145. diagonal_variables = [Symbol('x%s' % i) for i in range(n)]
  146. M = [[0 for i in range(n)] for j in range(n)]
  147. for i in range(n):
  148. M[i][i] = diagonal_variables[i]
  149. M = Matrix(M)
  150. eigenvals = M.eigenvals()
  151. assert len(eigenvals) == n
  152. for i in range(n):
  153. assert eigenvals[diagonal_variables[i]] == 1
  154. eigenvals = M.eigenvals(multiple=True)
  155. assert set(eigenvals) == set(diagonal_variables)
  156. # with multiplicity
  157. M = Matrix([[x, 0, 0], [1, y, 0], [2, 3, x]])
  158. eigenvals = M.eigenvals()
  159. assert eigenvals == {x: 2, y: 1}
  160. eigenvals = M.eigenvals(multiple=True)
  161. assert len(eigenvals) == 3
  162. assert eigenvals.count(x) == 2
  163. assert eigenvals.count(y) == 1
  164. def test_eigenvals():
  165. M = Matrix([[0, 1, 1],
  166. [1, 0, 0],
  167. [1, 1, 1]])
  168. assert M.eigenvals() == {2*S.One: 1, -S.One: 1, S.Zero: 1}
  169. m = Matrix([
  170. [3, 0, 0, 0, -3],
  171. [0, -3, -3, 0, 3],
  172. [0, 3, 0, 3, 0],
  173. [0, 0, 3, 0, 3],
  174. [3, 0, 0, 3, 0]])
  175. # XXX Used dry-run test because arbitrary symbol that appears in
  176. # CRootOf may not be unique.
  177. assert m.eigenvals()
  178. def test_eigenvects():
  179. M = Matrix([[0, 1, 1],
  180. [1, 0, 0],
  181. [1, 1, 1]])
  182. vecs = M.eigenvects()
  183. for val, mult, vec_list in vecs:
  184. assert len(vec_list) == 1
  185. assert M*vec_list[0] == val*vec_list[0]
  186. def test_left_eigenvects():
  187. M = Matrix([[0, 1, 1],
  188. [1, 0, 0],
  189. [1, 1, 1]])
  190. vecs = M.left_eigenvects()
  191. for val, mult, vec_list in vecs:
  192. assert len(vec_list) == 1
  193. assert vec_list[0]*M == val*vec_list[0]
  194. @slow
  195. def test_bidiagonalize():
  196. M = Matrix([[1, 0, 0],
  197. [0, 1, 0],
  198. [0, 0, 1]])
  199. assert M.bidiagonalize() == M
  200. assert M.bidiagonalize(upper=False) == M
  201. assert M.bidiagonalize() == M
  202. assert M.bidiagonal_decomposition() == (M, M, M)
  203. assert M.bidiagonal_decomposition(upper=False) == (M, M, M)
  204. assert M.bidiagonalize() == M
  205. import random
  206. #Real Tests
  207. for real_test in range(2):
  208. test_values = []
  209. row = 2
  210. col = 2
  211. for _ in range(row * col):
  212. value = random.randint(-1000000000, 1000000000)
  213. test_values = test_values + [value]
  214. # L -> Lower Bidiagonalization
  215. # M -> Mutable Matrix
  216. # N -> Immutable Matrix
  217. # 0 -> Bidiagonalized form
  218. # 1,2,3 -> Bidiagonal_decomposition matrices
  219. # 4 -> Product of 1 2 3
  220. M = Matrix(row, col, test_values)
  221. N = ImmutableMatrix(M)
  222. N1, N2, N3 = N.bidiagonal_decomposition()
  223. M1, M2, M3 = M.bidiagonal_decomposition()
  224. M0 = M.bidiagonalize()
  225. N0 = N.bidiagonalize()
  226. N4 = N1 * N2 * N3
  227. M4 = M1 * M2 * M3
  228. N2.simplify()
  229. N4.simplify()
  230. N0.simplify()
  231. M0.simplify()
  232. M2.simplify()
  233. M4.simplify()
  234. LM0 = M.bidiagonalize(upper=False)
  235. LM1, LM2, LM3 = M.bidiagonal_decomposition(upper=False)
  236. LN0 = N.bidiagonalize(upper=False)
  237. LN1, LN2, LN3 = N.bidiagonal_decomposition(upper=False)
  238. LN4 = LN1 * LN2 * LN3
  239. LM4 = LM1 * LM2 * LM3
  240. LN2.simplify()
  241. LN4.simplify()
  242. LN0.simplify()
  243. LM0.simplify()
  244. LM2.simplify()
  245. LM4.simplify()
  246. assert M == M4
  247. assert M2 == M0
  248. assert N == N4
  249. assert N2 == N0
  250. assert M == LM4
  251. assert LM2 == LM0
  252. assert N == LN4
  253. assert LN2 == LN0
  254. #Complex Tests
  255. for complex_test in range(2):
  256. test_values = []
  257. size = 2
  258. for _ in range(size * size):
  259. real = random.randint(-1000000000, 1000000000)
  260. comp = random.randint(-1000000000, 1000000000)
  261. value = real + comp * I
  262. test_values = test_values + [value]
  263. M = Matrix(size, size, test_values)
  264. N = ImmutableMatrix(M)
  265. # L -> Lower Bidiagonalization
  266. # M -> Mutable Matrix
  267. # N -> Immutable Matrix
  268. # 0 -> Bidiagonalized form
  269. # 1,2,3 -> Bidiagonal_decomposition matrices
  270. # 4 -> Product of 1 2 3
  271. N1, N2, N3 = N.bidiagonal_decomposition()
  272. M1, M2, M3 = M.bidiagonal_decomposition()
  273. M0 = M.bidiagonalize()
  274. N0 = N.bidiagonalize()
  275. N4 = N1 * N2 * N3
  276. M4 = M1 * M2 * M3
  277. N2.simplify()
  278. N4.simplify()
  279. N0.simplify()
  280. M0.simplify()
  281. M2.simplify()
  282. M4.simplify()
  283. LM0 = M.bidiagonalize(upper=False)
  284. LM1, LM2, LM3 = M.bidiagonal_decomposition(upper=False)
  285. LN0 = N.bidiagonalize(upper=False)
  286. LN1, LN2, LN3 = N.bidiagonal_decomposition(upper=False)
  287. LN4 = LN1 * LN2 * LN3
  288. LM4 = LM1 * LM2 * LM3
  289. LN2.simplify()
  290. LN4.simplify()
  291. LN0.simplify()
  292. LM0.simplify()
  293. LM2.simplify()
  294. LM4.simplify()
  295. assert M == M4
  296. assert M2 == M0
  297. assert N == N4
  298. assert N2 == N0
  299. assert M == LM4
  300. assert LM2 == LM0
  301. assert N == LN4
  302. assert LN2 == LN0
  303. M = Matrix(18, 8, range(1, 145))
  304. M = M.applyfunc(lambda i: Float(i))
  305. assert M.bidiagonal_decomposition()[1] == M.bidiagonalize()
  306. assert M.bidiagonal_decomposition(upper=False)[1] == M.bidiagonalize(upper=False)
  307. a, b, c = M.bidiagonal_decomposition()
  308. diff = a * b * c - M
  309. assert abs(max(diff)) < 10**-12
  310. def test_diagonalize():
  311. m = Matrix(2, 2, [0, -1, 1, 0])
  312. raises(MatrixError, lambda: m.diagonalize(reals_only=True))
  313. P, D = m.diagonalize()
  314. assert D.is_diagonal()
  315. assert D == Matrix([
  316. [-I, 0],
  317. [ 0, I]])
  318. # make sure we use floats out if floats are passed in
  319. m = Matrix(2, 2, [0, .5, .5, 0])
  320. P, D = m.diagonalize()
  321. assert all(isinstance(e, Float) for e in D.values())
  322. assert all(isinstance(e, Float) for e in P.values())
  323. _, D2 = m.diagonalize(reals_only=True)
  324. assert D == D2
  325. m = Matrix(
  326. [[0, 1, 0, 0], [1, 0, 0, 0.002], [0.002, 0, 0, 1], [0, 0, 1, 0]])
  327. P, D = m.diagonalize()
  328. assert allclose(P*D, m*P)
  329. def test_is_diagonalizable():
  330. a, b, c = symbols('a b c')
  331. m = Matrix(2, 2, [a, c, c, b])
  332. assert m.is_symmetric()
  333. assert m.is_diagonalizable()
  334. assert not Matrix(2, 2, [1, 1, 0, 1]).is_diagonalizable()
  335. m = Matrix(2, 2, [0, -1, 1, 0])
  336. assert m.is_diagonalizable()
  337. assert not m.is_diagonalizable(reals_only=True)
  338. def test_jordan_form():
  339. m = Matrix(3, 2, [-3, 1, -3, 20, 3, 10])
  340. raises(NonSquareMatrixError, lambda: m.jordan_form())
  341. # the next two tests test the cases where the old
  342. # algorithm failed due to the fact that the block structure can
  343. # *NOT* be determined from algebraic and geometric multiplicity alone
  344. # This can be seen most easily when one lets compute the J.c.f. of a matrix that
  345. # is in J.c.f already.
  346. m = Matrix(4, 4, [2, 1, 0, 0,
  347. 0, 2, 1, 0,
  348. 0, 0, 2, 0,
  349. 0, 0, 0, 2
  350. ])
  351. P, J = m.jordan_form()
  352. assert m == J
  353. m = Matrix(4, 4, [2, 1, 0, 0,
  354. 0, 2, 0, 0,
  355. 0, 0, 2, 1,
  356. 0, 0, 0, 2
  357. ])
  358. P, J = m.jordan_form()
  359. assert m == J
  360. A = Matrix([[ 2, 4, 1, 0],
  361. [-4, 2, 0, 1],
  362. [ 0, 0, 2, 4],
  363. [ 0, 0, -4, 2]])
  364. P, J = A.jordan_form()
  365. assert simplify(P*J*P.inv()) == A
  366. assert Matrix(1, 1, [1]).jordan_form() == (Matrix([1]), Matrix([1]))
  367. assert Matrix(1, 1, [1]).jordan_form(calc_transform=False) == Matrix([1])
  368. # If we have eigenvalues in CRootOf form, raise errors
  369. m = Matrix([[3, 0, 0, 0, -3], [0, -3, -3, 0, 3], [0, 3, 0, 3, 0], [0, 0, 3, 0, 3], [3, 0, 0, 3, 0]])
  370. raises(MatrixError, lambda: m.jordan_form())
  371. # make sure that if the input has floats, the output does too
  372. m = Matrix([
  373. [ 0.6875, 0.125 + 0.1875*sqrt(3)],
  374. [0.125 + 0.1875*sqrt(3), 0.3125]])
  375. P, J = m.jordan_form()
  376. assert all(isinstance(x, Float) or x == 0 for x in P)
  377. assert all(isinstance(x, Float) or x == 0 for x in J)
  378. def test_singular_values():
  379. x = Symbol('x', real=True)
  380. A = Matrix([[0, 1*I], [2, 0]])
  381. # if singular values can be sorted, they should be in decreasing order
  382. assert A.singular_values() == [2, 1]
  383. A = eye(3)
  384. A[1, 1] = x
  385. A[2, 2] = 5
  386. vals = A.singular_values()
  387. # since Abs(x) cannot be sorted, test set equality
  388. assert set(vals) == {5, 1, Abs(x)}
  389. A = Matrix([[sin(x), cos(x)], [-cos(x), sin(x)]])
  390. vals = [sv.trigsimp() for sv in A.singular_values()]
  391. assert vals == [S.One, S.One]
  392. A = Matrix([
  393. [2, 4],
  394. [1, 3],
  395. [0, 0],
  396. [0, 0]
  397. ])
  398. assert A.singular_values() == \
  399. [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221))]
  400. assert A.T.singular_values() == \
  401. [sqrt(sqrt(221) + 15), sqrt(15 - sqrt(221)), 0, 0]
  402. def test___eq__():
  403. assert (Matrix(
  404. [[0, 1, 1],
  405. [1, 0, 0],
  406. [1, 1, 1]]) == {}) is False
  407. def test_definite():
  408. # Examples from Gilbert Strang, "Introduction to Linear Algebra"
  409. # Positive definite matrices
  410. m = Matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
  411. assert m.is_positive_definite == True
  412. assert m.is_positive_semidefinite == True
  413. assert m.is_negative_definite == False
  414. assert m.is_negative_semidefinite == False
  415. assert m.is_indefinite == False
  416. m = Matrix([[5, 4], [4, 5]])
  417. assert m.is_positive_definite == True
  418. assert m.is_positive_semidefinite == True
  419. assert m.is_negative_definite == False
  420. assert m.is_negative_semidefinite == False
  421. assert m.is_indefinite == False
  422. # Positive semidefinite matrices
  423. m = Matrix([[2, -1, -1], [-1, 2, -1], [-1, -1, 2]])
  424. assert m.is_positive_definite == False
  425. assert m.is_positive_semidefinite == True
  426. assert m.is_negative_definite == False
  427. assert m.is_negative_semidefinite == False
  428. assert m.is_indefinite == False
  429. m = Matrix([[1, 2], [2, 4]])
  430. assert m.is_positive_definite == False
  431. assert m.is_positive_semidefinite == True
  432. assert m.is_negative_definite == False
  433. assert m.is_negative_semidefinite == False
  434. assert m.is_indefinite == False
  435. # Examples from Mathematica documentation
  436. # Non-hermitian positive definite matrices
  437. m = Matrix([[2, 3], [4, 8]])
  438. assert m.is_positive_definite == True
  439. assert m.is_positive_semidefinite == True
  440. assert m.is_negative_definite == False
  441. assert m.is_negative_semidefinite == False
  442. assert m.is_indefinite == False
  443. # Hermetian matrices
  444. m = Matrix([[1, 2*I], [-I, 4]])
  445. assert m.is_positive_definite == True
  446. assert m.is_positive_semidefinite == True
  447. assert m.is_negative_definite == False
  448. assert m.is_negative_semidefinite == False
  449. assert m.is_indefinite == False
  450. # Symbolic matrices examples
  451. a = Symbol('a', positive=True)
  452. b = Symbol('b', negative=True)
  453. m = Matrix([[a, 0, 0], [0, a, 0], [0, 0, a]])
  454. assert m.is_positive_definite == True
  455. assert m.is_positive_semidefinite == True
  456. assert m.is_negative_definite == False
  457. assert m.is_negative_semidefinite == False
  458. assert m.is_indefinite == False
  459. m = Matrix([[b, 0, 0], [0, b, 0], [0, 0, b]])
  460. assert m.is_positive_definite == False
  461. assert m.is_positive_semidefinite == False
  462. assert m.is_negative_definite == True
  463. assert m.is_negative_semidefinite == True
  464. assert m.is_indefinite == False
  465. m = Matrix([[a, 0], [0, b]])
  466. assert m.is_positive_definite == False
  467. assert m.is_positive_semidefinite == False
  468. assert m.is_negative_definite == False
  469. assert m.is_negative_semidefinite == False
  470. assert m.is_indefinite == True
  471. m = Matrix([
  472. [0.0228202735623867, 0.00518748979085398,
  473. -0.0743036351048907, -0.00709135324903921],
  474. [0.00518748979085398, 0.0349045359786350,
  475. 0.0830317991056637, 0.00233147902806909],
  476. [-0.0743036351048907, 0.0830317991056637,
  477. 1.15859676366277, 0.340359081555988],
  478. [-0.00709135324903921, 0.00233147902806909,
  479. 0.340359081555988, 0.928147644848199]
  480. ])
  481. assert m.is_positive_definite == True
  482. assert m.is_positive_semidefinite == True
  483. assert m.is_indefinite == False
  484. # test for issue 19547: https://github.com/sympy/sympy/issues/19547
  485. m = Matrix([
  486. [0, 0, 0],
  487. [0, 1, 2],
  488. [0, 2, 1]
  489. ])
  490. assert not m.is_positive_definite
  491. assert not m.is_positive_semidefinite
  492. def test_positive_semidefinite_cholesky():
  493. from sympy.matrices.eigen import _is_positive_semidefinite_cholesky
  494. m = Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
  495. assert _is_positive_semidefinite_cholesky(m) == True
  496. m = Matrix([[0, 0, 0], [0, 5, -10*I], [0, 10*I, 5]])
  497. assert _is_positive_semidefinite_cholesky(m) == False
  498. m = Matrix([[1, 0, 0], [0, 0, 0], [0, 0, -1]])
  499. assert _is_positive_semidefinite_cholesky(m) == False
  500. m = Matrix([[0, 1], [1, 0]])
  501. assert _is_positive_semidefinite_cholesky(m) == False
  502. # https://www.value-at-risk.net/cholesky-factorization/
  503. m = Matrix([[4, -2, -6], [-2, 10, 9], [-6, 9, 14]])
  504. assert _is_positive_semidefinite_cholesky(m) == True
  505. m = Matrix([[9, -3, 3], [-3, 2, 1], [3, 1, 6]])
  506. assert _is_positive_semidefinite_cholesky(m) == True
  507. m = Matrix([[4, -2, 2], [-2, 1, -1], [2, -1, 5]])
  508. assert _is_positive_semidefinite_cholesky(m) == True
  509. m = Matrix([[1, 2, -1], [2, 5, 1], [-1, 1, 9]])
  510. assert _is_positive_semidefinite_cholesky(m) == False
  511. def test_issue_20582():
  512. A = Matrix([
  513. [5, -5, -3, 2, -7],
  514. [-2, -5, 0, 2, 1],
  515. [-2, -7, -5, -2, -6],
  516. [7, 10, 3, 9, -2],
  517. [4, -10, 3, -8, -4]
  518. ])
  519. # XXX Used dry-run test because arbitrary symbol that appears in
  520. # CRootOf may not be unique.
  521. assert A.eigenvects()
  522. def test_issue_19210():
  523. t = Symbol('t')
  524. H = Matrix([[3, 0, 0, 0], [0, 1 , 2, 0], [0, 2, 2, 0], [0, 0, 0, 4]])
  525. A = (-I * H * t).jordan_form()
  526. assert A == (Matrix([
  527. [0, 1, 0, 0],
  528. [0, 0, -4/(-1 + sqrt(17)), 4/(1 + sqrt(17))],
  529. [0, 0, 1, 1],
  530. [1, 0, 0, 0]]), Matrix([
  531. [-4*I*t, 0, 0, 0],
  532. [ 0, -3*I*t, 0, 0],
  533. [ 0, 0, t*(-3*I/2 + sqrt(17)*I/2), 0],
  534. [ 0, 0, 0, t*(-sqrt(17)*I/2 - 3*I/2)]]))
  535. def test_issue_20275():
  536. # XXX We use complex expansions because complex exponentials are not
  537. # recognized by polys.domains
  538. A = DFT(3).as_explicit().expand(complex=True)
  539. eigenvects = A.eigenvects()
  540. assert eigenvects[0] == (
  541. -1, 1,
  542. [Matrix([[1 - sqrt(3)], [1], [1]])]
  543. )
  544. assert eigenvects[1] == (
  545. 1, 1,
  546. [Matrix([[1 + sqrt(3)], [1], [1]])]
  547. )
  548. assert eigenvects[2] == (
  549. -I, 1,
  550. [Matrix([[0], [-1], [1]])]
  551. )
  552. A = DFT(4).as_explicit().expand(complex=True)
  553. eigenvects = A.eigenvects()
  554. assert eigenvects[0] == (
  555. -1, 1,
  556. [Matrix([[-1], [1], [1], [1]])]
  557. )
  558. assert eigenvects[1] == (
  559. 1, 2,
  560. [Matrix([[1], [0], [1], [0]]), Matrix([[2], [1], [0], [1]])]
  561. )
  562. assert eigenvects[2] == (
  563. -I, 1,
  564. [Matrix([[0], [-1], [0], [1]])]
  565. )
  566. # XXX We skip test for some parts of eigenvectors which are very
  567. # complicated and fragile under expression tree changes
  568. A = DFT(5).as_explicit().expand(complex=True)
  569. eigenvects = A.eigenvects()
  570. assert eigenvects[0] == (
  571. -1, 1,
  572. [Matrix([[1 - sqrt(5)], [1], [1], [1], [1]])]
  573. )
  574. assert eigenvects[1] == (
  575. 1, 2,
  576. [Matrix([[S(1)/2 + sqrt(5)/2], [0], [1], [1], [0]]),
  577. Matrix([[S(1)/2 + sqrt(5)/2], [1], [0], [0], [1]])]
  578. )
  579. def test_issue_20752():
  580. b = symbols('b', nonzero=True)
  581. m = Matrix([[0, 0, 0], [0, b, 0], [0, 0, b]])
  582. assert m.is_positive_semidefinite is None
  583. def test_issue_25282():
  584. dd = sd = [0] * 11 + [1]
  585. ds = [2, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0]
  586. ss = ds.copy()
  587. ss[8] = 2
  588. def rotate(x, i):
  589. return x[i:] + x[:i]
  590. mat = []
  591. for i in range(12):
  592. mat.append(rotate(ss, i) + rotate(sd, i))
  593. for i in range(12):
  594. mat.append(rotate(ds, i) + rotate(dd, i))
  595. assert sum(Matrix(mat).eigenvals().values()) == 24