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""" 

2Holds files for l1 regularization of LikelihoodModel, using 

3scipy.optimize.slsqp 

4""" 

5import numpy as np 

6from scipy.optimize import fmin_slsqp 

7import statsmodels.base.l1_solvers_common as l1_solvers_common 

8 

9 

10def fit_l1_slsqp( 

11 f, score, start_params, args, kwargs, disp=False, maxiter=1000, 

12 callback=None, retall=False, full_output=False, hess=None): 

13 """ 

14 Solve the l1 regularized problem using scipy.optimize.fmin_slsqp(). 

15 

16 Specifically: We convert the convex but non-smooth problem 

17 

18 .. math:: \\min_\\beta f(\\beta) + \\sum_k\\alpha_k |\\beta_k| 

19 

20 via the transformation to the smooth, convex, constrained problem in twice 

21 as many variables (adding the "added variables" :math:`u_k`) 

22 

23 .. math:: \\min_{\\beta,u} f(\\beta) + \\sum_k\\alpha_k u_k, 

24 

25 subject to 

26 

27 .. math:: -u_k \\leq \\beta_k \\leq u_k. 

28 

29 Parameters 

30 ---------- 

31 All the usual parameters from LikelhoodModel.fit 

32 alpha : non-negative scalar or numpy array (same size as parameters) 

33 The weight multiplying the l1 penalty term 

34 trim_mode : 'auto, 'size', or 'off' 

35 If not 'off', trim (set to zero) parameters that would have been zero 

36 if the solver reached the theoretical minimum. 

37 If 'auto', trim params using the Theory above. 

38 If 'size', trim params if they have very small absolute value 

39 size_trim_tol : float or 'auto' (default = 'auto') 

40 For use when trim_mode === 'size' 

41 auto_trim_tol : float 

42 For sue when trim_mode == 'auto'. Use 

43 qc_tol : float 

44 Print warning and do not allow auto trim when (ii) in "Theory" (above) 

45 is violated by this much. 

46 qc_verbose : bool 

47 If true, print out a full QC report upon failure 

48 acc : float (default 1e-6) 

49 Requested accuracy as used by slsqp 

50 """ 

51 start_params = np.array(start_params).ravel('F') 

52 

53 ### Extract values 

54 # k_params is total number of covariates, 

55 # possibly including a leading constant. 

56 k_params = len(start_params) 

57 # The start point 

58 x0 = np.append(start_params, np.fabs(start_params)) 

59 # alpha is the regularization parameter 

60 alpha = np.array(kwargs['alpha_rescaled']).ravel('F') 

61 # Make sure it's a vector 

62 alpha = alpha * np.ones(k_params) 

63 assert alpha.min() >= 0 

64 # Convert display parameters to scipy.optimize form 

65 disp_slsqp = _get_disp_slsqp(disp, retall) 

66 # Set/retrieve the desired accuracy 

67 acc = kwargs.setdefault('acc', 1e-10) 

68 

69 ### Wrap up for use in fmin_slsqp 

70 func = lambda x_full: _objective_func(f, x_full, k_params, alpha, *args) 

71 f_ieqcons_wrap = lambda x_full: _f_ieqcons(x_full, k_params) 

72 fprime_wrap = lambda x_full: _fprime(score, x_full, k_params, alpha) 

73 fprime_ieqcons_wrap = lambda x_full: _fprime_ieqcons(x_full, k_params) 

74 

75 ### Call the solver 

76 results = fmin_slsqp( 

77 func, x0, f_ieqcons=f_ieqcons_wrap, fprime=fprime_wrap, acc=acc, 

78 iter=maxiter, disp=disp_slsqp, full_output=full_output, 

79 fprime_ieqcons=fprime_ieqcons_wrap) 

80 params = np.asarray(results[0][:k_params]) 

81 

82 ### Post-process 

83 # QC 

84 qc_tol = kwargs['qc_tol'] 

85 qc_verbose = kwargs['qc_verbose'] 

86 passed = l1_solvers_common.qc_results( 

87 params, alpha, score, qc_tol, qc_verbose) 

88 # Possibly trim 

89 trim_mode = kwargs['trim_mode'] 

90 size_trim_tol = kwargs['size_trim_tol'] 

91 auto_trim_tol = kwargs['auto_trim_tol'] 

92 params, trimmed = l1_solvers_common.do_trim_params( 

93 params, k_params, alpha, score, passed, trim_mode, size_trim_tol, 

94 auto_trim_tol) 

95 

96 ### Pack up return values for statsmodels optimizers 

97 # TODO These retvals are returned as mle_retvals...but the fit was not ML. 

98 # This could be confusing someday. 

99 if full_output: 

100 x_full, fx, its, imode, smode = results 

101 fopt = func(np.asarray(x_full)) 

102 converged = (imode == 0) 

103 warnflag = str(imode) + ' ' + smode 

104 iterations = its 

105 gopt = float('nan') # Objective is non-differentiable 

106 hopt = float('nan') 

107 retvals = { 

108 'fopt': fopt, 'converged': converged, 'iterations': iterations, 

109 'gopt': gopt, 'hopt': hopt, 'trimmed': trimmed, 

110 'warnflag': warnflag} 

111 

112 ### Return 

113 if full_output: 

114 return params, retvals 

115 else: 

116 return params 

117 

118 

119def _get_disp_slsqp(disp, retall): 

120 if disp or retall: 

121 if disp: 

122 disp_slsqp = 1 

123 if retall: 

124 disp_slsqp = 2 

125 else: 

126 disp_slsqp = 0 

127 return disp_slsqp 

128 

129 

130def _objective_func(f, x_full, k_params, alpha, *args): 

131 """ 

132 The regularized objective function 

133 """ 

134 x_params = x_full[:k_params] 

135 x_added = x_full[k_params:] 

136 ## Return 

137 return f(x_params, *args) + (alpha * x_added).sum() 

138 

139 

140def _fprime(score, x_full, k_params, alpha): 

141 """ 

142 The regularized derivative 

143 """ 

144 x_params = x_full[:k_params] 

145 # The derivative just appends a vector of constants 

146 return np.append(score(x_params), alpha) 

147 

148 

149def _f_ieqcons(x_full, k_params): 

150 """ 

151 The inequality constraints. 

152 """ 

153 x_params = x_full[:k_params] 

154 x_added = x_full[k_params:] 

155 # All entries in this vector must be \geq 0 in a feasible solution 

156 return np.append(x_params + x_added, x_added - x_params) 

157 

158 

159def _fprime_ieqcons(x_full, k_params): 

160 """ 

161 Derivative of the inequality constraints 

162 """ 

163 I = np.eye(k_params) # noqa:E741 

164 A = np.concatenate((I, I), axis=1) 

165 B = np.concatenate((-I, I), axis=1) 

166 C = np.concatenate((A, B), axis=0) 

167 ## Return 

168 return C