36 const std::string& update_type)
39 , learning_rate_decr(0.5)
40 , update_type{ update_type }
46 void run(Eigen::VectorXd& beta0,
47 Eigen::VectorXd& beta,
49 const Eigen::ArrayXd& lambda,
50 const std::unique_ptr<Loss>& loss,
52 const Eigen::VectorXd& gradient,
53 const std::vector<int>& working_set,
54 const Eigen::MatrixXd& x,
55 const Eigen::VectorXd& x_centers,
56 const Eigen::VectorXd& x_scales,
57 const Eigen::MatrixXd& y)
override;
60 void run(Eigen::VectorXd& beta0,
61 Eigen::VectorXd& beta,
63 const Eigen::ArrayXd& lambda,
64 const std::unique_ptr<Loss>& loss,
66 const Eigen::VectorXd& gradient,
67 const std::vector<int>& working_set,
68 const Eigen::SparseMatrix<double>& x,
69 const Eigen::VectorXd& x_centers,
70 const Eigen::VectorXd& x_scales,
71 const Eigen::MatrixXd& y)
override;
74 template<
typename MatrixType>
75 void runImpl(Eigen::VectorXd& beta0,
76 Eigen::VectorXd& beta,
78 const Eigen::ArrayXd& lambda,
79 const std::unique_ptr<Loss>& loss,
81 const Eigen::VectorXd& gradient,
82 const std::vector<int>& working_set,
84 const Eigen::VectorXd& x_centers,
85 const Eigen::VectorXd& x_scales,
86 const Eigen::MatrixXd& y)
89 using Eigen::MatrixXd;
90 using Eigen::VectorXd;
92 int n_working = working_set.size();
94 Eigen::VectorXd beta_old = beta(working_set);
95 Eigen::VectorXd beta_diff(n_working);
96 Eigen::VectorXd beta0_old = beta0;
98 double g_old = loss->loss(eta, y);
101 Eigen::MatrixXd residual = loss->residual(eta, y);
102 Eigen::VectorXd intercept_grad = residual.colwise().mean();
103 Eigen::VectorXd grad_working = gradient(working_set);
105 int line_search_iter = 0;
106 const int max_line_search_iter = 100;
110 beta0 = beta0_old - this->learning_rate * intercept_grad;
114 penalty.
prox(beta_old - this->learning_rate * grad_working,
115 this->learning_rate * lambda.head(n_working));
117 beta_diff = beta(working_set) - beta_old;
118 double beta_diff_norm =
119 beta_diff.dot(grad_working) +
120 (1.0 / (2 * this->learning_rate)) * beta_diff.squaredNorm();
123 Eigen::VectorXd beta0_diff = beta0 - beta0_old;
126 intercept_grad.dot(beta0_diff) +
127 (1.0 / (2 * this->learning_rate)) * beta0_diff.squaredNorm();
139 double g = loss->loss(eta, y);
140 double q = g_old + beta_diff_norm;
142 if (q >= g * (1 - 1e-12) || this->learning_rate < 1e-12) {
143 this->learning_rate *= 1.1;
146 this->learning_rate *= this->learning_rate_decr;
150 if (line_search_iter >= max_line_search_iter) {
155 if (update_type ==
"fista") {
156 if (beta_prev.size() == 0) {
160 this->t = 0.5 * (1.0 + std::sqrt(1.0 + 4.0 * t_old * t_old));
161 Eigen::VectorXd beta_current = beta(working_set);
164 (beta_current - beta_prev(working_set)) * (t_old - 1.0) / this->t;
165 beta_prev(working_set) = beta_current;
178 double learning_rate;
179 double learning_rate_decr;
180 std::string update_type;
182 Eigen::VectorXd beta_prev;
void run(Eigen::VectorXd &beta0, Eigen::VectorXd &beta, Eigen::MatrixXd &eta, const Eigen::ArrayXd &lambda, const std::unique_ptr< Loss > &loss, const SortedL1Norm &penalty, const Eigen::VectorXd &gradient, const std::vector< int > &working_set, const Eigen::MatrixXd &x, const Eigen::VectorXd &x_centers, const Eigen::VectorXd &x_scales, const Eigen::MatrixXd &y) override
Pure virtual function defining the solver's optimization routine.
Eigen::MatrixXd linearPredictor(const T &x, const std::vector< int > &active_set, const Eigen::VectorXd &beta0, const Eigen::VectorXd &beta, const Eigen::VectorXd &x_centers, const Eigen::VectorXd &x_scales, const JitNormalization jit_normalization, const bool intercept)