0

我有以下难题。以下代码从向量池中获取一个向量，将向量绑定到矩阵并在新矩阵上执行函数并返回一个标量结果。

In2 <- diag(nXtr+1)
  mu <- array(1,c(dim(Xcal)[1],1))
  Y.hat.calib <- array(0,c(nC,1))  
  alpha <- array(0,c(nC,1))
  P = c()

  for (i in 1:dim(Xcal)[1]){
    Xtr2 <- rbind(Xtr,Xcal[i,])
    K2 <-(Xtr2%*%t(Xtr2)+1)^2
    rowCnt <- dim(Xtr2)[1]
    mu[i] <- sqrt(1 + t(c(rep(1,(rowCnt-1)),0))%*%solve(K2+a*In2)%*%K2%*%c(rep(0,(rowCnt-1)),1))

    #---------------------------------------------------------------------
    Y.hat.calib[i] <- kCal[,i]%*%solve(K + a*In)%*%Ytr
    alpha[i] <- (abs(Y.hat.calib[i] - Ycal[i]))/mu[i]
    P <- c(P,alpha[i])  
    #---------------------------------------------------------------------

我已经在需要的地方预先分配了，但真的需要摆脱循环，因为它太耗时了。我玩过各种想法，但无法想出办法来做到这一点。

任何帮助将一如既往地感激不尽。如果有什么我错过了，请告诉我。

4

1 回答 1

0

摆脱for循环不会自动使事情变得更快。您可以对代码进行的最大更改是solve尽可能避免计算；solve计算量非常大。

我没有尝试确保这没有错误，因为您没有提供示例数据。但是您可以遵循一般的想法：不要solve在每个循环中都执行，在逻辑上分离您对mufrom的计算alpha，并在可能的情况下用矩阵乘法替换逐列函数。

# Your mu boils down to this:
get.mu <- function(Xcal.i, Xtr2=Xtr2, a.times.In2=a.times.In2) {  
  Xtr2 <- rbind(Xtr, Xcal.i) 
  K2 <-(Xtr2 %*% t(Xtr2) + 1)^2
  first.solve  <- solve(K2 + a.times.In2) %*% K2
  bread <- c(rep(1, nrow(Xtr2-1)), 0)
  sqrt(1 + t(bread) %*% first.solve  %*% bread) # mu
}

In2 <- diag(nXtr+1)
# Constants that you recalculated every loop.
a.times.In2 <- a*In2
second.solve <- solve(K  + a*In ) %*% Ytr
# Y.hat.calib can be fully calculated in one matrix multiplication.
Y.hat.calib <- kCal %*% second.solve
# Which means that the difference is a constant:
Y.diff <- abs(Y.hat.calib - Ycal)

# So alpha and mu could be calculated like:
mu <- apply(X.cal, 2, get.mu)
alpha <- t(t(Y.diff) / mu.i)
于 2014-07-26T17:11:57.497 回答