paralelizar una regresión ventana deslizante en I

Question

Estoy corriendo una regresión rodadura muy similar al siguiente código:

library(PerformanceAnalytics) 
library(quantmod) 
data(managers) 

FL <- as.formula(Next(HAM1)~HAM1+HAM2+HAM3+HAM4) 
MyRegression <- function(df,FL) { 
    df <- as.data.frame(df) 
    model <- lm(FL,data=df[1:30,]) 
    predict(model,newdata=df[31,]) 
} 

system.time(Result <- rollapply(managers, 31, FUN="MyRegression",FL, 
    by.column = FALSE, align = "right", na.pad = TRUE)) 

Tengo algunos procesadores adicionales, así que estoy tratando de encontrar una manera de poner en paralelo la ventana rodante. Si esto fuera una regresión no rodante, podría paralelizarlo fácilmente usando la familia de funciones apply ...

Answer 1

La más obvia es usar lm.fit() en lugar de lm() para no incurrir en todos los gastos generales al procesar la fórmula, etc. .

actualización: Así que cuando dije obvio lo que quería decir era salta a la vista, pero engañosamente difícil de implementar!

Después de un poco de tocar el violín alrededor, me ocurrió con esta

library(PerformanceAnalytics) 
library(quantmod) 
data(managers) 

La primera etapa es darse cuenta de que la matriz de modelo puede ser creado previamente, por lo que hacer eso y convertirlo de nuevo en un objeto Zoológico de uso con rollapply():

mmat2 <- model.frame(Next(HAM1) ~ HAM1 + HAM2 + HAM3 + HAM4, data = managers, 
        na.action = na.pass) 
mmat2 <- cbind.data.frame(mmat2[,1], Intercept = 1, mmat2[,-1]) 
mmatZ <- as.zoo(mmat2) 

Ahora necesitamos una función que dará empleo a lm.fit() a hacer el trabajo pesado sin tener que crear matrices de diseño en cada iteración:

MyRegression2 <- function(Z) { 
    ## store value we want to predict for 
    pred <- Z[31, -1, drop = FALSE] 
    ## get rid of any rows with NA in training data 
    Z <- Z[1:30, ][!rowSums(is.na(Z[1:30,])) > 0, ] 
    ## Next() would lag and leave NA in row 30 for response 
    ## but we precomputed model matrix, so drop last row still in Z 
    Z <- Z[-nrow(Z),] 
    ## fit the model 
    fit <- lm.fit(Z[, -1, drop = FALSE], Z[,1]) 
    ## get things we need to predict, in case pivoting turned on in lm.fit 
    p <- fit$rank 
    p1 <- seq_len(p) 
    piv <- fit$qr$pivot[p1] 
    ## model coefficients 
    beta <- fit$coefficients 
    ## this gives the predicted value for row 31 of data passed in 
    drop(pred[, piv, drop = FALSE] %*% beta[piv]) 
} 

Una comparación de los tiempos:

> system.time(Result <- rollapply(managers, 31, FUN="MyRegression",FL, 
+         by.column = FALSE, align = "right", 
+         na.pad = TRUE)) 
    user system elapsed 
    0.925 0.002 1.020 
> 
> system.time(Result2 <- rollapply(mmatZ, 31, FUN = MyRegression2, 
+         by.column = FALSE, align = "right", 
+         na.pad = TRUE)) 
    user system elapsed 
    0.048 0.000 0.05 

donde se aprecia una mejora bastante razonable con respecto al original. Y ahora compruebe que los objetos resultantes son los mismos:

> all.equal(Result, Result2) 
[1] TRUE 

¡Disfrútelo!

Answer 2

Puede reducir el tiempo de ejecución al actualizar una descomposición. Esto produce un costo de en cada iteración en lugar de donde n es el ancho de la ventana. A continuación hay un código para comparar los dos. Es probable que sea mucho más rápido hacerlo en C++ pero los LINPACK dchud y dchdd no están incluidos con R, por lo que tendría que escribir un paquete para hacerlo. Además, Recuerdo haber leído que usted puede hacer más rápido con otras implementaciones que la LINPACK dchud y dchdd para la actualización R

library(SamplerCompare) # for LINPACK `chdd` and `chud` 
roll_forcast <- function(X, y, width){ 
    n <- nrow(X) 
    p <- ncol(X) 
    out <- rep(NA_real_, n) 

    is_first <- TRUE 
    i <- width 
    while(i < n){ 
    if(is_first){ 
     is_first <- FALSE 
     qr. <- qr(X[1:width, ]) 
     R <- qr.R(qr.) 

     # Use X^T for the rest 
     X <- t(X) 

     XtY <- drop(tcrossprod(y[1:width], X[, 1:width])) 
    } else { 
     x_new <- X[, i] 
     x_old <- X[, i - width] 

     # update R 
     R <- .Fortran(
     "dchud", R, p, p, x_new, 0., 0L, 0L, 
     0., 0., numeric(p), numeric(p), 
     PACKAGE = "SamplerCompare")[[1]] 

     # downdate R 
     R <- .Fortran(
     "dchdd", R, p, p, x_old, 0., 0L, 0L, 
     0., 0., numeric(p), numeric(p), integer(1), 
     PACKAGE = "SamplerCompare")[[1]] 

     # update XtY 
     XtY <- XtY + y[i] * x_new - y[i - width] * x_old 
    } 

    coef. <- .Internal(backsolve(R, XtY, p, TRUE, TRUE)) 
    coef. <- .Internal(backsolve(R, coef., p, TRUE, FALSE)) 

    i <- i + 1 
    out[i] <- X[, i] %*% coef. 
    } 

    out 
} 

# simulate data 
set.seed(101) 
n <- 10000 
wdth = 50 
X <- matrix(rnorm(10 * n), n, 10) 
y <- drop(X %*% runif(10)) + rnorm(n) 
Z <- cbind(y, X) 

# assign other function 
lm_version <- function(Z, width = wdth) { 
    pred <- Z[width + 1, -1, drop = FALSE] 
    ## fit the model 
    Z <- Z[-nrow(Z), ] 
    fit <- lm.fit(Z[, -1, drop = FALSE], Z[,1]) 
    ## get things we need to predict, in case pivoting turned on in lm.fit 
    p <- fit$rank 
    p1 <- seq_len(p) 
    piv <- fit$qr$pivot[p1] 
    ## model coefficients 
    beta <- fit$coefficients 
    ## this gives the predicted value for row 31 of data passed in 
    drop(pred[, piv, drop = FALSE] %*% beta[piv]) 
} 

# show that they yield the same 
library(zoo) 
all.equal(
    rollapply(Z, wdth + 1, FUN = lm_version, 
      by.column = FALSE, align = "right", fill = NA_real_), 
    roll_forcast(X, y, wdth)) 
#R> [1] TRUE 

# benchmark 
library(compiler) 
roll_forcast <- cmpfun(roll_forcast) 
lm_version <- cmpfun(lm_version) 
microbenchmark::microbenchmark(
    new = roll_forcast(X, y, wdth), 
    prev = rollapply(Z, wdth + 1, FUN = lm_version, 
        by.column = FALSE, align = "right", fill = NA_real_), 
    times = 10) 
#R> Unit: milliseconds 
#R> expr  min  lq  mean median  uq  max neval cld 
#R> new 113.7637 115.4498 129.6562 118.6540 122.4930 230.3414 10 a 
#R> prev 639.6499 674.1677 682.1996 678.6195 686.8816 763.8034 10 b