test_inverse.py 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193
  1. from sympy import ZZ, Matrix
  2. from sympy.polys.matrices import DM, DomainMatrix
  3. from sympy.polys.matrices.dense import ddm_iinv
  4. from sympy.polys.matrices.exceptions import DMNonInvertibleMatrixError
  5. from sympy.matrices.exceptions import NonInvertibleMatrixError
  6. import pytest
  7. from sympy.testing.pytest import raises
  8. from sympy.core.numbers import all_close
  9. from sympy.abc import x
  10. # Examples are given as adjugate matrix and determinant adj_det should match
  11. # these exactly but inv_den only matches after cancel_denom.
  12. INVERSE_EXAMPLES = [
  13. (
  14. 'zz_1',
  15. DomainMatrix([], (0, 0), ZZ),
  16. DomainMatrix([], (0, 0), ZZ),
  17. ZZ(1),
  18. ),
  19. (
  20. 'zz_2',
  21. DM([[2]], ZZ),
  22. DM([[1]], ZZ),
  23. ZZ(2),
  24. ),
  25. (
  26. 'zz_3',
  27. DM([[2, 0],
  28. [0, 2]], ZZ),
  29. DM([[2, 0],
  30. [0, 2]], ZZ),
  31. ZZ(4),
  32. ),
  33. (
  34. 'zz_4',
  35. DM([[1, 2],
  36. [3, 4]], ZZ),
  37. DM([[ 4, -2],
  38. [-3, 1]], ZZ),
  39. ZZ(-2),
  40. ),
  41. (
  42. 'zz_5',
  43. DM([[2, 2, 0],
  44. [0, 2, 2],
  45. [0, 0, 2]], ZZ),
  46. DM([[4, -4, 4],
  47. [0, 4, -4],
  48. [0, 0, 4]], ZZ),
  49. ZZ(8),
  50. ),
  51. (
  52. 'zz_6',
  53. DM([[1, 2, 3],
  54. [4, 5, 6],
  55. [7, 8, 9]], ZZ),
  56. DM([[-3, 6, -3],
  57. [ 6, -12, 6],
  58. [-3, 6, -3]], ZZ),
  59. ZZ(0),
  60. ),
  61. ]
  62. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  63. def test_Matrix_inv(name, A, A_inv, den):
  64. def _check(**kwargs):
  65. if den != 0:
  66. assert A.inv(**kwargs) == A_inv
  67. else:
  68. raises(NonInvertibleMatrixError, lambda: A.inv(**kwargs))
  69. K = A.domain
  70. A = A.to_Matrix()
  71. A_inv = A_inv.to_Matrix() / K.to_sympy(den)
  72. _check()
  73. for method in ['GE', 'LU', 'ADJ', 'CH', 'LDL', 'QR']:
  74. _check(method=method)
  75. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  76. def test_dm_inv_den(name, A, A_inv, den):
  77. if den != 0:
  78. A_inv_f, den_f = A.inv_den()
  79. assert A_inv_f.cancel_denom(den_f) == A_inv.cancel_denom(den)
  80. else:
  81. raises(DMNonInvertibleMatrixError, lambda: A.inv_den())
  82. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  83. def test_dm_inv(name, A, A_inv, den):
  84. A = A.to_field()
  85. if den != 0:
  86. A_inv = A_inv.to_field() / den
  87. assert A.inv() == A_inv
  88. else:
  89. raises(DMNonInvertibleMatrixError, lambda: A.inv())
  90. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  91. def test_ddm_inv(name, A, A_inv, den):
  92. A = A.to_field().to_ddm()
  93. if den != 0:
  94. A_inv = (A_inv.to_field() / den).to_ddm()
  95. assert A.inv() == A_inv
  96. else:
  97. raises(DMNonInvertibleMatrixError, lambda: A.inv())
  98. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  99. def test_sdm_inv(name, A, A_inv, den):
  100. A = A.to_field().to_sdm()
  101. if den != 0:
  102. A_inv = (A_inv.to_field() / den).to_sdm()
  103. assert A.inv() == A_inv
  104. else:
  105. raises(DMNonInvertibleMatrixError, lambda: A.inv())
  106. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  107. def test_dense_ddm_iinv(name, A, A_inv, den):
  108. A = A.to_field().to_ddm().copy()
  109. K = A.domain
  110. A_result = A.copy()
  111. if den != 0:
  112. A_inv = (A_inv.to_field() / den).to_ddm()
  113. ddm_iinv(A_result, A, K)
  114. assert A_result == A_inv
  115. else:
  116. raises(DMNonInvertibleMatrixError, lambda: ddm_iinv(A_result, A, K))
  117. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  118. def test_Matrix_adjugate(name, A, A_inv, den):
  119. A = A.to_Matrix()
  120. A_inv = A_inv.to_Matrix()
  121. assert A.adjugate() == A_inv
  122. for method in ["bareiss", "berkowitz", "bird", "laplace", "lu"]:
  123. assert A.adjugate(method=method) == A_inv
  124. @pytest.mark.parametrize('name, A, A_inv, den', INVERSE_EXAMPLES)
  125. def test_dm_adj_det(name, A, A_inv, den):
  126. assert A.adj_det() == (A_inv, den)
  127. def test_inverse_inexact():
  128. M = Matrix([[x-0.3, -0.06, -0.22],
  129. [-0.46, x-0.48, -0.41],
  130. [-0.14, -0.39, x-0.64]])
  131. Mn = Matrix([[1.0*x**2 - 1.12*x + 0.1473, 0.06*x + 0.0474, 0.22*x - 0.081],
  132. [0.46*x - 0.237, 1.0*x**2 - 0.94*x + 0.1612, 0.41*x - 0.0218],
  133. [0.14*x + 0.1122, 0.39*x - 0.1086, 1.0*x**2 - 0.78*x + 0.1164]])
  134. d = 1.0*x**3 - 1.42*x**2 + 0.4249*x - 0.0546540000000002
  135. Mi = Mn / d
  136. M_dm = M.to_DM()
  137. M_dmd = M_dm.to_dense()
  138. M_dm_num, M_dm_den = M_dm.inv_den()
  139. M_dmd_num, M_dmd_den = M_dmd.inv_den()
  140. # XXX: We don't check M_dm().to_field().inv() which currently uses division
  141. # and produces a more complicate result from gcd cancellation failing.
  142. # DomainMatrix.inv() over RR(x) should be changed to clear denominators and
  143. # use DomainMatrix.inv_den().
  144. Minvs = [
  145. M.inv(),
  146. (M_dm_num.to_field() / M_dm_den).to_Matrix(),
  147. (M_dmd_num.to_field() / M_dmd_den).to_Matrix(),
  148. M_dm_num.to_Matrix() / M_dm_den.as_expr(),
  149. M_dmd_num.to_Matrix() / M_dmd_den.as_expr(),
  150. ]
  151. for Minv in Minvs:
  152. for Mi1, Mi2 in zip(Minv.flat(), Mi.flat()):
  153. assert all_close(Mi2, Mi1)