test_lsq_linear.py 4.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150
  1. import numpy as np
  2. from numpy.linalg import lstsq
  3. from numpy.testing import assert_allclose, assert_equal, assert_
  4. from pytest import raises as assert_raises
  5. from scipy.sparse import rand
  6. from scipy.sparse.linalg import aslinearoperator
  7. from scipy.optimize import lsq_linear
  8. A = np.array([
  9. [0.171, -0.057],
  10. [-0.049, -0.248],
  11. [-0.166, 0.054],
  12. ])
  13. b = np.array([0.074, 1.014, -0.383])
  14. class BaseMixin(object):
  15. def setup_method(self):
  16. self.rnd = np.random.RandomState(0)
  17. def test_dense_no_bounds(self):
  18. for lsq_solver in self.lsq_solvers:
  19. res = lsq_linear(A, b, method=self.method, lsq_solver=lsq_solver)
  20. assert_allclose(res.x, lstsq(A, b, rcond=-1)[0])
  21. def test_dense_bounds(self):
  22. # Solutions for comparison are taken from MATLAB.
  23. lb = np.array([-1, -10])
  24. ub = np.array([1, 0])
  25. for lsq_solver in self.lsq_solvers:
  26. res = lsq_linear(A, b, (lb, ub), method=self.method,
  27. lsq_solver=lsq_solver)
  28. assert_allclose(res.x, lstsq(A, b, rcond=-1)[0])
  29. lb = np.array([0.0, -np.inf])
  30. for lsq_solver in self.lsq_solvers:
  31. res = lsq_linear(A, b, (lb, np.inf), method=self.method,
  32. lsq_solver=lsq_solver)
  33. assert_allclose(res.x, np.array([0.0, -4.084174437334673]),
  34. atol=1e-6)
  35. lb = np.array([-1, 0])
  36. for lsq_solver in self.lsq_solvers:
  37. res = lsq_linear(A, b, (lb, np.inf), method=self.method,
  38. lsq_solver=lsq_solver)
  39. assert_allclose(res.x, np.array([0.448427311733504, 0]),
  40. atol=1e-15)
  41. ub = np.array([np.inf, -5])
  42. for lsq_solver in self.lsq_solvers:
  43. res = lsq_linear(A, b, (-np.inf, ub), method=self.method,
  44. lsq_solver=lsq_solver)
  45. assert_allclose(res.x, np.array([-0.105560998682388, -5]))
  46. ub = np.array([-1, np.inf])
  47. for lsq_solver in self.lsq_solvers:
  48. res = lsq_linear(A, b, (-np.inf, ub), method=self.method,
  49. lsq_solver=lsq_solver)
  50. assert_allclose(res.x, np.array([-1, -4.181102129483254]))
  51. lb = np.array([0, -4])
  52. ub = np.array([1, 0])
  53. for lsq_solver in self.lsq_solvers:
  54. res = lsq_linear(A, b, (lb, ub), method=self.method,
  55. lsq_solver=lsq_solver)
  56. assert_allclose(res.x, np.array([0.005236663400791, -4]))
  57. def test_dense_rank_deficient(self):
  58. A = np.array([[-0.307, -0.184]])
  59. b = np.array([0.773])
  60. lb = [-0.1, -0.1]
  61. ub = [0.1, 0.1]
  62. for lsq_solver in self.lsq_solvers:
  63. res = lsq_linear(A, b, (lb, ub), method=self.method,
  64. lsq_solver=lsq_solver)
  65. assert_allclose(res.x, [-0.1, -0.1])
  66. A = np.array([
  67. [0.334, 0.668],
  68. [-0.516, -1.032],
  69. [0.192, 0.384],
  70. ])
  71. b = np.array([-1.436, 0.135, 0.909])
  72. lb = [0, -1]
  73. ub = [1, -0.5]
  74. for lsq_solver in self.lsq_solvers:
  75. res = lsq_linear(A, b, (lb, ub), method=self.method,
  76. lsq_solver=lsq_solver)
  77. assert_allclose(res.optimality, 0, atol=1e-11)
  78. def test_full_result(self):
  79. lb = np.array([0, -4])
  80. ub = np.array([1, 0])
  81. res = lsq_linear(A, b, (lb, ub), method=self.method)
  82. assert_allclose(res.x, [0.005236663400791, -4])
  83. r = A.dot(res.x) - b
  84. assert_allclose(res.cost, 0.5 * np.dot(r, r))
  85. assert_allclose(res.fun, r)
  86. assert_allclose(res.optimality, 0.0, atol=1e-12)
  87. assert_equal(res.active_mask, [0, -1])
  88. assert_(res.nit < 15)
  89. assert_(res.status == 1 or res.status == 3)
  90. assert_(isinstance(res.message, str))
  91. assert_(res.success)
  92. class SparseMixin(object):
  93. def test_sparse_and_LinearOperator(self):
  94. m = 5000
  95. n = 1000
  96. A = rand(m, n, random_state=0)
  97. b = self.rnd.randn(m)
  98. res = lsq_linear(A, b)
  99. assert_allclose(res.optimality, 0, atol=1e-6)
  100. A = aslinearoperator(A)
  101. res = lsq_linear(A, b)
  102. assert_allclose(res.optimality, 0, atol=1e-6)
  103. def test_sparse_bounds(self):
  104. m = 5000
  105. n = 1000
  106. A = rand(m, n, random_state=0)
  107. b = self.rnd.randn(m)
  108. lb = self.rnd.randn(n)
  109. ub = lb + 1
  110. res = lsq_linear(A, b, (lb, ub))
  111. assert_allclose(res.optimality, 0.0, atol=1e-6)
  112. res = lsq_linear(A, b, (lb, ub), lsmr_tol=1e-13)
  113. assert_allclose(res.optimality, 0.0, atol=1e-6)
  114. res = lsq_linear(A, b, (lb, ub), lsmr_tol='auto')
  115. assert_allclose(res.optimality, 0.0, atol=1e-6)
  116. class TestTRF(BaseMixin, SparseMixin):
  117. method = 'trf'
  118. lsq_solvers = ['exact', 'lsmr']
  119. class TestBVLS(BaseMixin):
  120. method = 'bvls'
  121. lsq_solvers = ['exact']