_norm.py 5.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193
  1. """Sparse matrix norms.
  2. """
  3. import numpy as np
  4. from scipy.sparse import issparse
  5. from scipy.sparse.linalg import svds
  6. import scipy.sparse as sp
  7. from numpy import Inf, sqrt, abs
  8. __all__ = ['norm']
  9. def _sparse_frobenius_norm(x):
  10. data = sp._sputils._todata(x)
  11. return np.linalg.norm(data)
  12. def norm(x, ord=None, axis=None):
  13. """
  14. Norm of a sparse matrix
  15. This function is able to return one of seven different matrix norms,
  16. depending on the value of the ``ord`` parameter.
  17. Parameters
  18. ----------
  19. x : a sparse matrix
  20. Input sparse matrix.
  21. ord : {non-zero int, inf, -inf, 'fro'}, optional
  22. Order of the norm (see table under ``Notes``). inf means numpy's
  23. `inf` object.
  24. axis : {int, 2-tuple of ints, None}, optional
  25. If `axis` is an integer, it specifies the axis of `x` along which to
  26. compute the vector norms. If `axis` is a 2-tuple, it specifies the
  27. axes that hold 2-D matrices, and the matrix norms of these matrices
  28. are computed. If `axis` is None then either a vector norm (when `x`
  29. is 1-D) or a matrix norm (when `x` is 2-D) is returned.
  30. Returns
  31. -------
  32. n : float or ndarray
  33. Notes
  34. -----
  35. Some of the ord are not implemented because some associated functions like,
  36. _multi_svd_norm, are not yet available for sparse matrix.
  37. This docstring is modified based on numpy.linalg.norm.
  38. https://github.com/numpy/numpy/blob/main/numpy/linalg/linalg.py
  39. The following norms can be calculated:
  40. ===== ============================
  41. ord norm for sparse matrices
  42. ===== ============================
  43. None Frobenius norm
  44. 'fro' Frobenius norm
  45. inf max(sum(abs(x), axis=1))
  46. -inf min(sum(abs(x), axis=1))
  47. 0 abs(x).sum(axis=axis)
  48. 1 max(sum(abs(x), axis=0))
  49. -1 min(sum(abs(x), axis=0))
  50. 2 Spectral norm (the largest singular value)
  51. -2 Not implemented
  52. other Not implemented
  53. ===== ============================
  54. The Frobenius norm is given by [1]_:
  55. :math:`||A||_F = [\\sum_{i,j} abs(a_{i,j})^2]^{1/2}`
  56. References
  57. ----------
  58. .. [1] G. H. Golub and C. F. Van Loan, *Matrix Computations*,
  59. Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15
  60. Examples
  61. --------
  62. >>> from scipy.sparse import *
  63. >>> import numpy as np
  64. >>> from scipy.sparse.linalg import norm
  65. >>> a = np.arange(9) - 4
  66. >>> a
  67. array([-4, -3, -2, -1, 0, 1, 2, 3, 4])
  68. >>> b = a.reshape((3, 3))
  69. >>> b
  70. array([[-4, -3, -2],
  71. [-1, 0, 1],
  72. [ 2, 3, 4]])
  73. >>> b = csr_matrix(b)
  74. >>> norm(b)
  75. 7.745966692414834
  76. >>> norm(b, 'fro')
  77. 7.745966692414834
  78. >>> norm(b, np.inf)
  79. 9
  80. >>> norm(b, -np.inf)
  81. 2
  82. >>> norm(b, 1)
  83. 7
  84. >>> norm(b, -1)
  85. 6
  86. The matrix 2-norm or the spectral norm is the largest singular
  87. value, computed approximately and with limitations.
  88. >>> b = diags([-1, 1], [0, 1], shape=(9, 10))
  89. >>> norm(b, 2)
  90. 1.9753...
  91. """
  92. if not issparse(x):
  93. raise TypeError("input is not sparse. use numpy.linalg.norm")
  94. # Check the default case first and handle it immediately.
  95. if axis is None and ord in (None, 'fro', 'f'):
  96. return _sparse_frobenius_norm(x)
  97. # Some norms require functions that are not implemented for all types.
  98. x = x.tocsr()
  99. if axis is None:
  100. axis = (0, 1)
  101. elif not isinstance(axis, tuple):
  102. msg = "'axis' must be None, an integer or a tuple of integers"
  103. try:
  104. int_axis = int(axis)
  105. except TypeError as e:
  106. raise TypeError(msg) from e
  107. if axis != int_axis:
  108. raise TypeError(msg)
  109. axis = (int_axis,)
  110. nd = 2
  111. if len(axis) == 2:
  112. row_axis, col_axis = axis
  113. if not (-nd <= row_axis < nd and -nd <= col_axis < nd):
  114. raise ValueError('Invalid axis %r for an array with shape %r' %
  115. (axis, x.shape))
  116. if row_axis % nd == col_axis % nd:
  117. raise ValueError('Duplicate axes given.')
  118. if ord == 2:
  119. # Only solver="lobpcg" supports all numpy dtypes
  120. _, s, _ = svds(x, k=1, solver="lobpcg")
  121. return s[0]
  122. elif ord == -2:
  123. raise NotImplementedError
  124. #return _multi_svd_norm(x, row_axis, col_axis, amin)
  125. elif ord == 1:
  126. return abs(x).sum(axis=row_axis).max(axis=col_axis)[0,0]
  127. elif ord == Inf:
  128. return abs(x).sum(axis=col_axis).max(axis=row_axis)[0,0]
  129. elif ord == -1:
  130. return abs(x).sum(axis=row_axis).min(axis=col_axis)[0,0]
  131. elif ord == -Inf:
  132. return abs(x).sum(axis=col_axis).min(axis=row_axis)[0,0]
  133. elif ord in (None, 'f', 'fro'):
  134. # The axis order does not matter for this norm.
  135. return _sparse_frobenius_norm(x)
  136. else:
  137. raise ValueError("Invalid norm order for matrices.")
  138. elif len(axis) == 1:
  139. a, = axis
  140. if not (-nd <= a < nd):
  141. raise ValueError('Invalid axis %r for an array with shape %r' %
  142. (axis, x.shape))
  143. if ord == Inf:
  144. M = abs(x).max(axis=a)
  145. elif ord == -Inf:
  146. M = abs(x).min(axis=a)
  147. elif ord == 0:
  148. # Zero norm
  149. M = (x != 0).sum(axis=a)
  150. elif ord == 1:
  151. # special case for speedup
  152. M = abs(x).sum(axis=a)
  153. elif ord in (2, None):
  154. M = sqrt(abs(x).power(2).sum(axis=a))
  155. else:
  156. try:
  157. ord + 1
  158. except TypeError as e:
  159. raise ValueError('Invalid norm order for vectors.') from e
  160. M = np.power(abs(x).power(ord).sum(axis=a), 1 / ord)
  161. if hasattr(M, 'toarray'):
  162. return M.toarray().ravel()
  163. elif hasattr(M, 'A'):
  164. return M.A.ravel()
  165. else:
  166. return M.ravel()
  167. else:
  168. raise ValueError("Improper number of dimensions to norm.")