hybrid.h

Go to the documentation of this file.

 
#pragma once
 
#include "../clusters.h"
#include "../losses/loss.h"
#include "../sorted_l1_norm.h"
#include "hybrid_cd.h"
#include "pgd.h"
#include "solver.h"
#include <memory>
 
namespace slope {
 
class Hybrid : public SolverBase
{
public:
  Hybrid(JitNormalization jit_normalization,
         bool intercept,
         bool update_clusters,
         int cd_iterations)
    : SolverBase(jit_normalization, intercept)
    , update_clusters(update_clusters)
    , cd_iterations(cd_iterations)
  {
  }
 
  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;
 
  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::SparseMatrix<double>& x,
           const Eigen::VectorXd& x_centers,
           const Eigen::VectorXd& x_scales,
           const Eigen::MatrixXd& y) override;
 
private:
  template<typename MatrixType>
  void runImpl(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_in,
               const std::vector<int>& working_set,
               const MatrixType& x,
               const Eigen::VectorXd& x_centers,
               const Eigen::VectorXd& x_scales,
               const Eigen::VectorXd& y)
  {
    using Eigen::MatrixXd;
    using Eigen::VectorXd;
 
    const int n = x.rows();
 
    PGD pgd_solver(jit_normalization, intercept, "pgd");
 
    // Run proximal gradient descent
    pgd_solver.run(beta0,
                   beta,
                   eta,
                   lambda,
                   loss,
                   penalty,
                   gradient_in,
                   working_set,
                   x,
                   x_centers,
                   x_scales,
                   y);
 
    Clusters clusters(beta);
 
    VectorXd w = VectorXd::Ones(n);
    VectorXd z = y;
    loss->updateWeightsAndWorkingResponse(w, z, eta, y);
 
    VectorXd residual = eta - z;
 
    for (int it = 0; it < this->cd_iterations; ++it) {
      coordinateDescent(beta0,
                        beta,
                        residual,
                        clusters,
                        lambda,
                        x,
                        w,
                        x_centers,
                        x_scales,
                        this->intercept,
                        this->jit_normalization,
                        this->update_clusters);
    }
 
    // The residual is kept up to date, but not eta. So we need to compute
    // it here.
    eta = residual + z;
    // TODO: register convergence status
  }
 
  double pgd_learning_rate =
    1.0; 
  double pgd_learning_rate_decr =
    0.5; 
  bool update_clusters = false; 
  int cd_iterations = 10;       
};
 
} // namespace slope