wighted least square algorithm (wls.cpp)

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(RcppEigen)]]
#include <RcppArmadillo.h>
#include <RcppEigen.h>
using namespace arma;

// [[Rcpp::export]]
Eigen::VectorXd eigen_fullPivHHQR(const Eigen::Map<Eigen::MatrixXd> & X,
                                  const Eigen::Map<Eigen::VectorXd> & w,
                                  const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::VectorXd beta = (X.transpose() * w.asDiagonal() * X).fullPivHouseholderQr().solve(X.transpose() * w.asDiagonal() * y);;
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_colPivHHQR(const Eigen::Map<Eigen::MatrixXd> & X,
                                 const Eigen::Map<Eigen::VectorXd> & w,
                                 const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::VectorXd beta = (X.transpose() * w.asDiagonal() * X).colPivHouseholderQr().solve(X.transpose() * w.asDiagonal() * y);;
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_HHQR(const Eigen::Map<Eigen::MatrixXd> & X,
                           const Eigen::Map<Eigen::VectorXd> & w,
                           const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::VectorXd beta = (X.transpose() * w.asDiagonal() * X).householderQr().solve(X.transpose() * w.asDiagonal() * y);
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_fullLU(const Eigen::Map<Eigen::MatrixXd> & X,
                             const Eigen::Map<Eigen::VectorXd> & w,
                             const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::VectorXd beta = (X.transpose() * w.asDiagonal() * X).fullPivLu().solve(X.transpose() * w.asDiagonal() * y);
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_llt(const Eigen::Map<Eigen::MatrixXd> & X,
                          const Eigen::Map<Eigen::VectorXd> & w,
                          const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::VectorXd beta = (X.transpose() * w.asDiagonal() * X).llt().solve(X.transpose() * w.asDiagonal() * y);
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_ldlt(const Eigen::Map<Eigen::MatrixXd> & X,
                           const Eigen::Map<Eigen::VectorXd> & w,
                           const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::VectorXd beta = (X.transpose() * w.asDiagonal() * X).ldlt().solve(X.transpose() * w.asDiagonal() * y);
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_chol_llt1(const Eigen::Map<Eigen::MatrixXd> & X,
                                const Eigen::Map<Eigen::VectorXd> & w,
                                const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::MatrixXd XWX = X.transpose() * w.asDiagonal() * X;
  Eigen::MatrixXd R = XWX.llt().matrixU();
  Eigen::VectorXd XWY = X.transpose() * (w.array() * y.array()).matrix();
  Eigen::VectorXd beta = R.householderQr().solve(R.transpose().householderQr().solve(XWY));
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_chol_llt2(const Eigen::Map<Eigen::MatrixXd> & X,
                                const Eigen::Map<Eigen::VectorXd> & w,
                                const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::MatrixXd XW = X.transpose() * w.asDiagonal();
  Eigen::MatrixXd R = (XW * X).llt().matrixU();
  Eigen::VectorXd beta = R.householderQr().solve(R.transpose().householderQr().solve(XW * y));
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_chol_llt3(const Eigen::Map<Eigen::MatrixXd> & X,
                                const Eigen::Map<Eigen::VectorXd> & w,
                                const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::MatrixXd XW(X.cols(), X.rows());
  for (unsigned int i = 0; i < X.cols(); ++i)
    XW.row(i) = X.col(i).array() * w.array();
  Eigen::MatrixXd R = (XW * X).llt().matrixU();
  Eigen::VectorXd beta = R.householderQr().solve(R.transpose().householderQr().solve(XW * y));
  return beta;
}

// [[Rcpp::export]]
Eigen::VectorXd eigen_colPivHHQR2(const Eigen::Map<Eigen::MatrixXd> & X,
                                  const Eigen::Map<Eigen::VectorXd> & w,
                                  const Eigen::Map<Eigen::VectorXd> & y) {
  Eigen::VectorXd sw = w.cwiseSqrt();
  Eigen::MatrixXd XW(X.rows(), X.cols());
  for (unsigned int i = 0; i < X.cols(); ++i)
    XW.col(i) = X.col(i).array() * sw.array();
  Eigen::VectorXd beta = XW.colPivHouseholderQr().solve(y.cwiseProduct(sw));
  return beta;
}

// [[Rcpp::export]]
arma::vec arma_qr(const arma::mat& X, const arma::vec& w, const arma::vec& y) {
  mat Q, R;
  qr(Q, R, X.each_col() % sqrt(w));
  vec p = solve(R.head_rows(X.n_cols), Q.head_cols(X.n_cols).t() * (y % sqrt(w)));
  return p;
}

// [[Rcpp::export]]
arma::vec arma_pinv(const arma::mat& X, const arma::vec& w, const arma::vec& y) {
  vec p = pinv(X.t() * (repmat(w, 1, X.n_cols) % X)) * X.t() * (w % y);
  return p;
}

// [[Rcpp::export]]
arma::vec arma_chol1(const arma::mat& X, const arma::vec& w, const arma::vec& y) {
  mat R = chol((X.each_col() % w).t() * X);
  vec p = solve(R, solve(R.t(), X.t() * (w % y), solve_opts::fast), solve_opts::fast);
  return p;
}

// [[Rcpp::export]]
arma::vec arma_chol2(const arma::mat& X, const arma::vec& w, const arma::vec& y) {
  mat XW = (X.each_col() % w).t();
  mat R = chol(XW * X);
  vec p = solve(R, solve(R.t(), XW * y, solve_opts::fast), solve_opts::fast);
  return p;
}

// [[Rcpp::export]]
arma::vec arma_direct1(const arma::mat& X, const arma::vec& w, const arma::vec& y) {
  vec sw = sqrt(w);
  vec p = solve(X.each_col() % sw, y % sw);
  return p;
}

// [[Rcpp::export]]
arma::vec arma_direct2(const arma::mat& X, const arma::vec& w, const arma::vec& y) {
  vec p = solve((X.each_col() % w).t() * X, X.t() * (w % y));
  return p;
}