Package Examples

[coatless]

The usual fastLM() approach for Armadillo is albeit problematic as we do not have support for non-square matrices in solve() yet.

That is, we wanted RcppArmadillo::fastLm():

#include <RcppArmadillo.H>

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
Rcpp::List fastLm_impl(const arma::mat& X, const arma::colvec& y) {
    int n = X.n_rows, k = X.n_cols;

    arma::colvec coef = arma::solve(X, y);     // fit model y ~ X
    arma::colvec res  = y - X*coef;            // residuals
    double s2 = arma::dot(res, res) / (n - k); // std.errors of coefficients
    arma::colvec std_err = arma::sqrt(s2 * arma::diagvec(arma::pinv(arma::trans(X)*X)));

    return Rcpp::List::create(Rcpp::Named("coefficients") = coef,
                              Rcpp::Named("stderr")       = std_err,
                              Rcpp::Named("df.residual")  = n - k);
}

To go to:

#include <RcppBandicoot.h>

// [[Rcpp::depends(RcppBandicoot)]]

// [[Rcpp::export]]
Rcpp::List fastLm_coot(const coot::mat& X, const coot::colvec& y) {
    coot::uword n = X.n_rows, k = X.n_cols;

    coot::colvec coef = coot::solve(X, y);        // fit model y ~ X
    coot::colvec res = y - X * coef;              // residuals

    double s2 = coot::dot(res, res) / (n - k);    // std.errors of coefficients

    coot::colvec std_err = coot::sqrt(s2 * coot::diagvec(coot::pinv(coot::trans(X) * X)));

    // Stresses the custom as<>() and wrap() functions.
    return Rcpp::List::create(
        Rcpp::Named("coefficients") = coef,
        Rcpp::Named("stderr")       = std_err,
        Rcpp::Named("df.residual")  = n - k
    );
}

In the interim, some quick ideas:

1. Gradient descent for OLS or Logistic
   • Unlikely to beat project speed for OLS,
   • See demo for Logistic (requires Ensmallen) https://coot.sourceforge.io/speed.html
2. PCA via Covariance Eigendecomposition: X' * X gives a p x p covariance (square!)
   • Task 4-ish
3. Large matrix ops (axpy / multiplication)
   • Task 1-3

Still need the "killer demo".