提问人：spellard 提问时间：7/27/2023 最后编辑：ThomasIsCodingspellard 更新时间：8/10/2023 访问量：96
在 r 中快速扩展系列。
Fast series expansion in r
问：
我想生成一个平滑函数的级数展开，类似于泰勒级数（但使用所谓的切比雪夫多项式）。但是，我还没有找到一种快速计算给定系数的级数展开的方法。特别是，对（切比雪夫/泰勒）多项式进行硬编码似乎比使用现有函数快得多。
显然，当近似阶数很高时，对多项式进行硬编码是不可行的，但问题是，如何用更快的代码实现这一点。这是我尝试过的：
library(mpoly)
library(foreach)
n=5
# Hard code the Chebyshev polynomials of order 1 to 4.
cheb = function(x) {x} # first order Cheb polynomial
cheb2 = function(x) {-1 + 2*x^2} # second order Chebyshev polynomial
cheb3 = function(x) {-3*x+4*x^3} # etc...
cheb4 = function(x) {1-8*x^2+8*x^4}
# Generate Chebyshev polynomials up to order n-1 using chebyshev() from library mpoly
cheb_funs = foreach(n_ind = 1:n) %do% as.function(chebyshev(n_ind-1))
# Calculate the series expansion using hard coded polynomials
Cheb_app1 = function(x)
{
p = x$p
v = x$v_cheb
v[1] + v[2]*cheb(p) + v[3]*cheb2(p) + v[4]*cheb3(p) + v[5]*cheb4(p)
}
# Calculate the series expansion using foreach() and cheb_funs
Cheb_app2 = function(x)
{
p = x$p;
v = x$v_cheb;
foreach(p_ind = p, .combine = "c") %do%
(foreach(n_ind = 1:n, .combine = "sum") %do%
(v[n_ind]*(cheb_funs[[n_ind]](p_ind))))
}
# Microbenchmark
z = list(p = runif(10, -1,1), v_cheb = runif(n, -1,1))
microbenchmark(Cheb_app1(z), Cheb_app2(z))
Unit: microseconds
expr min lq mean median uq max neval
Cheb_app1(z) 3.473 4.011 120.6585 7.955 12.944 11217.96 100
Cheb_app2(z) 30249.810 33243.943 35605.0237 34127.602 38195.954 47322.77 100
答：
0赞
spellard
7/27/2023
#1
最后，这是我想出的一个解决方案。
library(pracma)
library(microbenchmark)
rm(list = ls())
n=5
# Hard code the Chebyshev polynomials of order 1 to 4.
cheb = function(x) {x} # first order Cheb polynomial
cheb2 = function(x) {-1 + 2*x^2} # second order Chebyshev polynomial
cheb3 = function(x) {-3*x+4*x^3} # etc...
cheb4 = function(x) {1-8*x^2+8*x^4}
# Generate Chebyshev polynomials up to order n-1 using chebyshev() from library mpoly
cheb_funs = foreach(n_ind = 1:n) %do% as.function(chebyshev(n_ind-1))
# Calculate the series expansion using hard coded polynomials
Cheb_app1 = function(x)
{
p = x$p
v = x$v_cheb
v[1] + v[2]*cheb(p) + v[3]*cheb2(p) + v[4]*cheb3(p) + v[5]*cheb4(p)
}
## Another solution ##
power <- function(x, y) sign(x)^y * abs(x)^y
polyMat = chebPoly(n-1)
Cheb_app3 = function(x) {colSums(polyMat %*% t(outer(x$p, (n-1):0, power))*x$v_cheb)}
# Microbenchmark
z = list(p = runif(10, -1,1), v_cheb = runif(n, -1,1))
microbenchmark(Cheb_app1(z),Cheb_app3(z))
Unit: microseconds
expr min lq mean median uq max neval
Cheb_app1(z) 3.366 3.6930 115.26432 4.0750 4.3250 11087.282 100
Cheb_app3(z) 14.305 15.1085 66.43251 15.4965 15.8955 4943.902 100
0赞
jblood94
7/27/2023
#2
do.call增加了一些额外的开销，但时间非常相似：Cheb_app2
library(mpoly)
n <- 5
# Hard code the Chebyshev polynomials of order 1 to 4.
cheb <- function(x) {x} # first order Cheb polynomial
cheb2 <- function(x) {-1 + 2*x^2} # second order Chebyshev polynomial
cheb3 <- function(x) {-3*x+4*x^3} # etc...
cheb4 <- function(x) {1-8*x^2+8*x^4}
cheb_all <- as.function(
Reduce("+", lapply(0:(n - 1), \(i) chebyshev(i)*mp(paste0("v", i)))),
varorder = c("x", paste0("v", 0:(n - 1))),
vector = FALSE
)
#> f(x, v0, v1, v2, v3, v4)
cheb_allc <- compiler::cmpfun(cheb_all)
cheb_funs <- as.function(chebyshev(0:(n - 1)))
Cheb_app1 <- function(x) {
p <- x$p
v <- x$v_cheb
v[1] + v[2]*cheb(p) + v[3]*cheb2(p) + v[4]*cheb3(p) + v[5]*cheb4(p)
}
Cheb_app2 <- function(x) do.call("cheb_all", c(list(x$p), as.list(x$v_cheb)))
Cheb_app2c <- function(x) do.call("cheb_allc", c(list(x$p), as.list(x$v_cheb)))
Cheb_app3 <- function(x) c(cheb_funs(x$p) %*% x$v_cheb)
# Microbenchmark
z <- list(p = runif(10, -1,1), v_cheb = runif(n, -1,1))
microbenchmark::microbenchmark(
Cheb_app1(z),
Cheb_app2(z),
Cheb_app2c(z),
Cheb_app3(z),
check = "equal",
times = 1e3
)
#> Unit: microseconds
#> expr min lq mean median uq max neval
#> Cheb_app1(z) 3.2 3.5 14.1300 3.6 3.8 8363.5 1000
#> Cheb_app2(z) 8.6 8.9 11.0045 9.1 9.3 1515.1 1000
#> Cheb_app2c(z) 6.4 6.8 8.5209 7.0 7.3 1292.8 1000
#> Cheb_app3(z) 41.8 43.2 48.3163 44.0 45.3 2089.1 1000
评论
0赞
spellard
7/28/2023
谢谢！这似乎对我不起作用，但我在函数中出现以下错误：“错误：”Reduce（“+”， lapply（0：（n - 1）， \“” 中的意外输入cheb_all
0赞
jblood94
7/28/2023
也许你有一个旧版本的 R，请尝试将反斜杠更改为 .function
0赞
ThomasIsCoding
7/29/2023
#3
如果你不介意速度：你可以为切布舍夫多项式定义一个函数，如下所示（例如，两个选项）：
- 递归表示
chebpoly <- Vectorize(function(x, N) {
if (N %in% c(0, 1)) {
return(x^N)
}
2 * x * Recall(x, N - 1) - Recall(x, N - 2)
}, vectorize.args = "N")
- 三角表示
chebpoly <- Vectorize(function(x, N) {
ifelse(abs(x) <= 1,
cos(N * acos(x)),
ifelse(x > 1,
cosh(N * acosh(x)),
(-1)^N * cosh(N * acosh((-x)))
)
)
}, vectorize.args = "N")
然后，你可以通过Cheb_app0
Cheb_app <- function(x) {
p <- x$p
v <- x$v_cheb
colSums(v * t(chebpoly(p, seq_along(v) - 1)))
}
标杆
set.seed(0)
n <- 5
z <- list(p = runif(10, -1, 1), v_cheb = runif(n, -1, 1))
microbenchmark(
Cheb_app0(z),
Cheb_app1(z),
check = "equivalent"
)
给
Unit: microseconds
expr min lq mean median uq max neval
Cheb_app0(z) 67.0 69.85 148.407 73.45 82.85 6503.1 100
Cheb_app1(z) 3.4 3.70 174.139 4.20 5.00 16914.2 100
评论
library(foreach)Unit: nanoseconds expr min lq mean median uq max neval Cheb_app1(z) 3900 4000 4536.5 4100 4300 23800 1000 Cheb_app2(z) 200 200 341.0 200 300 15600 1000