3
$\begingroup$

I apologize if this is the wrong place for this question; there are a number of potential points of failure each of which suggest either Math StackExchange or StackOverflow or here, but since the application is statistics, I think this might be the most appropriate subreddit.

I'm developing change point methods for linear models that require computing an estimate of the coefficients at each point in the (ordered) sample. The code also needs to be fast since I'm using it in extensive simulations. This suggests that I need to use Rcpp and RcppArmadillo to perform computations efficiently.

Estimating linear models via least squares should generally be done using QR decompositions, but I don't want to recompute from scratch a QR decomposition every time I add/remove a row from a matrix (which, considering that I'm computing a new linear model as I move through the ordered sample, a row will be added/removed from the pre-post split sample all the time). Turns out there are algorithms for updating the QR decomposition. I found nicely laid out algorithms in this paper, and I made a C++ implementation. For brevity I will provide code only for updating the QR decomposition when a row is added. I'm using C++ Armadillo and Rcpp:

Header:

#define QR_ADD_ROW(Q, R, d, u, mu) \
{ \
    (Q).insert_rows((Q).n_rows, 1, true); \
    (Q).insert_cols((Q).n_cols, 1, true); \
    (Q)((Q).n_rows - 1, (Q).n_cols - 1) = 1; \
    double cs[2]; \
    double muu = (mu); \
    for (unsigned int j = 0; j < (R).n_cols; ++j) { \
        givens((R)(j, j), (u)(j), cs); \
        (R)(j, j) = cs[0] * (R)(j, j) - cs[1] * (u)(j); \
        if (j < (R).n_cols - 1) { \
            (R).submat(j, j + 1, j, (R).n_cols - 1) = \
                  (R).submat(j, j + 1, j, (R).n_cols - 1) * cs[0] - \
                  (u).cols(j + 1, (R).n_cols - 1) * cs[1]; \
            (u).cols(j + 1, (u).n_cols - 1) = \
                  (R).submat(j, j + 1, j, (R).n_cols - 1) * cs[1] + \
                  (u).cols(j + 1, (u).n_cols - 1) * cs[0]; \
        } \
        arma::vec t1 = (Q).col(j); \
        arma::vec t2 = (Q).col((Q).n_cols - 1); \
        double tt1 = (d)(j); \
        double tt2 = muu; \
        (d)(j) = cs[0] * tt1 - cs[1] * tt2; \
        muu = cs[1] * tt1 + cs[0] * tt2; \
        (Q).col(j) = cs[0] * t1 - cs[1] * t2; \
        (Q).col((Q).n_cols - 1) = cs[1] * t1 + cs[0] * t2; \
    } \
    (R).insert_rows((R).n_rows, 1, true); \
    (d).insert_rows((d).n_rows, 1); \
    (d)((d).n_rows - 1) = muu; \
}

void givens(const double& a, const double& b, double arr[]);

Code

#include "header.h"
#include <RcppArmadillo.h>

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

/* Credit for these algorithms to:
 * Sven Hammarling and Craig Lucas. "Updating the $QR$ factorization and the 
 * least squares problem". MIMS Preprint:2008.111 (2008).
 */

void givens(const double& a, const double& b, double arr[]) {
    //                c,   s
    double res[] = {1.0, 0.0};
    double t;

    if (b != 0) {
        if (std::abs(b) >= std::abs(a)) {
            t = -a / b;
            res[1] = 1 / std::sqrt(1 + t * t);
            res[0] = res[1] * t;
        } else {
            t = -b / a;
            res[0] = 1 / std::sqrt(1 + t * t);
            res[1] = res[0] * t;
        }
    }

    arr[0] = res[0];
    arr[1] = res[1];
}

// [[Rcpp::export]]
Rcpp::List qr_add_row_seq(arma::mat X, arma::vec Y) {
    const unsigned int d = X.n_cols;
    const unsigned int n = X.n_rows;
    Rcpp::List res;
    if (n < d) {
        throw std::domain_error("Bad X passed; cannot have fewer rows than "
                                "columns");
    }
    unsigned int i;
    arma::rowvec u;
    double mu;

    arma::mat Q;
    arma::mat R;
    arma::vec r;

    qr(Q, R, X.rows(0, d - 1));
    r = Q.t() * Y.rows(0, d - 1);
    res.push_back(Rcpp::List::create(Rcpp::Named("Q") = Q,
                  Rcpp::Named("R") = R,
                  Rcpp::Named("d") = r));
    for (i = d; i < n; ++i) {
        u = X.row(i);
        mu = Y(i);
        QR_ADD_ROW(Q, R, r, u, mu);
        res.push_back(Rcpp::List::create(Rcpp::Named("Q") = Q,
                                         Rcpp::Named("R") = R,
                                         Rcpp::Named("d") = r));
    }

    return res;
}

Now we can run some R test code. In particular, let's see how solving the least squares problem using the last QR decomposition obtained with the above C++ code to the solution found with R.

library(Rcpp)
library(RcppArmadillo)

x <- rnorm(100)
y <- 1 + 2 * x + rnorm(100)
(fit <- lm(y ~ x))
X <- model.matrix(fit)
# Sum of square errors
sum((y - predict(fit))^2)

sourceCpp("code.cpp")
res <- qr_add_row_seq(X, y)
# Slightly different
(tmp <- with(res[[99]], solve(R[1:2,],d[1:2,])))
names(tmp) <- names(coef(fit))
fit2 <- fit
fit2$coefficients <- tmp
sum((y - predict(fit))^2)  # Larger than before

with(res[[99]], Q %*% R) - X  # Significant numerical errors

Running the above code reveals the following:

  1. My code with updating QR decompositions does not get the least-squares solution.
  2. We should have X == Q %*% R but my code is off in about the hundredths' place often.
  3. The two fits are close; it seems that my algorithm is essentially doing the right thing but is suffering numerical issues.

I also checked what happens with we do a complete QR decomposition using arma::qr() and it looks like Armadillo's algorithm gets a QR decomposition that is much better at recovering X than the end result of successive updates (error on the order of $10^{-16}$ as opposed to $10^{-2}$). So I'm doubting that the problem is different numerical accuracies in R and Armadillo.

So what is going wrong?

(I also appologize for writing a macro rather than a function; C++ is more foreign to me and I always wondered why C++ programmers like header-only libraries or why much functionality is in the macros, and I've since learned that actually macros are bad for more complex functionality at least in C++, while in C they may still have some place.)

$\endgroup$
1
  • $\begingroup$ Are you sure QR decomposition is the best way to go? Many people do least squares with stochastic gradient descent. With all the round-off error people sometimes get in the exact methods, sometimes it turns out that the iterative methods are actually more accurate, and they're definitely faster. $\endgroup$ May 1, 2020 at 19:48

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.