math.h

Go to the documentation of this file.

 
#pragma once
 
#include "clusters.h"
#include "eigen_compat.h"
#include "jit_normalization.h"
#include "utils.h"
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <algorithm>
#include <cassert>
#include <limits>
#include <numeric>
#include <vector>
 
#ifdef _OPENMP
#include "threads.h"
#endif
 
namespace slope {

template<typename T>
int
sign(T val)
{
  return (T(0) < val) - (val < T(0));
}

template<typename T>
Eigen::ArrayXd
cumSum(const T& x, const bool leading_zero = false)
{
  const size_t start = leading_zero ? 1 : 0;
  Eigen::ArrayXd out(x.size() + start);
 
  if (leading_zero) {
    out(0) = 0.0;
  }
 
  std::partial_sum(x.begin(), x.end(), out.begin() + start);
 
  return out;
}

template<typename T>
T
sigmoid(const T& x)
{
  return 1.0 / (1.0 + std::exp(-x));
}

template<typename T>
T
logit(const T& x)
{
  assert(x > 0 && x < 1 && "Input must be in (0, 1)");
 
  return std::log(x) - std::log1p(-x);
}

template<typename T>
T
clamp(const T& x, const T& lo, const T& hi)
{
  return x < lo ? lo : x > hi ? hi : x;
}

Eigen::VectorXd
logSumExp(const Eigen::MatrixXd& a);

Eigen::MatrixXd
softmax(const Eigen::MatrixXd& x);
 
/*
 * Computes the gradient of the loss with respect to \f(\beta\f).
 *
 * @tparam T The type of the input matrix.
 * @param x The input matrix.
 * @param active_set Indicies for active set
 * @param beta0 Intercept
 * @param beta Coefficients
 * @param x_centers The vector of center values for each column of x.
 * @param x_scales The vector of scale values for each column of x.
 * @param jit_normalization Type of JIT normalization.
 * @param intercept Whether to fit an intercept.
 * @return The computed gradient vector.
 */
template<typename T>
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)
{
  int n = x.rows();
  int p = x.cols();
  int m = beta0.size();
 
  Eigen::MatrixXd eta = Eigen::MatrixXd::Zero(n, m);
 
#ifdef _OPENMP
  bool large_problem = active_set.size() > 100 && n * active_set.size() > 1e7;
#pragma omp parallel num_threads(Threads::get()) if (large_problem)
#endif
  {
    Eigen::MatrixXd eta_local = Eigen::MatrixXd::Zero(n, m);
 
#ifdef _OPENMP
#pragma omp for nowait
#endif
    for (int i = 0; i < static_cast<int>(active_set.size()); ++i) {
      int ind = active_set[i];
      auto [k, j] = std::div(ind, p);
 
      switch (jit_normalization) {
        case JitNormalization::Both:
          eta_local.col(k) += x.col(j) * beta(ind) / x_scales(j);
          eta_local.col(k).array() -= beta(ind) * x_centers(j) / x_scales(j);
          break;
 
        case JitNormalization::Center:
          eta_local.col(k) += x.col(j) * beta(ind);
          eta_local.col(k).array() -= beta(ind) * x_centers(j);
          break;
 
        case JitNormalization::Scale:
          eta_local.col(k) += x.col(j) * beta(ind) / x_scales(j);
          break;
 
        case JitNormalization::None:
          eta_local.col(k) += x.col(j) * beta(ind);
          break;
      }
    }
 
#ifdef _OPENMP
#pragma omp critical
#endif
    {
      eta += eta_local;
    }
  }
 
  if (intercept) {
    eta.rowwise() += beta0.transpose();
  }
 
  return eta;
}

namespace detail {
 
template<typename T, typename IndexAt>
void
updateGradientImpl(Eigen::VectorXd& gradient,
                   const T& x,
                   const Eigen::MatrixXd& residual,
                   const int active_size,
                   const IndexAt& index_at,
                   const Eigen::VectorXd& x_centers,
                   const Eigen::VectorXd& x_scales,
                   const Eigen::VectorXd& w,
                   const JitNormalization jit_normalization)
{
  const int n = x.rows();
  const int p = x.cols();
  const int m = residual.cols();
 
  assert(gradient.size() == p * m &&
         "Gradient matrix has incorrect dimensions");
 
  Eigen::MatrixXd weighted_residual(n, m);
  Eigen::ArrayXd wr_sums(m);
 
#ifdef _OPENMP
  bool large_problem =
    active_size > 100 && static_cast<long long>(n) * active_size > 1e5;
#pragma omp parallel for num_threads(Threads::get()) if (large_problem)
#endif
  for (int k = 0; k < m; ++k) {
    weighted_residual.col(k) = residual.col(k).cwiseProduct(w);
    for (int i = 0; i < n; ++i) {
      if (w(i) == 0.0) {
        weighted_residual(i, k) = 0.0;
      }
    }
    wr_sums(k) = weighted_residual.col(k).sum();
  }
 
#ifdef _OPENMP
#pragma omp parallel for num_threads(Threads::get()) if (large_problem)
#endif
  for (int i = 0; i < active_size; ++i) {
    int ind = index_at(i);
    auto [k, j] = std::div(ind, p);
 
    switch (jit_normalization) {
      case JitNormalization::Both:
        gradient(ind) =
          (x.col(j).dot(weighted_residual.col(k)) - x_centers(j) * wr_sums(k)) /
          (x_scales(j) * n);
        break;
      case JitNormalization::Center:
        gradient(ind) =
          (x.col(j).dot(weighted_residual.col(k)) - x_centers(j) * wr_sums(k)) /
          n;
        break;
      case JitNormalization::Scale:
        gradient(ind) =
          x.col(j).dot(weighted_residual.col(k)) / (x_scales(j) * n);
        break;
      case JitNormalization::None:
        gradient(ind) = x.col(j).dot(weighted_residual.col(k)) / n;
        break;
    }
  }
}
 
} // namespace detail

template<typename T>
void
updateGradient(Eigen::VectorXd& gradient,
               const T& x,
               const Eigen::MatrixXd& residual,
               const std::vector<int>& active_set,
               const Eigen::VectorXd& x_centers,
               const Eigen::VectorXd& x_scales,
               const Eigen::VectorXd& w,
               const JitNormalization jit_normalization)
{
  detail::updateGradientImpl(
    gradient,
    x,
    residual,
    static_cast<int>(active_set.size()),
    [&active_set](int i) { return active_set[i]; },
    x_centers,
    x_scales,
    w,
    jit_normalization);
}

template<typename T>
void
updateGradient(Eigen::VectorXd& gradient,
               const T& x,
               const Eigen::MatrixXd& residual,
               const Eigen::VectorXd& x_centers,
               const Eigen::VectorXd& x_scales,
               const Eigen::VectorXd& w,
               const JitNormalization jit_normalization)
{
  detail::updateGradientImpl(
    gradient,
    x,
    residual,
    static_cast<int>(gradient.size()),
    [](int i) { return i; },
    x_centers,
    x_scales,
    w,
    jit_normalization);
}
 
/*
 * Computes the gradient of the loss with respect to \f(\beta\f).
 *
 * @tparam T The type of the input matrix.
 * @param gradient The residual vector.
 * @param x The input matrix.
 * @param offset Gradient offset
 * @param active_set Indices for the active_set
 * @param x_centers The vector of center values for each column of x.
 * @param x_scales The vector of scale values for each column of x.
 * @param jit_normalization Type of JIT normalization
 * just-in-time.
 */
namespace detail {
 
template<typename T, typename IndexAt>
void
offsetGradientImpl(Eigen::VectorXd& gradient,
                   const T& x,
                   const Eigen::VectorXd& offset,
                   const int active_size,
                   const IndexAt& index_at,
                   const Eigen::VectorXd& x_centers,
                   const Eigen::VectorXd& x_scales,
                   const JitNormalization jit_normalization)
{
  const int n = x.rows();
  const int p = x.cols();
 
  for (int i = 0; i < active_size; ++i) {
    int ind = index_at(i);
    auto [k, j] = std::div(ind, p);
 
    switch (jit_normalization) {
      case JitNormalization::Both:
        gradient(ind) -=
          offset(k) * (x.col(j).sum() / n - x_centers(j)) / x_scales(j);
        break;
      case JitNormalization::Center:
        gradient(ind) -= offset(k) * (x.col(j).sum() / n - x_centers(j));
        break;
      case JitNormalization::Scale:
        gradient(ind) -= offset(k) * x.col(j).sum() / (n * x_scales(j));
        break;
      case JitNormalization::None:
        gradient(ind) -= offset(k) * x.col(j).sum() / n;
        break;
    }
  }
}
 
} // namespace detail

template<typename T>
void
offsetGradient(Eigen::VectorXd& gradient,
               const T& x,
               const Eigen::VectorXd& offset,
               const std::vector<int>& active_set,
               const Eigen::VectorXd& x_centers,
               const Eigen::VectorXd& x_scales,
               const JitNormalization jit_normalization)
{
  detail::offsetGradientImpl(
    gradient,
    x,
    offset,
    static_cast<int>(active_set.size()),
    [&active_set](int i) { return active_set[i]; },
    x_centers,
    x_scales,
    jit_normalization);
}

template<typename T>
void
offsetGradient(Eigen::VectorXd& gradient,
               const T& x,
               const Eigen::VectorXd& offset,
               const Eigen::VectorXd& x_centers,
               const Eigen::VectorXd& x_scales,
               const JitNormalization jit_normalization)
{
  detail::offsetGradientImpl(
    gradient,
    x,
    offset,
    static_cast<int>(gradient.size()),
    [](int i) { return i; },
    x_centers,
    x_scales,
    jit_normalization);
}

std::vector<int>
setUnion(const std::vector<int>& a, const std::vector<int>& b);

std::vector<int>
setDiff(const std::vector<int>& a, const std::vector<int>& b);

template<typename T>
int
whichMax(const T& x)
{
  return std::distance(x.begin(), std::max_element(x.begin(), x.end()));
}

template<typename T>
int
whichMin(const T& x)
{
  return std::distance(x.begin(), std::min_element(x.begin(), x.end()));
}

template<typename T, typename Comparator>
int
whichBest(const T& x, const Comparator& comp)
{
  if (x.size() == 0)
    return -1;
 
  return std::distance(x.begin(), std::max_element(x.begin(), x.end(), comp));
}

Eigen::ArrayXd
geomSpace(const double start, const double end, const int n);

template<typename T>
Eigen::VectorXd
l1Norms(const T& x)
{
  const int p = x.cols();
 
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    out(j) = x.col(j).cwiseAbs().sum();
  }
 
  return out;
}

template<typename T>
Eigen::VectorXd
l2Norms(const Eigen::SparseMatrixBase<T>& x)
{
  const int p = x.cols();
 
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    out(j) = x.col(j).norm();
  }
 
  return out;
}

template<typename T>
Eigen::VectorXd
l2Norms(const Eigen::MatrixBase<T>& x)
{
  return x.colwise().norm();
}

template<typename T>
Eigen::VectorXd
maxAbs(const Eigen::SparseMatrixBase<T>& x)
{
  const int p = x.cols();
 
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    double x_j_maxabs = 0.0;
 
    for (typename T::InnerIterator it(x.derived(), j); it; ++it) {
      x_j_maxabs = std::max(x_j_maxabs, std::abs(it.value()));
    }
 
    out(j) = x_j_maxabs;
  }
 
  return out;
}

template<typename T>
Eigen::VectorXd
maxAbs(const Eigen::MatrixBase<T>& x)
{
  return x.cwiseAbs().colwise().maxCoeff();
}

template<typename T>
Eigen::VectorXd
means(const Eigen::SparseMatrixBase<T>& x)
{
  const int n = x.rows();
  const int p = x.cols();
 
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    out(j) = x.col(j).sum() / n;
  }
 
  return out;
}

template<typename T>
Eigen::VectorXd
means(const Eigen::MatrixBase<T>& x)
{
  return x.colwise().mean();
}

template<typename T>
Eigen::VectorXd
stdDevs(const Eigen::SparseMatrixBase<T>& x)
{
  const int n = x.rows();
  const int p = x.cols();
 
  Eigen::VectorXd x_means = means(x);
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    double sum_sq_diff = 0.0;
    const double mean = x_means(j);
    double minimum = std::numeric_limits<double>::infinity();
    double maximum = -std::numeric_limits<double>::infinity();
 
    // Process non-zero elements
    for (typename T::InnerIterator it(x.derived(), j); it; ++it) {
      double diff = it.value() - mean;
      sum_sq_diff += diff * diff;
      minimum = std::min(minimum, it.value());
      maximum = std::max(maximum, it.value());
    }
 
    // Account for zeros
    int nz_count = slope::nonZeros(x.col(j));
    if (nz_count < n) {
      sum_sq_diff += (n - nz_count) * mean * mean;
      minimum = std::min(minimum, 0.0);
      maximum = std::max(maximum, 0.0);
    }
 
    // Rounding in the mean must not give constant columns a nonzero scale.
    if (minimum == maximum) {
      out(j) = 0.0;
      continue;
    }
 
    // Standard deviation is sqrt of the average of squared differences
    out(j) = std::sqrt(sum_sq_diff / n);
  }
 
  return out;
}

template<typename T>
Eigen::VectorXd
stdDevs(const Eigen::MatrixBase<T>& x)
{
  int n = x.rows();
  int p = x.cols();
 
  Eigen::VectorXd x_means = means(x);
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    // Rounding in the mean must not give constant columns a nonzero scale.
    if (n > 0 && x.col(j).minCoeff() == x.col(j).maxCoeff()) {
      out(j) = 0.0;
    } else {
      out(j) = (x.col(j).array() - x_means(j)).matrix().norm();
    }
  }
 
  out.array() /= std::sqrt(n);
 
  return out;
}

template<typename T>
Eigen::VectorXd
ranges(const Eigen::SparseMatrixBase<T>& x)
{
  const int p = x.cols();
 
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    double x_j_max = 0.0;
    double x_j_min = 0.0;
 
    for (typename T::InnerIterator it(x.derived(), j); it; ++it) {
      x_j_max = std::max(x_j_max, it.value());
      x_j_min = std::min(x_j_min, it.value());
    }
 
    out(j) = x_j_max - x_j_min;
  }
 
  return out;
}

template<typename T>
Eigen::VectorXd
ranges(const Eigen::MatrixBase<T>& x)
{
  return x.colwise().maxCoeff() - x.colwise().minCoeff();
}

template<typename T>
Eigen::VectorXd
mins(const Eigen::SparseMatrixBase<T>& x)
{
  const int p = x.cols();
 
  Eigen::VectorXd out(p);
 
  for (int j = 0; j < p; ++j) {
    double x_j_min = 0.0;
 
    for (typename T::InnerIterator it(x.derived(), j); it; ++it) {
      x_j_min = std::min(x_j_min, it.value());
    }
 
    out(j) = x_j_min;
  }
 
  return out;
}

template<typename T>
Eigen::VectorXd
mins(const Eigen::MatrixBase<T>& x)
{
  return x.colwise().minCoeff();
}

template<typename T>
Eigen::VectorXd
clusterGradient(Eigen::VectorXd& beta,
                Eigen::MatrixXd& residual,
                Clusters& clusters,
                const T& x,
                const Eigen::MatrixXd& w,
                const Eigen::VectorXd& x_centers,
                const Eigen::VectorXd& x_scales,
                const JitNormalization jit_normalization)
{
  using namespace Eigen;
 
  const int n = x.rows();
  const int p = x.cols();
  const int n_clusters = clusters.size();
 
  Eigen::VectorXd gradient = Eigen::VectorXd::Zero(n_clusters);
 
  for (int j = 0; j < n_clusters; ++j) {
    double c_old = clusters.coeff(j);
 
    if (c_old == 0) {
      gradient(j) = 0;
      continue;
    }
 
    int cluster_size = clusters.cluster_size(j);
    std::vector<int> s;
    s.reserve(cluster_size);
 
    for (auto c_it = clusters.cbegin(j); c_it != clusters.cend(j); ++c_it) {
      int ind = *c_it;
      double s_k = sign(beta(ind));
      s.emplace_back(s_k);
    }
 
    double hess = 1;
    double grad = 0;
 
    if (cluster_size == 1) {
      int k = *clusters.cbegin(j);
      std::tie(grad, hess) = computeGradientAndHessian(
        x, k, w, residual, x_centers, x_scales, s[0], jit_normalization, n);
    } else {
      std::tie(hess, grad) = computeClusterGradientAndHessian(
        x, j, s, clusters, w, residual, x_centers, x_scales, jit_normalization);
    }
 
    gradient(j) = grad;
  }
 
  return gradient;
}
 
} // namespace slope