Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.8/site-packages/scipy/optimize/_lsq/trf_linear.py: 8%
143 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-12 06:31 +0000
1"""The adaptation of Trust Region Reflective algorithm for a linear
2least-squares problem."""
3import numpy as np
4from numpy.linalg import norm
5from scipy.linalg import qr, solve_triangular
6from scipy.sparse.linalg import lsmr
7from scipy.optimize import OptimizeResult
9from .givens_elimination import givens_elimination
10from .common import (
11 EPS, step_size_to_bound, find_active_constraints, in_bounds,
12 make_strictly_feasible, build_quadratic_1d, evaluate_quadratic,
13 minimize_quadratic_1d, CL_scaling_vector, reflective_transformation,
14 print_header_linear, print_iteration_linear, compute_grad,
15 regularized_lsq_operator, right_multiplied_operator)
18def regularized_lsq_with_qr(m, n, R, QTb, perm, diag, copy_R=True):
19 """Solve regularized least squares using information from QR-decomposition.
21 The initial problem is to solve the following system in a least-squares
22 sense::
24 A x = b
25 D x = 0
27 where D is diagonal matrix. The method is based on QR decomposition
28 of the form A P = Q R, where P is a column permutation matrix, Q is an
29 orthogonal matrix and R is an upper triangular matrix.
31 Parameters
32 ----------
33 m, n : int
34 Initial shape of A.
35 R : ndarray, shape (n, n)
36 Upper triangular matrix from QR decomposition of A.
37 QTb : ndarray, shape (n,)
38 First n components of Q^T b.
39 perm : ndarray, shape (n,)
40 Array defining column permutation of A, such that ith column of
41 P is perm[i]-th column of identity matrix.
42 diag : ndarray, shape (n,)
43 Array containing diagonal elements of D.
45 Returns
46 -------
47 x : ndarray, shape (n,)
48 Found least-squares solution.
49 """
50 if copy_R:
51 R = R.copy()
52 v = QTb.copy()
54 givens_elimination(R, v, diag[perm])
56 abs_diag_R = np.abs(np.diag(R))
57 threshold = EPS * max(m, n) * np.max(abs_diag_R)
58 nns, = np.nonzero(abs_diag_R > threshold)
60 R = R[np.ix_(nns, nns)]
61 v = v[nns]
63 x = np.zeros(n)
64 x[perm[nns]] = solve_triangular(R, v)
66 return x
69def backtracking(A, g, x, p, theta, p_dot_g, lb, ub):
70 """Find an appropriate step size using backtracking line search."""
71 alpha = 1
72 while True:
73 x_new, _ = reflective_transformation(x + alpha * p, lb, ub)
74 step = x_new - x
75 cost_change = -evaluate_quadratic(A, g, step)
76 if cost_change > -0.1 * alpha * p_dot_g:
77 break
78 alpha *= 0.5
80 active = find_active_constraints(x_new, lb, ub)
81 if np.any(active != 0):
82 x_new, _ = reflective_transformation(x + theta * alpha * p, lb, ub)
83 x_new = make_strictly_feasible(x_new, lb, ub, rstep=0)
84 step = x_new - x
85 cost_change = -evaluate_quadratic(A, g, step)
87 return x, step, cost_change
90def select_step(x, A_h, g_h, c_h, p, p_h, d, lb, ub, theta):
91 """Select the best step according to Trust Region Reflective algorithm."""
92 if in_bounds(x + p, lb, ub):
93 return p
95 p_stride, hits = step_size_to_bound(x, p, lb, ub)
96 r_h = np.copy(p_h)
97 r_h[hits.astype(bool)] *= -1
98 r = d * r_h
100 # Restrict step, such that it hits the bound.
101 p *= p_stride
102 p_h *= p_stride
103 x_on_bound = x + p
105 # Find the step size along reflected direction.
106 r_stride_u, _ = step_size_to_bound(x_on_bound, r, lb, ub)
108 # Stay interior.
109 r_stride_l = (1 - theta) * r_stride_u
110 r_stride_u *= theta
112 if r_stride_u > 0:
113 a, b, c = build_quadratic_1d(A_h, g_h, r_h, s0=p_h, diag=c_h)
114 r_stride, r_value = minimize_quadratic_1d(
115 a, b, r_stride_l, r_stride_u, c=c)
116 r_h = p_h + r_h * r_stride
117 r = d * r_h
118 else:
119 r_value = np.inf
121 # Now correct p_h to make it strictly interior.
122 p_h *= theta
123 p *= theta
124 p_value = evaluate_quadratic(A_h, g_h, p_h, diag=c_h)
126 ag_h = -g_h
127 ag = d * ag_h
128 ag_stride_u, _ = step_size_to_bound(x, ag, lb, ub)
129 ag_stride_u *= theta
130 a, b = build_quadratic_1d(A_h, g_h, ag_h, diag=c_h)
131 ag_stride, ag_value = minimize_quadratic_1d(a, b, 0, ag_stride_u)
132 ag *= ag_stride
134 if p_value < r_value and p_value < ag_value:
135 return p
136 elif r_value < p_value and r_value < ag_value:
137 return r
138 else:
139 return ag
142def trf_linear(A, b, x_lsq, lb, ub, tol, lsq_solver, lsmr_tol,
143 max_iter, verbose, *, lsmr_maxiter=None):
144 m, n = A.shape
145 x, _ = reflective_transformation(x_lsq, lb, ub)
146 x = make_strictly_feasible(x, lb, ub, rstep=0.1)
148 if lsq_solver == 'exact':
149 QT, R, perm = qr(A, mode='economic', pivoting=True)
150 QT = QT.T
152 if m < n:
153 R = np.vstack((R, np.zeros((n - m, n))))
155 QTr = np.zeros(n)
156 k = min(m, n)
157 elif lsq_solver == 'lsmr':
158 r_aug = np.zeros(m + n)
159 auto_lsmr_tol = False
160 if lsmr_tol is None:
161 lsmr_tol = 1e-2 * tol
162 elif lsmr_tol == 'auto':
163 auto_lsmr_tol = True
165 r = A.dot(x) - b
166 g = compute_grad(A, r)
167 cost = 0.5 * np.dot(r, r)
168 initial_cost = cost
170 termination_status = None
171 step_norm = None
172 cost_change = None
174 if max_iter is None:
175 max_iter = 100
177 if verbose == 2:
178 print_header_linear()
180 for iteration in range(max_iter):
181 v, dv = CL_scaling_vector(x, g, lb, ub)
182 g_scaled = g * v
183 g_norm = norm(g_scaled, ord=np.inf)
184 if g_norm < tol:
185 termination_status = 1
187 if verbose == 2:
188 print_iteration_linear(iteration, cost, cost_change,
189 step_norm, g_norm)
191 if termination_status is not None:
192 break
194 diag_h = g * dv
195 diag_root_h = diag_h ** 0.5
196 d = v ** 0.5
197 g_h = d * g
199 A_h = right_multiplied_operator(A, d)
200 if lsq_solver == 'exact':
201 QTr[:k] = QT.dot(r)
202 p_h = -regularized_lsq_with_qr(m, n, R * d[perm], QTr, perm,
203 diag_root_h, copy_R=False)
204 elif lsq_solver == 'lsmr':
205 lsmr_op = regularized_lsq_operator(A_h, diag_root_h)
206 r_aug[:m] = r
207 if auto_lsmr_tol:
208 eta = 1e-2 * min(0.5, g_norm)
209 lsmr_tol = max(EPS, min(0.1, eta * g_norm))
210 p_h = -lsmr(lsmr_op, r_aug, maxiter=lsmr_maxiter,
211 atol=lsmr_tol, btol=lsmr_tol)[0]
213 p = d * p_h
215 p_dot_g = np.dot(p, g)
216 if p_dot_g > 0:
217 termination_status = -1
219 theta = 1 - min(0.005, g_norm)
220 step = select_step(x, A_h, g_h, diag_h, p, p_h, d, lb, ub, theta)
221 cost_change = -evaluate_quadratic(A, g, step)
223 # Perhaps almost never executed, the idea is that `p` is descent
224 # direction thus we must find acceptable cost decrease using simple
225 # "backtracking", otherwise the algorithm's logic would break.
226 if cost_change < 0:
227 x, step, cost_change = backtracking(
228 A, g, x, p, theta, p_dot_g, lb, ub)
229 else:
230 x = make_strictly_feasible(x + step, lb, ub, rstep=0)
232 step_norm = norm(step)
233 r = A.dot(x) - b
234 g = compute_grad(A, r)
236 if cost_change < tol * cost:
237 termination_status = 2
239 cost = 0.5 * np.dot(r, r)
241 if termination_status is None:
242 termination_status = 0
244 active_mask = find_active_constraints(x, lb, ub, rtol=tol)
246 return OptimizeResult(
247 x=x, fun=r, cost=cost, optimality=g_norm, active_mask=active_mask,
248 nit=iteration + 1, status=termination_status,
249 initial_cost=initial_cost)