subspaces.py 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174
  1. from .utilities import _iszero
  2. def _columnspace(M, simplify=False):
  3. """Returns a list of vectors (Matrix objects) that span columnspace of ``M``
  4. Examples
  5. ========
  6. >>> from sympy import Matrix
  7. >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
  8. >>> M
  9. Matrix([
  10. [ 1, 3, 0],
  11. [-2, -6, 0],
  12. [ 3, 9, 6]])
  13. >>> M.columnspace()
  14. [Matrix([
  15. [ 1],
  16. [-2],
  17. [ 3]]), Matrix([
  18. [0],
  19. [0],
  20. [6]])]
  21. See Also
  22. ========
  23. nullspace
  24. rowspace
  25. """
  26. reduced, pivots = M.echelon_form(simplify=simplify, with_pivots=True)
  27. return [M.col(i) for i in pivots]
  28. def _nullspace(M, simplify=False, iszerofunc=_iszero):
  29. """Returns list of vectors (Matrix objects) that span nullspace of ``M``
  30. Examples
  31. ========
  32. >>> from sympy import Matrix
  33. >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
  34. >>> M
  35. Matrix([
  36. [ 1, 3, 0],
  37. [-2, -6, 0],
  38. [ 3, 9, 6]])
  39. >>> M.nullspace()
  40. [Matrix([
  41. [-3],
  42. [ 1],
  43. [ 0]])]
  44. See Also
  45. ========
  46. columnspace
  47. rowspace
  48. """
  49. reduced, pivots = M.rref(iszerofunc=iszerofunc, simplify=simplify)
  50. free_vars = [i for i in range(M.cols) if i not in pivots]
  51. basis = []
  52. for free_var in free_vars:
  53. # for each free variable, we will set it to 1 and all others
  54. # to 0. Then, we will use back substitution to solve the system
  55. vec = [M.zero] * M.cols
  56. vec[free_var] = M.one
  57. for piv_row, piv_col in enumerate(pivots):
  58. vec[piv_col] -= reduced[piv_row, free_var]
  59. basis.append(vec)
  60. return [M._new(M.cols, 1, b) for b in basis]
  61. def _rowspace(M, simplify=False):
  62. """Returns a list of vectors that span the row space of ``M``.
  63. Examples
  64. ========
  65. >>> from sympy import Matrix
  66. >>> M = Matrix(3, 3, [1, 3, 0, -2, -6, 0, 3, 9, 6])
  67. >>> M
  68. Matrix([
  69. [ 1, 3, 0],
  70. [-2, -6, 0],
  71. [ 3, 9, 6]])
  72. >>> M.rowspace()
  73. [Matrix([[1, 3, 0]]), Matrix([[0, 0, 6]])]
  74. """
  75. reduced, pivots = M.echelon_form(simplify=simplify, with_pivots=True)
  76. return [reduced.row(i) for i in range(len(pivots))]
  77. def _orthogonalize(cls, *vecs, normalize=False, rankcheck=False):
  78. """Apply the Gram-Schmidt orthogonalization procedure
  79. to vectors supplied in ``vecs``.
  80. Parameters
  81. ==========
  82. vecs
  83. vectors to be made orthogonal
  84. normalize : bool
  85. If ``True``, return an orthonormal basis.
  86. rankcheck : bool
  87. If ``True``, the computation does not stop when encountering
  88. linearly dependent vectors.
  89. If ``False``, it will raise ``ValueError`` when any zero
  90. or linearly dependent vectors are found.
  91. Returns
  92. =======
  93. list
  94. List of orthogonal (or orthonormal) basis vectors.
  95. Examples
  96. ========
  97. >>> from sympy import I, Matrix
  98. >>> v = [Matrix([1, I]), Matrix([1, -I])]
  99. >>> Matrix.orthogonalize(*v)
  100. [Matrix([
  101. [1],
  102. [I]]), Matrix([
  103. [ 1],
  104. [-I]])]
  105. See Also
  106. ========
  107. MatrixBase.QRdecomposition
  108. References
  109. ==========
  110. .. [1] https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
  111. """
  112. from .decompositions import _QRdecomposition_optional
  113. if not vecs:
  114. return []
  115. all_row_vecs = (vecs[0].rows == 1)
  116. vecs = [x.vec() for x in vecs]
  117. M = cls.hstack(*vecs)
  118. Q, R = _QRdecomposition_optional(M, normalize=normalize)
  119. if rankcheck and Q.cols < len(vecs):
  120. raise ValueError("GramSchmidt: vector set not linearly independent")
  121. ret = []
  122. for i in range(Q.cols):
  123. if all_row_vecs:
  124. col = cls(Q[:, i].T)
  125. else:
  126. col = cls(Q[:, i])
  127. ret.append(col)
  128. return ret