_norm.py 5.7 KB

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