Ускорение итеративной функции с 2 входами/выходами

Question

Использование R, мне интересно, какой лучший способ итеративно оценить функцию нескольких входов и выходов. Я мотивирован участков увиденных на: http://paulbourke.net/fractals/clifford/

Основные уравнения:

x_{n+1} = sin(A* y_n) + C* cos(A* x_n) 
y_{n+1} = sin(B* x_n) + D* cos(B* y_n) 

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

#Parameters 
A <- -1.4 
B <- 1.6 
C <- 1.0 
D <- 0.7 

n_iter <- 10000000 

#Initial values 
x0 <- 0 
y0 <- 0 

#function to calculate n+1 points 
cliff <- function(x,y){ 
    c(sin(A*y) + C*cos(A*x), sin(B*x) + D*cos(B*y)) 
} 

#matrix to store results 
res_mat <- matrix(0,nrow=n_iter,ncol=2) 

#recursive loop (definitely not the fastest way to do this?) 
for (i in 2:n_iter){ 
    res_mat[i,] <- cliff(res_mat[i-1,1],res_mat[i-1,2]) 
} 

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

Answer 1

Один вариант использования Rcpp; для таких итеративных функций, как это, где каждое новое значение является сложной функцией значения предыдущей итерации, это часто дает неплохие ускорения.

library(Rcpp) 
cliff.rcpp = cppFunction(" 
NumericMatrix cliff(int nIter, double A, double B, double C, double D) { 
    NumericMatrix x(nIter, 2); 
    for (int i=1; i < nIter; ++i) { 
    x(i,0) = sin(A*x(i-1,1)) + C*cos(A*x(i-1,0)); 
    x(i,1) = sin(B*x(i-1,0)) + D*cos(B*x(i-1,1)); 
    } 
    return x; 
}") 
cliff.rcpp(10, 1, 2, 3, 4) 
#    [,1]  [,2] 
# [1,] 0.0000000 0.0000000 
# [2,] 3.0000000 4.0000000 
# [3,] -3.7267800 -0.8614156 
# [4,] -3.2595913 -1.5266964 
# [5,] -3.9781665 -4.2182644 
# [6,] -1.1296464 -3.1953775 
# [7,] 1.3346977 3.2046776 
# [8,] 0.6386906 4.4230487 
# [9,] 1.4501988 -2.3914781 
# [10,] -0.3208062 0.5208984 

Мы можем видеть, что это возвращает одинаковые результаты к коду в вопросе:

cliff.orig <- function(n_iter, A, B, C, D) { 
    #function to calculate n+1 points 
    cliff <- function(x,y){ 
    c(sin(A*y) + C*cos(A*x), sin(B*x) + D*cos(B*y)) 
    } 

    #matrix to store results 
    res_mat <- matrix(0,nrow=n_iter,ncol=2) 

    #recursive loop (definitely not the fastest way to do this?) 
    for (i in 2:n_iter){ 
    res_mat[i,] <- cliff(res_mat[i-1,1],res_mat[i-1,2]) 
    } 
    res_mat 
} 
identical(cliff.rcpp(10, 1, 2, 3, 4), cliff.orig(10, 1, 2, 3, 4)) 
# [1] TRUE 

Для входа в оригинальный вопрос, подход Rcpp дает ~ 50 раз SpeedUp:

system.time(cliff.rcpp(10000000, -1.4, 1.6, 1.0, 0.7)) 
# user system elapsed 
# 0.661 0.046 0.717 
system.time(cliff.orig(10000000, -1.4, 1.6, 1.0, 0.7)) 
# user system elapsed 
# 34.591 0.245 35.040