Efficient routines for a regression with orthogonal regressors? I have a standard OLS regression setup, where (sets of) the regresors are orthogonal to each other. I am looking for a fast low-level way (using qr() instead of lm() ) to do this in R. So far, I can think of at least three approaches:


*

*Simple: Run the standard OLS: should be slow as does not make use of the orthogonal structure

*Separate: Run regressions separately: we know that with orthogonal regressors, we could get the same estimates by estimating OLS on each (sets of ) regressors separately. 

*Sparse: Use a sparse routine:  the X matrix contains many zeros, so could be considered sparse. Or use a "block diagonal" routine in R for inverting the $X^\prime X$ noting that with (sets of ) orthogonal regressors, it is (block) diagonal.
Against my expectations, in a single experiment (code below) I ran, the simplest method (1, simple) was the fastest, while the fancy sparse method the slowest.  Any ideas of why this would happen? Any idea of a better way/different routines that make use of the block-diagonal structure?
Thanks!
library(Matrix)
library(microbenchmark)


n <- 1000
k <- 5
x <- matrix(rnorm(n*k), ncol=k)
eps <- rnorm(n)
y <- as.matrix(ifelse(1:n<(n/2), 
                  x %*% seq(0.8, by=0.05,length.out=k)+eps, 
                  x %*% seq(0.5, by=0.05,length.out=k)+eps))


X_s <- bdiag(x[1:(n/2),], x[(n/2+1):n,])
X <- as.matrix(X_s)
microbenchmark(Simple=qr.coef(qr(X), y),
           Separate=c(qr.coef(qr(X[,1]), y), qr.coef(qr(X[,2]), y)),
           Sparse=qr.coef(qr(X_s), y))

 A: The reason the sparse method is so much slower is that you are using a rather dense matrix as a sparse matrix. 50% non-zero elements are just too many for a matrix to be considered sparse. You really should aim for ~5% so you can enjoy obvious and consistent speed-ups. In general, I would not expect a sparse method to be competitive to a dense method unless the zero entries were at least around the 85%-90% of the matrix.
Having said that what you do is correct to use QR. I would suggest using qr.solve directly so I save one extra assignment but other than that you are using the standard way to solve the normal equations.
If you really want to be aggressive, try to use Cholesky $LL^T$ decomposition. Using the triangular matrix $L$ you can back/forwardsolve the system at hand. This is a very fast operation. So something like in R:
cholCoef <- function(X,y){
              U <- chol(crossprod(X)); 
              backsolve( U, backsolve( transpose=TRUE, U, x= crossprod(X,y))) 
            }

Please note that this can be numerically unstable and fail horribly if the matrix crossprod(X) is not positive definite or the condition number associated with $X$ is already high but hey... we want speed no stability. Generally speaking the Cholesky's complexity is of the order $\frac{1}{3} n^3$ FLOPS while QR's complexity is of the order  $\frac{4}{3} n^3$. (See Trefethen & Bau, 1997 Numerical linear algebra for more details).
Finally it worth noting that the package RcppEigen has fastLM function you may want to consider. The argument method defines the actual algorithm to use. Probably the fastLM( ...,  method = 2) is good enough to almost general use. It will do an $LDL^T$ decomposition instead of the standard Cholesky: $LL^T$ so it rather stable. In my system CholSolve is the fastest using the following microbenchmark.
microbenchmark( 
  Simple= qr.solve(X, y), 
  CholSolve = cholCoef( X = X, y),  
  Separate=c(qr.coef(qr(X[,1]), y), qr.coef(qr(X[,2]), y)), 
  Sparse=qr.coef(qr(X_s), y),  
  FastLM = fastLm(X,y, method= 2L))

As mentioned by @Cliff the matrices you deal with are rather small so you would probably want to use something bigger to evaluate the differences between different algorithms' numerical performance. Finally note that the BLAS variant use by your R installation will make a very big difference. I would strongly urge you to use OpenBLAS or ATLAS if you are not using them already. (See my comment to check which BLAS you use.)
