Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

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 

8 

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) 

16 

17 

18def regularized_lsq_with_qr(m, n, R, QTb, perm, diag, copy_R=True): 

19 """Solve regularized least squares using information from QR-decomposition. 

20 

21 The initial problem is to solve the following system in a least-squares 

22 sense: 

23 :: 

24 

25 A x = b 

26 D x = 0 

27 

28 where D is diagonal matrix. The method is based on QR decomposition 

29 of the form A P = Q R, where P is a column permutation matrix, Q is an 

30 orthogonal matrix and R is an upper triangular matrix. 

31 

32 Parameters 

33 ---------- 

34 m, n : int 

35 Initial shape of A. 

36 R : ndarray, shape (n, n) 

37 Upper triangular matrix from QR decomposition of A. 

38 QTb : ndarray, shape (n,) 

39 First n components of Q^T b. 

40 perm : ndarray, shape (n,) 

41 Array defining column permutation of A, such that ith column of 

42 P is perm[i]-th column of identity matrix. 

43 diag : ndarray, shape (n,) 

44 Array containing diagonal elements of D. 

45 

46 Returns 

47 ------- 

48 x : ndarray, shape (n,) 

49 Found least-squares solution. 

50 """ 

51 if copy_R: 

52 R = R.copy() 

53 v = QTb.copy() 

54 

55 givens_elimination(R, v, diag[perm]) 

56 

57 abs_diag_R = np.abs(np.diag(R)) 

58 threshold = EPS * max(m, n) * np.max(abs_diag_R) 

59 nns, = np.nonzero(abs_diag_R > threshold) 

60 

61 R = R[np.ix_(nns, nns)] 

62 v = v[nns] 

63 

64 x = np.zeros(n) 

65 x[perm[nns]] = solve_triangular(R, v) 

66 

67 return x 

68 

69 

70def backtracking(A, g, x, p, theta, p_dot_g, lb, ub): 

71 """Find an appropriate step size using backtracking line search.""" 

72 alpha = 1 

73 while True: 

74 x_new, _ = reflective_transformation(x + alpha * p, lb, ub) 

75 step = x_new - x 

76 cost_change = -evaluate_quadratic(A, g, step) 

77 if cost_change > -0.1 * alpha * p_dot_g: 

78 break 

79 alpha *= 0.5 

80 

81 active = find_active_constraints(x_new, lb, ub) 

82 if np.any(active != 0): 

83 x_new, _ = reflective_transformation(x + theta * alpha * p, lb, ub) 

84 x_new = make_strictly_feasible(x_new, lb, ub, rstep=0) 

85 step = x_new - x 

86 cost_change = -evaluate_quadratic(A, g, step) 

87 

88 return x, step, cost_change 

89 

90 

91def select_step(x, A_h, g_h, c_h, p, p_h, d, lb, ub, theta): 

92 """Select the best step according to Trust Region Reflective algorithm.""" 

93 if in_bounds(x + p, lb, ub): 

94 return p 

95 

96 p_stride, hits = step_size_to_bound(x, p, lb, ub) 

97 r_h = np.copy(p_h) 

98 r_h[hits.astype(bool)] *= -1 

99 r = d * r_h 

100 

101 # Restrict step, such that it hits the bound. 

102 p *= p_stride 

103 p_h *= p_stride 

104 x_on_bound = x + p 

105 

106 # Find the step size along reflected direction. 

107 r_stride_u, _ = step_size_to_bound(x_on_bound, r, lb, ub) 

108 

109 # Stay interior. 

110 r_stride_l = (1 - theta) * r_stride_u 

111 r_stride_u *= theta 

112 

113 if r_stride_u > 0: 

114 a, b, c = build_quadratic_1d(A_h, g_h, r_h, s0=p_h, diag=c_h) 

115 r_stride, r_value = minimize_quadratic_1d( 

116 a, b, r_stride_l, r_stride_u, c=c) 

117 r_h = p_h + r_h * r_stride 

118 r = d * r_h 

119 else: 

120 r_value = np.inf 

121 

122 # Now correct p_h to make it strictly interior. 

123 p_h *= theta 

124 p *= theta 

125 p_value = evaluate_quadratic(A_h, g_h, p_h, diag=c_h) 

126 

127 ag_h = -g_h 

128 ag = d * ag_h 

129 ag_stride_u, _ = step_size_to_bound(x, ag, lb, ub) 

130 ag_stride_u *= theta 

131 a, b = build_quadratic_1d(A_h, g_h, ag_h, diag=c_h) 

132 ag_stride, ag_value = minimize_quadratic_1d(a, b, 0, ag_stride_u) 

133 ag *= ag_stride 

134 

135 if p_value < r_value and p_value < ag_value: 

136 return p 

137 elif r_value < p_value and r_value < ag_value: 

138 return r 

139 else: 

140 return ag 

141 

142 

143def trf_linear(A, b, x_lsq, lb, ub, tol, lsq_solver, lsmr_tol, max_iter, 

144 verbose): 

145 m, n = A.shape 

146 x, _ = reflective_transformation(x_lsq, lb, ub) 

147 x = make_strictly_feasible(x, lb, ub, rstep=0.1) 

148 

149 if lsq_solver == 'exact': 

150 QT, R, perm = qr(A, mode='economic', pivoting=True) 

151 QT = QT.T 

152 

153 if m < n: 

154 R = np.vstack((R, np.zeros((n - m, n)))) 

155 

156 QTr = np.zeros(n) 

157 k = min(m, n) 

158 elif lsq_solver == 'lsmr': 

159 r_aug = np.zeros(m + n) 

160 auto_lsmr_tol = False 

161 if lsmr_tol is None: 

162 lsmr_tol = 1e-2 * tol 

163 elif lsmr_tol == 'auto': 

164 auto_lsmr_tol = True 

165 

166 r = A.dot(x) - b 

167 g = compute_grad(A, r) 

168 cost = 0.5 * np.dot(r, r) 

169 initial_cost = cost 

170 

171 termination_status = None 

172 step_norm = None 

173 cost_change = None 

174 

175 if max_iter is None: 

176 max_iter = 100 

177 

178 if verbose == 2: 

179 print_header_linear() 

180 

181 for iteration in range(max_iter): 

182 v, dv = CL_scaling_vector(x, g, lb, ub) 

183 g_scaled = g * v 

184 g_norm = norm(g_scaled, ord=np.inf) 

185 if g_norm < tol: 

186 termination_status = 1 

187 

188 if verbose == 2: 

189 print_iteration_linear(iteration, cost, cost_change, 

190 step_norm, g_norm) 

191 

192 if termination_status is not None: 

193 break 

194 

195 diag_h = g * dv 

196 diag_root_h = diag_h ** 0.5 

197 d = v ** 0.5 

198 g_h = d * g 

199 

200 A_h = right_multiplied_operator(A, d) 

201 if lsq_solver == 'exact': 

202 QTr[:k] = QT.dot(r) 

203 p_h = -regularized_lsq_with_qr(m, n, R * d[perm], QTr, perm, 

204 diag_root_h, copy_R=False) 

205 elif lsq_solver == 'lsmr': 

206 lsmr_op = regularized_lsq_operator(A_h, diag_root_h) 

207 r_aug[:m] = r 

208 if auto_lsmr_tol: 

209 eta = 1e-2 * min(0.5, g_norm) 

210 lsmr_tol = max(EPS, min(0.1, eta * g_norm)) 

211 p_h = -lsmr(lsmr_op, r_aug, atol=lsmr_tol, btol=lsmr_tol)[0] 

212 

213 p = d * p_h 

214 

215 p_dot_g = np.dot(p, g) 

216 if p_dot_g > 0: 

217 termination_status = -1 

218 

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) 

222 

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) 

231 

232 step_norm = norm(step) 

233 r = A.dot(x) - b 

234 g = compute_grad(A, r) 

235 

236 if cost_change < tol * cost: 

237 termination_status = 2 

238 

239 cost = 0.5 * np.dot(r, r) 

240 

241 if termination_status is None: 

242 termination_status = 0 

243 

244 active_mask = find_active_constraints(x, lb, ub, rtol=tol) 

245 

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)