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>
#include <optional>
 
namespace slope {
 
class Hybrid : public SolverBase
{
public:
  Hybrid(JitNormalization jit_normalization,
         bool intercept,
         bool update_clusters,
         int cd_iterations,
         const std::string& cd_type,
         std::optional<int> random_seed = std::nullopt)
    : SolverBase(jit_normalization, intercept)
    , update_clusters(update_clusters)
    , cd_iterations(cd_iterations)
    , cd_type(cd_type)
    , rng(random_seed.has_value() ? std::mt19937(*random_seed)
                                  : std::mt19937(std::random_device{}()))
  {
  }
 
  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;
 
  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::Map<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::Map<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::MatrixXd& y)
  {
    using Eigen::MatrixXd;
    using Eigen::VectorXd;
 
    const int n = x.rows();
    const int m = eta.cols();
 
    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);
 
    // TODO: Make these parameters and initialize once
    MatrixXd w = MatrixXd::Ones(n, m);
    MatrixXd z = y;
 
    loss->updateWeightsAndWorkingResponse(w, z, eta, y);
 
    MatrixXd residual = eta - z;
    const VectorXd weight_sums = w.colwise().sum().transpose();
 
    Eigen::ArrayXd lambda_cumsum(lambda.size() + 1);
    lambda_cumsum(0) = 0.0;
    std::partial_sum(lambda.begin(), lambda.end(), lambda_cumsum.begin() + 1);
 
    for (int it = 0; it < this->cd_iterations; ++it) {
      double old_obj =
        computeObjective(penalty, beta, residual, w, lambda, working_set);
 
      // Store old values to revert if no progress is made
      Clusters old_clusters = clusters;
      Eigen::MatrixXd old_residual = residual;
      Eigen::VectorXd old_beta = beta;
      Eigen::VectorXd old_beta0 = beta0;
 
      coordinateDescent(beta0,
                        beta,
                        residual,
                        clusters,
                        lambda_cumsum,
                        x,
                        w,
                        weight_sums,
                        x_centers,
                        x_scales,
                        this->intercept,
                        this->jit_normalization,
                        this->update_clusters,
                        rng,
                        this->cd_type);
 
      double new_obj =
        computeObjective(penalty, beta, residual, w, lambda, working_set);
 
      if (!std::isfinite(new_obj) || new_obj > old_obj) {
        // No progress, revert to previous state
        clusters = old_clusters;
        residual = old_residual;
        beta = old_beta;
        beta0 = old_beta0;
 
        break;
      }
    }
 
    // 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 computeObjective(const SortedL1Norm& penalty,
                          const Eigen::VectorXd& beta,
                          const Eigen::MatrixXd& residual,
                          const Eigen::MatrixXd& w,
                          const Eigen::ArrayXd& lambda,
                          const std::vector<int>& working_set)
  {
    double val =
      0.5 * (residual.array().square() * w.array()).sum() / residual.rows() +
      penalty.eval(beta(working_set), lambda.head(working_set.size()));
 
    return val;
  }
 
  // TODO: These should be used in the PGD solver and taken as arguments to the
  // Hybrid solver and not just set and ignored here.
  double pgd_learning_rate =
    1.0; 
  double pgd_learning_rate_decr =
    0.5; 
 
  bool update_clusters = false; 
  int cd_iterations = 10;       
  std::string cd_type =
    "cyclical"; 
  std::mt19937 rng{
    std::random_device{}()
  }; 
};
 
} // namespace slope