R: Нужны советы по оптимизации алгоритма обратного выбора

Question

Я пытаюсь оптимизировать небольшой алгоритм, который я сделал для выполнения обратного выбора в области регрессии сплайна. Основное, что делает алгоритм:

Возьмите k вектор узлов с n компонентами.

Устранить i-й компонент, i=1,...,n.

Оценить сплайн-регрессию, используя вектор узла k[-i], i=1,...,n.

Выберите вариант с меньшей Остаточной суммой квадратов (RSS) и оцените BIC для этой модели.

Пусть k=k[-which.min(RSS)].

Запустите алгоритм снова до n=1.

Мой код

library(splines)
prune<-function(K,y0,x,deg=3){
  KNOTS<-matrix(nrow = (length(K)),ncol=(length(K)-1))
  y<-y0
  BIC<-vector(length =(length(K)-1) )

  for(j in 1:(length(K)-1)){
    RSS<-vector(1:(length(K)))
    {for(i in 1:(length(K)))
      RSS[i]<-sum((y-lm(y~bs(x,knots = K[-i],Boundary.knots = c(min(x),max(x)),degree = deg,intercept = T)-1)$fitted.values)^2)
    }
    K<-K[-which.min(RSS)]
    for (k in 1:(length(K)))
      KNOTS[k,j]<-K[k]

    BIC[j]<-BIC(lm(y~bs(x,knots = K,Boundary.knots = c(min(x),max(x)),degree = deg,intercept = T)-1))
  }

  KSTAR<-as.vector(KNOTS[, which.min(BIC)])
  I<-which.min(BIC)
  KSTAR<-na.omit(KSTAR)  
  return(list(KSTAR,KNOTS,I))
}

, где K - вектор узла, y0 - зависимая переменная, x - независимая переменная и deg - степень B-сплайнов. Проблема в том, что когда у меня много наблюдений, скажем, 5000 или 10000, для вычисления результата требуется слишком много.

Мои навыки программирования на R довольно просты, и, вероятно, есть вещи, которые можно написать по-другому, чтобы увеличить скорость.

Не могли бы вы дать мне совет?

Я читал некоторые блоги об ускорении кода R, но, честно говоря, я не понимаю, как можно применить такие вещи, как предварительное распределение, векторизация и т. Д.

для примера можно начать с

x=(0:4000)/4000  
y=sin(4*x)+2*exp(-30*(4*x)^2)+rnorm(4000,0,0.1)
K=seq(0,4000,by=100)

Я не прошу вас делать мою работу, просто ищу совет относительно структуры алгоритма. Заранее спасибо!

Answer 1

Большая часть времени, затраченного на ваш код, приходится на повторные вызовы lm. Это можно увидеть, если вы попробуете:

N <- 4000
x=(0:N)/N
y=sin(4*x)+2*exp(-30*(4*x)^2)+rnorm(N+1,0,0.1)
K=seq(0,N,by=100)
library(profvis)
profvis(prune1(K,x,y))

Если вы прервете вызов функции,

lm(y~bs(x,knots = K[-i],Boundary.knots = c(min(x),max(x)),degree = deg,intercept = T)-1)

как:

x1 <- bs(x,knots = K,degree = 3,intercept = T)
  lmod <- lm(y~x1-1)

Вывод bs, x1 имеет значения только в первых нескольких столбцах. Остальные все ноль. Кроме того, поскольку вам на самом деле не нужны дополнительные сведения, предоставляемые lm, вы можете дополнительно сократить их до самого необходимого. Как только вы это сделаете, ваша функция будет намного быстрее:

library(rbenchmark)
library(MASS)
benchmark(a={
  x1 <- bs(x,knots = K,degree = 3,intercept = T)
  lmod <- lm(y~x1-1)
  RSS<-sum(lmod$residuals**2)
},
b={
  x1 <- bs(x,knots = K,degree = 3,intercept = T)
  x2 <- x1[,which(colSums(abs(x1))>0)]   # Removing zero columns
  x21 <- ginv(x2)            # Simplified lm. If you don't want to risk it,
  y1 <- x2 %*%(x21  %*% y)   # you can try: lmod <- lm(y~x2-1)
  RSS <- sum((y-y1)^2)       #              RSS<-sum(lmod$residuals**2)
},replications = 9)          
#   test replications elapsed relative user.self sys.self user.child sys.child
# 1    a            9   0.187    5.054     0.171    0.016          0         0
# 2    b            9   0.037    1.000     0.036    0.000          0         0
all.equal(lmod$residuals,as.numeric(y-y1),check.attributes = F)
# [1] TRUE

Это должно быть еще быстрее для больших объемов данных. Ваша последняя функция будет выглядеть примерно так:

prune1<-function(K,y0,x,deg=3){
  KNOTS<-matrix(nrow = (length(K)),ncol=(length(K)-1))
  y<-y0
  BIC<-numeric(length =(length(K)-1) )
  kmin <- 0
  for(j in 1:(length(K)-1)){
    min_rss=Inf
    for(i in 1:(length(K))){
      x1 <- bs(x,knots = K,degree = 3,intercept = T)
      x2 <- x1[,which(colSums(abs(x1))>0)]
      x21 <- ginv(x2)
      y1 <- x2 %*%(x21  %*% y)
      RSS <- sum((y-y1)^2)
      if(RSS<min_rss){
        min_rss <- RSS
        kmin <- i
      }
    }
    K<-K[-kmin]
    KNOTS[1:length(K),j]<-K
    x1 <- bs(x,knots = K,degree = deg,intercept = T)
    x2 <- x1[,which(colSums(abs(x1))>0)]
    lmod <- lm(y~x2-1)
    BIC[j]<-BIC(lmod)
  }

  KSTAR<-as.vector(KNOTS[, which.min(BIC)])
  I<-which.min(BIC)
  KSTAR<-na.omit(KSTAR)  
  return(list(KSTAR,KNOTS,I))
}

benchmark(a={l <- prune(K,y,x)},
          b={l1 <- prune1(K,y,x)},
          replications = 1)
# test replications elapsed relative user.self sys.self user.child sys.child
# 1    a            1  15.480    3.091    15.483        0          0         0
# 2    b            1   5.008    1.000     5.008        0          0         0

Я также опробовал функцию для 10000 наблюдений и 100 узлов и получил время выполнения 3 минуты