如何用apply或lapply替换多项式计算函数中的for循环？

首先看原函数代码：

est_poly <- function(x,y,slopes,n) {
  x_vals <- seq(x[1],x[2],length.out = n)
  p0 <- y[1]
  p1 <- slopes[1]
  p2 <- (3*(y[2]-y[1])/(x_vals[-1]-x_vals[1])-2*(slopes[1])-slopes[2])/(x_vals[-1]-x_vals[1])
  p3 <- (slopes[1]+slopes[2]-2*(y[2]-y[1])/(x_vals[-1]-x_vals[1]))/(x_vals[-1]-x_vals[1])^2
  x1 <- x_vals[1]
  result <- c()
  for (i in x_vals) {
    poly <- p0+p1*(i-x1)+p2*(i-x1)^2+p3*(i-x1)^3
    result <- append(result, poly[4])
  }
  return(matrix(data = c(x_vals,result), nrow = n))
}
x_eg2 <- c(0,1)
y_eg2 <- c(0,1)
slopes_eg <- c(0,2)
est_poly(x_eg2, y_eg2, slopes_eg, n=5)

先修正原函数的计算错误 #

原函数中p2和p3的计算逻辑有问题：使用x_vals[-1]-x_vals[1]会得到多个值（对应x_vals除第一个点外的所有点到起点的距离），导致p2和p3是向量，后续取poly[4]是依赖巧合得到结果。正确的计算应该基于区间总长度x[2]-x[1]。

用sapply替换for循环 #

修正后用sapply遍历x_vals计算每个点的多项式值：

est_poly_sapply <- function(x,y,slopes,n) {
  x_vals <- seq(x[1],x[2],length.out = n)
  # 计算区间总长度
  L <- x[2] - x[1]
  p0 <- y[1]
  p1 <- slopes[1]
  p2 <- (3*(y[2]-y[1])/L - 2*slopes[1] - slopes[2])/L
  p3 <- (slopes[1] + slopes[2] - 2*(y[2]-y[1])/L)/L^2
  x1 <- x_vals[1]
  
  # 用sapply替换for循环
  result <- sapply(x_vals, function(i) {
    p0 + p1*(i-x1) + p2*(i-x1)^2 + p3*(i-x1)^3
  })
  
  return(matrix(c(x_vals, result), nrow = n))
}

更高效的向量化实现（无需循环或apply） #

R天生支持向量运算，直接对整个x_vals向量计算，效率比apply更高：

est_poly_vectorized <- function(x,y,slopes,n) {
  x_vals <- seq(x[1],x[2],length.out = n)
  L <- x[2] - x[1]
  p0 <- y[1]
  p1 <- slopes[1]
  p2 <- (3*(y[2]-y[1])/L - 2*slopes[1] - slopes[2])/L
  p3 <- (slopes[1] + slopes[2] - 2*(y[2]-y[1])/L)/L^2
  x1 <- x_vals[1]
  
  # 直接向量运算，无需循环
  delta <- x_vals - x1
  result <- p0 + p1*delta + p2*delta^2 + p3*delta^3
  
  return(matrix(c(x_vals, result), nrow = n))
}

测试输出 #

运行以下代码会得到预期结果：

x_eg2 <- c(0,1)
y_eg2 <- c(0,1)
slopes_eg <- c(0,2)
est_poly_vectorized(x_eg2, y_eg2, slopes_eg, n=5)
# 输出：
     [,1] [,2]
[1,] 0.00 0.00
[2,] 0.25 0.06
[3,] 0.50 0.25
[4,] 0.75 0.56
[5,] 1.00 1.00

内容的提问来源于stack exchange，提问作者user19825372