_nnls.py 2.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596
  1. import numpy as np
  2. from ._slsqplib import nnls as _nnls
  3. from scipy._lib.deprecation import _deprecate_positional_args, _NoValue
  4. __all__ = ['nnls']
  5. @_deprecate_positional_args(version='1.18.0',
  6. deprecated_args={'atol'})
  7. def nnls(A, b, *, maxiter=None, atol=_NoValue):
  8. """
  9. Solve ``argmin_x || Ax - b ||_2^2`` for ``x>=0``.
  10. This problem, often called as NonNegative Least Squares, is a convex
  11. optimization problem with convex constraints. It typically arises when
  12. the ``x`` models quantities for which only nonnegative values are
  13. attainable; weight of ingredients, component costs and so on.
  14. Parameters
  15. ----------
  16. A : (m, n) ndarray
  17. Coefficient array
  18. b : (m,) ndarray, float
  19. Right-hand side vector.
  20. maxiter: int, optional
  21. Maximum number of iterations, optional. Default value is ``3 * n``.
  22. atol : float, optional
  23. .. deprecated:: 1.18.0
  24. This parameter is deprecated and will be removed in SciPy 1.18.0.
  25. It is not used in the implementation.
  26. Returns
  27. -------
  28. x : ndarray
  29. Solution vector.
  30. rnorm : float
  31. The 2-norm of the residual, ``|| Ax-b ||_2``.
  32. See Also
  33. --------
  34. lsq_linear : Linear least squares with bounds on the variables
  35. Notes
  36. -----
  37. The code is based on the classical algorithm of [1]_. It utilizes an active
  38. set method and solves the KKT (Karush-Kuhn-Tucker) conditions for the
  39. non-negative least squares problem.
  40. References
  41. ----------
  42. .. [1] : Lawson C., Hanson R.J., "Solving Least Squares Problems", SIAM,
  43. 1995, :doi:`10.1137/1.9781611971217`
  44. Examples
  45. --------
  46. >>> import numpy as np
  47. >>> from scipy.optimize import nnls
  48. ...
  49. >>> A = np.array([[1, 0], [1, 0], [0, 1]])
  50. >>> b = np.array([2, 1, 1])
  51. >>> nnls(A, b)
  52. (array([1.5, 1. ]), 0.7071067811865475)
  53. >>> b = np.array([-1, -1, -1])
  54. >>> nnls(A, b)
  55. (array([0., 0.]), 1.7320508075688772)
  56. """
  57. A = np.asarray_chkfinite(A, dtype=np.float64, order='C')
  58. b = np.asarray_chkfinite(b, dtype=np.float64)
  59. if len(A.shape) != 2:
  60. raise ValueError(f"Expected a 2D array, but the shape of A is {A.shape}")
  61. if (b.ndim > 2) or ((b.ndim == 2) and (b.shape[1] != 1)):
  62. raise ValueError("Expected a 1D array,(or 2D with one column), but the,"
  63. f" shape of b is {b.shape}")
  64. elif (b.ndim == 2) and (b.shape[1] == 1):
  65. b = b.ravel()
  66. m, n = A.shape
  67. if m != b.shape[0]:
  68. raise ValueError(
  69. "Incompatible dimensions. The first dimension of " +
  70. f"A is {m}, while the shape of b is {(b.shape[0], )}")
  71. if not maxiter:
  72. maxiter = 3*n
  73. x, rnorm, info = _nnls(A, b, maxiter)
  74. if info == 3:
  75. raise RuntimeError("Maximum number of iterations reached.")
  76. return x, rnorm