_norm.py 6.0 KB

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