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"""Bounded-variable least-squares algorithm.""" 

2import numpy as np 

3from numpy.linalg import norm, lstsq 

4from scipy.optimize import OptimizeResult 

5 

6from .common import print_header_linear, print_iteration_linear 

7 

8 

9def compute_kkt_optimality(g, on_bound): 

10 """Compute the maximum violation of KKT conditions.""" 

11 g_kkt = g * on_bound 

12 free_set = on_bound == 0 

13 g_kkt[free_set] = np.abs(g[free_set]) 

14 return np.max(g_kkt) 

15 

16 

17def bvls(A, b, x_lsq, lb, ub, tol, max_iter, verbose): 

18 m, n = A.shape 

19 

20 x = x_lsq.copy() 

21 on_bound = np.zeros(n) 

22 

23 mask = x < lb 

24 x[mask] = lb[mask] 

25 on_bound[mask] = -1 

26 

27 mask = x > ub 

28 x[mask] = ub[mask] 

29 on_bound[mask] = 1 

30 

31 free_set = on_bound == 0 

32 active_set = ~free_set 

33 free_set, = np.nonzero(free_set) 

34 

35 r = A.dot(x) - b 

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

37 initial_cost = cost 

38 g = A.T.dot(r) 

39 

40 cost_change = None 

41 step_norm = None 

42 iteration = 0 

43 

44 if verbose == 2: 

45 print_header_linear() 

46 

47 # This is the initialization loop. The requirement is that the 

48 # least-squares solution on free variables is feasible before BVLS starts. 

49 # One possible initialization is to set all variables to lower or upper 

50 # bounds, but many iterations may be required from this state later on. 

51 # The implemented ad-hoc procedure which intuitively should give a better 

52 # initial state: find the least-squares solution on current free variables, 

53 # if its feasible then stop, otherwise, set violating variables to 

54 # corresponding bounds and continue on the reduced set of free variables. 

55 

56 while free_set.size > 0: 

57 if verbose == 2: 

58 optimality = compute_kkt_optimality(g, on_bound) 

59 print_iteration_linear(iteration, cost, cost_change, step_norm, 

60 optimality) 

61 

62 iteration += 1 

63 x_free_old = x[free_set].copy() 

64 

65 A_free = A[:, free_set] 

66 b_free = b - A.dot(x * active_set) 

67 z = lstsq(A_free, b_free, rcond=-1)[0] 

68 

69 lbv = z < lb[free_set] 

70 ubv = z > ub[free_set] 

71 v = lbv | ubv 

72 

73 if np.any(lbv): 

74 ind = free_set[lbv] 

75 x[ind] = lb[ind] 

76 active_set[ind] = True 

77 on_bound[ind] = -1 

78 

79 if np.any(ubv): 

80 ind = free_set[ubv] 

81 x[ind] = ub[ind] 

82 active_set[ind] = True 

83 on_bound[ind] = 1 

84 

85 ind = free_set[~v] 

86 x[ind] = z[~v] 

87 

88 r = A.dot(x) - b 

89 cost_new = 0.5 * np.dot(r, r) 

90 cost_change = cost - cost_new 

91 cost = cost_new 

92 g = A.T.dot(r) 

93 step_norm = norm(x[free_set] - x_free_old) 

94 

95 if np.any(v): 

96 free_set = free_set[~v] 

97 else: 

98 break 

99 

100 if max_iter is None: 

101 max_iter = n 

102 max_iter += iteration 

103 

104 termination_status = None 

105 

106 # Main BVLS loop. 

107 

108 optimality = compute_kkt_optimality(g, on_bound) 

109 for iteration in range(iteration, max_iter): 

110 if verbose == 2: 

111 print_iteration_linear(iteration, cost, cost_change, 

112 step_norm, optimality) 

113 

114 if optimality < tol: 

115 termination_status = 1 

116 

117 if termination_status is not None: 

118 break 

119 

120 move_to_free = np.argmax(g * on_bound) 

121 on_bound[move_to_free] = 0 

122 free_set = on_bound == 0 

123 active_set = ~free_set 

124 free_set, = np.nonzero(free_set) 

125 

126 x_free = x[free_set] 

127 x_free_old = x_free.copy() 

128 lb_free = lb[free_set] 

129 ub_free = ub[free_set] 

130 

131 A_free = A[:, free_set] 

132 b_free = b - A.dot(x * active_set) 

133 z = lstsq(A_free, b_free, rcond=-1)[0] 

134 

135 lbv, = np.nonzero(z < lb_free) 

136 ubv, = np.nonzero(z > ub_free) 

137 v = np.hstack((lbv, ubv)) 

138 

139 if v.size > 0: 

140 alphas = np.hstack(( 

141 lb_free[lbv] - x_free[lbv], 

142 ub_free[ubv] - x_free[ubv])) / (z[v] - x_free[v]) 

143 

144 i = np.argmin(alphas) 

145 i_free = v[i] 

146 alpha = alphas[i] 

147 

148 x_free *= 1 - alpha 

149 x_free += alpha * z 

150 

151 if i < lbv.size: 

152 on_bound[free_set[i_free]] = -1 

153 else: 

154 on_bound[free_set[i_free]] = 1 

155 else: 

156 x_free = z 

157 

158 x[free_set] = x_free 

159 step_norm = norm(x_free - x_free_old) 

160 

161 r = A.dot(x) - b 

162 cost_new = 0.5 * np.dot(r, r) 

163 cost_change = cost - cost_new 

164 

165 if cost_change < tol * cost: 

166 termination_status = 2 

167 cost = cost_new 

168 

169 g = A.T.dot(r) 

170 optimality = compute_kkt_optimality(g, on_bound) 

171 

172 if termination_status is None: 

173 termination_status = 0 

174 

175 return OptimizeResult( 

176 x=x, fun=r, cost=cost, optimality=optimality, active_mask=on_bound, 

177 nit=iteration + 1, status=termination_status, 

178 initial_cost=initial_cost)