slope.h

Go to the documentation of this file.

 
#pragma once
 
#include "logger.h"
#include "screening.h"
#include "slope_fit.h"
#include "slope_path.h"
#include "solvers/hybrid_cd.h"
#include "timer.h"
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <cassert>
#include <optional>
 
namespace slope {
 
class Slope
{
public:
  Slope() = default;
 
  void setSolver(const std::string& solver);
 
  void setIntercept(bool intercept);
 
  void setNormalization(const std::string& type);
 
  void setUpdateClusters(bool update_clusters);
 
  void setReturnClusters(const bool return_clusters);
 
  void setAlphaMinRatio(double alpha_min_ratio);
 
  void setAlphaType(const std::string& alpha_type);
 
  void setLearningRateDecr(double learning_rate_decr);
 
  void setQ(double q);
 
  void setOscarParameters(const double theta1, const double theta2);
 
  void setTol(double tol);
 
  void setRelaxTol(double tol);
 
  void setRelaxMaxOuterIterations(int max_it);
 
  void setRelaxMaxInnerIterations(int max_it);
 
  void setMaxIterations(int max_it);
 
  void setPathLength(int path_length);
 
  void setHybridCdIterations(int cd_iterations);
 
  void setLambdaType(const std::string& lambda_type);
 
  void setLoss(const std::string& loss_type);
 
  void setScreening(const std::string& screening_type);
 
  void setModifyX(const bool modify_x);
 
  void setDevChangeTol(const double dev_change_tol);
 
  void setDevRatioTol(const double dev_ratio_tol);
 
  void setMaxClusters(const int max_clusters);
 
  void setCentering(const std::string& type);
 
  void setCentering(const Eigen::VectorXd& x_centers);
 
  void setScaling(const std::string& type);
 
  void setDiagnostics(const bool collect_diagnostics);
 
  void setScaling(const Eigen::VectorXd& x_scales);
 
  void setAlphaEstimationMaxIterations(const int alpha_est_maxit);
 
  int getAlphaEstimationMaxIterations() const;
 
  bool getFitIntercept() const;
 
  const std::string& getLossType();
 
  template<typename T>
  SlopePath path(T& x,
                 const Eigen::MatrixXd& y_in,
                 Eigen::ArrayXd alpha = Eigen::ArrayXd::Zero(0),
                 Eigen::ArrayXd lambda = Eigen::ArrayXd::Zero(0));
 
  template<typename T>
  SlopeFit fit(T& x,
               const Eigen::MatrixXd& y_in,
               const double alpha = 1.0,
               Eigen::ArrayXd lambda = Eigen::ArrayXd::Zero(0));
 
  template<typename T>
  SlopeFit relax(const SlopeFit& fit,
                 T& x,
                 const Eigen::VectorXd& y_in,
                 const double gamma = 0.0,
                 Eigen::VectorXd beta0 = Eigen::VectorXd(0),
                 Eigen::VectorXd beta = Eigen::VectorXd(0))
  {
    using Eigen::MatrixXd;
    using Eigen::VectorXd;
 
    int n = x.rows();
    int p = x.cols();
 
    if (beta0.size() == 0) {
      beta0 = fit.getIntercepts(false);
    }
 
    if (beta.size() == 0) {
      beta = fit.getCoefs(false);
    }
 
    double alpha = 0;
 
    Timer timer;
 
    std::vector<double> primals, duals, time;
    timer.start();
 
    auto jit_normalization =
      normalize(x, x_centers, x_scales, centering_type, scaling_type, modify_x);
 
    bool update_clusters = false;
 
    std::unique_ptr<Loss> loss = setupLoss(this->loss_type);
 
    MatrixXd y = loss->preprocessResponse(y_in);
 
    int m = y.cols();
 
    Eigen::ArrayXd lambda_relax = Eigen::ArrayXd::Zero(p * m);
 
    auto working_set = activeSet(beta);
 
    Eigen::MatrixXd eta = linearPredictor(x,
                                          working_set,
                                          beta0,
                                          beta,
                                          x_centers,
                                          x_scales,
                                          jit_normalization,
                                          intercept);
    VectorXd gradient = VectorXd::Zero(p * m);
    VectorXd residual(n);
    VectorXd working_residual(n);
 
    VectorXd w = VectorXd::Ones(n);
    VectorXd w_ones = VectorXd::Ones(n);
    VectorXd z = y;
 
    Clusters clusters = fit.getClusters();
 
    int passes = 0;
 
    for (int irls_it = 0; irls_it < max_it_outer_relax; irls_it++) {
      residual = loss->residual(eta, y);
 
      if (collect_diagnostics) {
        primals.push_back(loss->loss(eta, y));
        duals.push_back(0.0);
        time.push_back(timer.elapsed());
      }
 
      Eigen::VectorXd cluster_gradient = clusterGradient(beta,
                                                         residual,
                                                         clusters,
                                                         x,
                                                         w_ones,
                                                         x_centers,
                                                         x_scales,
                                                         jit_normalization);
 
      double norm_grad = cluster_gradient.lpNorm<Eigen::Infinity>();
 
      if (norm_grad < tol_relax) {
        break;
      }
 
      loss->updateWeightsAndWorkingResponse(w, z, eta, y);
      working_residual = eta - z;
 
      for (int inner_it = 0; inner_it < max_it_inner_relax; ++inner_it) {
        passes++;
 
        double max_abs_gradient = coordinateDescent(beta0,
                                                    beta,
                                                    working_residual,
                                                    clusters,
                                                    lambda_relax,
                                                    x,
                                                    w,
                                                    x_centers,
                                                    x_scales,
                                                    intercept,
                                                    jit_normalization,
                                                    update_clusters);
 
        if (max_abs_gradient < tol_relax) {
          break;
        }
      }
 
      eta = working_residual + z;
 
      if (irls_it == max_it_outer_relax) {
        WarningLogger::addWarning(WarningCode::MAXIT_REACHED,
                                  "Maximum number of IRLS iterations reached.");
      }
    }
 
    double dev = loss->deviance(eta, y);
 
    if (gamma > 0) {
      Eigen::VectorXd old_coefs = fit.getCoefs(false);
      Eigen::VectorXd old_intercept = fit.getIntercepts(false);
      beta = (1 - gamma) * beta + gamma * old_coefs;
    }
 
    SlopeFit fit_out{ beta0,
                      beta.reshaped(p, m).sparseView(),
                      clusters,
                      alpha,
                      lambda_relax,
                      dev,
                      fit.getNullDeviance(),
                      primals,
                      duals,
                      time,
                      passes,
                      centering_type,
                      scaling_type,
                      intercept,
                      x_centers,
                      x_scales };
 
    return fit_out;
  }
 
  template<typename T>
  SlopePath relax(const SlopePath& path,
                  T& x,
                  const Eigen::VectorXd& y,
                  const double gamma = 0.0)
  {
    std::vector<SlopeFit> fits;
 
    Eigen::VectorXd beta0 = path(0).getIntercepts(false);
    Eigen::VectorXd beta = path(0).getCoefs(false);
 
    for (size_t i = 0; i < path.size(); i++) {
      auto relaxed_fit = relax(path(i), x, y, gamma, beta0, beta);
 
      fits.emplace_back(relaxed_fit);
 
      // Update warm starts
      // TODO: Maybe be more clever about whether to use the
      // previous values or the regularized estimates and warm starts.
      // Maybe just pick the solution with larger coefficients?
      beta0 = relaxed_fit.getIntercepts(false);
      beta = relaxed_fit.getCoefs(false);
    }
 
    return fits;
  }
 
private:
  // Parameters
  bool collect_diagnostics = false;
  bool intercept = true;
  bool modify_x = false;
  bool return_clusters = true;
  bool update_clusters = false;
  double alpha_min_ratio = -1;
  double dev_change_tol = 1e-5;
  double dev_ratio_tol = 0.999;
  double learning_rate_decr = 0.5;
  double q = 0.1;
  double theta1 = 1.0;
  double theta2 = 0.5;
  double tol = 1e-4;
  double tol_relax = 1e-4;
  int alpha_est_maxit = 1000;
  int cd_iterations = 10;
  int max_it = 1e4;
  int max_it_inner_relax = 1e5;
  int max_it_outer_relax = 50;
  int path_length = 100;
  std::optional<int> max_clusters = std::nullopt;
  std::string alpha_type = "path";
  std::string centering_type = "mean";
  std::string lambda_type = "bh";
  std::string loss_type = "quadratic";
  std::string scaling_type = "sd";
  std::string screening_type = "strong";
  std::string solver_type = "auto";
 
  // Data
  Eigen::VectorXd x_centers;
  Eigen::VectorXd x_scales;
};
 
} // namespace slope