2

I need to compute a hat matrix (as from linear regression). Standard R code would be:

H <- tcrossprod(tcrossprod(X, solve(crossprod(X))), X)

with X being a relatively large matrix (i.e 1e5*100), and this line has to run thousands of times. I understand the most limiting part is the inverse computation, but the crossproducts may be time-consuming too. Is there any faster alternative to perform these matrix operations? I tried Rcpp and reviewed several posts but any alternative I tested was slower. Maybe I did not code properly my C++ code, as I am not an advanced C++ programmer.

Thanks!

DGMartin
  • 67
  • 7
  • 1
    `crossprod(X)` is a 100x100 matrix, the `solve()` part is therefore instantaneous. What takes times is to compute `tcrossprod(same_size_as_X, X)` at the end. Are you sure this is really what you need to do? Are you not doing anything with this matrix afterwards? E.g. multiplying it with a vector or getting the diagonal? – F. Privé Oct 23 '20 at 14:31
  • 2
    general comment: it's often a mistake to solve statistical problems by directly translating the linear algebra. For example, the 'hatvalues' (diagonal of the hat matrix) are computed in R starting from the results of the QR decomposition: see [here](https://github.com/wch/r-source/blob/3afaa607c48ab59ac8230ddda298cff326e58671/src/library/stats/src/lminfl.f) (if you need the whole hat matrix you can't use exactly that solution, but the general point that there's probably more efficient linear algebra that will get what you want holds ...) – Ben Bolker Oct 23 '20 at 14:32
  • For the general question of "how to do matrix operations fast?" it's single thread MKL + data parallel (e.g., if you're doing a simulation, run multiple iterations in parallel). – thc Oct 23 '20 at 16:45
  • I see, thank you all for the comments, also @eddelbuettel for the response below. As this will be part of an R package probably changing the 'reference BLAST' is not the way to go, as each user may have different libraries.. I'll give a try to find a solution that does not directly translate the linear algebra, as pointed out by Ben. Indeed this line is part of the calculation of the sums of squares crossproduct matrix in multivariate linear models. I asked again [here](https://stackoverflow.com/questions/64534242/faster-alternative-to-r-caranova-for-sum-of-square-crossproduct-matrix-calcula). – DGMartin Oct 26 '20 at 12:02

1 Answers1

13

Chasing the code for this line by line is a little difficult because the setup of R code is a little on the complicated side. But read on, pointers below.

The important part is that the topic has been discussed many times: what happens is that R dispatches this to the BLAS (Basic Linear Algebra Subprogram) and LAPACK (Linear Algebra PACKage) libraries. Which contain the most efficient code known to man for this. In general, you cannot gain on it by rewriting.

One can gain performance differences by switching one BLAS/LAPACK implementation for another---there are many, many posts on this online too. R itself comes with the so-called 'reference BLAS' known to be correct, but slowest. You can switch to Atlas, OpenBLAS, MKL, ... depending on your operating system; instructions on how to do so are in some of the R manuals that come with your installation.

For completeness, per file src/main/names.c the commands %*%, crossprod and tcrossprod all refer to do_matprod. This is in file src/main/array.c and does much argument checking and arranging and branching on types of arguments but e.g. one path then calls

F77_CALL(dsyrk)(uplo, trans, &nc, &nr, &one, x, &nr, &zero, z, &nc
        FCONE FCONE);

which is this LAPACK function. It is essentially the same for all others making this an unlikely venue for your optimisation.

Dirk Eddelbuettel
  • 360,940
  • 56
  • 644
  • 725