wrathematics

library(memuse)
library(hdfio)

f = nycflights13::flights

for (compression in c(0, 4, 9)){
  suppressWarnings(file.remove("flights.h5"))
  write_h5df(f, "flights.h5", compression=compression)
  print(Sys.filesize("flights.h5"))
}

wrathematics

library(inline)
csrc = "
  SEXP ret;
  PROTECT(ret = allocVector(LGLSXP, 2));

  LOGICAL(ret)[0] =  1;
  LOGICAL(ret)[1] = -1;

  UNPROTECT(1);
  return ret;

wrathematics

---
title: "A Bit on the Numerics Behind R's Linear Model Fitters"
output: html_document
---

The classical linear regression setup is that you want to solve a system of equations that looks something like:
$X_{m\times n}\beta_{n\times 1} = y_{m\times 1}$

In statistics/data analysis, it's pretty typical for $m>n$. If $X$ has ["full rank"](https://en.wikipedia.org/wiki/Rank_(linear_algebra)) (more on this later). In this case, [we can derive a closed form solution](https://en.wikipedia.org/wiki/Least_squares#Linear_least_squares) for $\beta$:

wrathematics

#include <stdio.h>
#include <papi.h>
#include <stdlib.h>

#define NUM_EVENTS 2

int main(int argc, char **argv)
{
  int a = atoi(argv[1]);
  int b = atoi(argv[2]);

wrathematics

#include <stdio.h>
#include <papi.h>

#define NUM_EVENTS 2

int main()
{
  int events[NUM_EVENTS] = {PAPI_L1_DCM, PAPI_L2_DCM};
  long long values[NUM_EVENTS];
  

wrathematics

// ----- count.c
// build with R CMD SHLIB count.c (or put it in the package, w/e)

#include <R.h>
#include <Rinternals.h>

SEXP colsum_which_eq0(SEXP nrows_, SEXP ncols_, SEXP J_)
{
  SEXP ret;
  const int nrows = INTEGER(nrows_)[0];

wrathematics

00000000001319f0 <Rf_unprotect@@Base>:
  1319f0:       48 8b 05 21 b1 39 00    mov    0x39b121(%rip),%rax        # 4ccb18 <R_PPStackTop@@Base-0x15b920>
  1319f7:       29 38                   sub    %edi,(%rax)
  1319f9:       c3                      retq   
  1319fa:       66 0f 1f 44 00 00       nopw   0x0(%rax,%rax,1)

wrathematics

static inline int get_rank(int m, int n, double *qr, double tol)
{
  int i;
  const double minval = fabs(tol*qr[0]);
  const int minmn = m > n ? n : m;
  
  for (i=1; i<minmn; i++)
  {
    if (fabs(qr[i + m*i]) < minval)
      return i;

wrathematics

### Keybase proof

I hereby claim:

  * I am wrathematics on github.
  * I am wrathematics (https://keybase.io/wrathematics) on keybase.
  * I have a public key ASBGgOTJJPQ52zRvzAj47Z_ZAN4LzTlqGyx8BiFn6fAQLQo

To claim this, I am signing this object:

wrathematics

PKG_LIBS = -L. -lgotest

OBJECTS = gotest/test.o

all: $(SHLIB)

$(SHLIB): $(OBJECTS)


GOLIB = libgotest.so