loss.h

Go to the documentation of this file.

 
#pragma once
 
#include <Eigen/Core>
 
namespace slope {

class Loss
{
public:
  virtual ~Loss() = default;

  virtual double loss(const Eigen::MatrixXd& eta, const Eigen::MatrixXd& y) = 0;

  virtual double loss(const Eigen::MatrixXd& eta,
                      const Eigen::MatrixXd& y,
                      const Eigen::VectorXd& observation_weights);

  virtual double dual(const Eigen::MatrixXd& theta,
                      const Eigen::MatrixXd& y,
                      const Eigen::VectorXd& w) = 0;

  Eigen::MatrixXd residual(const Eigen::MatrixXd& eta,
                           const Eigen::MatrixXd& y);
  // TODO: This function could be removed since it is just
  // g^{-1}(eta) - y for all families

  virtual Eigen::MatrixXd dualPoint(const Eigen::MatrixXd& eta,
                                    const Eigen::MatrixXd& y,
                                    bool fit_intercept);

  virtual Eigen::MatrixXd dualPoint(const Eigen::MatrixXd& eta,
                                    const Eigen::MatrixXd& y,
                                    bool fit_intercept,
                                    const Eigen::VectorXd& observation_weights);

  virtual Eigen::MatrixXd hessianDiagonal(const Eigen::MatrixXd& eta) = 0;

  virtual void updateWeightsAndWorkingResponse(Eigen::MatrixXd& w,
                                               Eigen::MatrixXd& z,
                                               const Eigen::MatrixXd& eta,
                                               const Eigen::MatrixXd& y);

  virtual void updateWeightsAndWorkingResponse(
    Eigen::MatrixXd& w,
    Eigen::MatrixXd& z,
    const Eigen::MatrixXd& eta,
    const Eigen::MatrixXd& y,
    const Eigen::VectorXd& observation_weights);

  virtual Eigen::MatrixXd preprocessResponse(const Eigen::MatrixXd& y) = 0;

  virtual void updateIntercept(Eigen::VectorXd& beta0,
                               const Eigen::MatrixXd& eta,
                               const Eigen::MatrixXd& y)
  {
    Eigen::MatrixXd residual = this->residual(eta, y);
    beta0 -= residual.colwise().mean() / this->lipschitz_constant;
  };
  // TODO: This function is not currently used anywhere and could
  // probably be removed.

  virtual Eigen::MatrixXd link(const Eigen::MatrixXd& mu) = 0;

  virtual Eigen::MatrixXd inverseLink(const Eigen::MatrixXd& eta) = 0;

  virtual Eigen::MatrixXd predict(const Eigen::MatrixXd& eta) = 0;

  virtual double deviance(const Eigen::MatrixXd& eta, const Eigen::MatrixXd& y)
  {
    return 2.0 * (loss(eta, y) - loss(link(y), y));
  }

  virtual double deviance(const Eigen::MatrixXd& eta,
                          const Eigen::MatrixXd& y,
                          const Eigen::VectorXd& observation_weights);
 
protected:
  Eigen::MatrixXd residualOnSupport(const Eigen::MatrixXd& eta,
                                    const Eigen::MatrixXd& y,
                                    const Eigen::VectorXd& observation_weights);

  explicit Loss(double lipschitz_constant)
    : lipschitz_constant(lipschitz_constant)
  {
  }

  Eigen::MatrixXd domainSafePoint(const Eigen::MatrixXd& mean,
                                  const Eigen::MatrixXd& y,
                                  double lower,
                                  double upper) const;
 
  Eigen::MatrixXd domainSafePoint(
    const Eigen::MatrixXd& mean,
    const Eigen::MatrixXd& y,
    double lower,
    double upper,
    const Eigen::VectorXd& observation_weights) const;
 
private:
  const double lipschitz_constant;
};
 
} // namespace slope