如何使用 Rmpfr::mpfr 减少计算时间？

问题描述

z1f 和 z1f.mpfr 函数除了使用 Rmpfr 外，其他功能相同。然而，即使我使用比没有 mpfr 的默认值更小的 precBits，z1f.mpfr 的运行时间也明显高于 z1f。以下代码：

library(Rmpfr)

DecimalDigits=20

z1f.mpfr <- function(sample_freq,precBits = DecimalDigits)
{
  khat<-length(sample_freq)
  n<-sum(sample_freq)
  prod<-mpfr(rep(1,khat),precBits)
  zf<-mpfr(rep(0,n),precBits)
  for (v in 1:n){
    for (k in 1:khat){
      if (sample_freq[k]>=1){
        zf[v] = zf[v]+mpfr(sample_freq[k],precBits)/mpfr(n,precBits)*prod[k];
        prod[k] = prod[k]*(mpfr(1,precBits)-mpfr(sample_freq[k]-1,precBits)/mpfr(n-v,precBits));
      }
    }
  }
  return(zf)
}


z1f <- function(sample_freq)
{
  khat<-length(sample_freq)
  n<-sum(sample_freq)
  prod<-rep(1,khat)
  zf<-rep(0,n)
  for (v in 1:n){
    for (k in 1:khat){
      if (sample_freq[k]>=1){
        zf[v] = zf[v]+sample_freq[k]/n*prod[k];
        prod[k] = prod[k]*(1-(sample_freq[k]-1)/(n-v));
      }
    }
  }
  return(zf)
}

sample <- c(44,30,19,10,13,7,9,6,8,11,4,3,1,2,1)

> system.time(z1f.mpfr(sample))
                      user                     system                    elapsed 
16.21299999999999386091076  0.02100000000000079580786 16.24599999999918509274721 

> system.time(z1f(sample))
                       user                      system                     elapsed 
0.0010000000000047748471843 0.0000000000000000000000000 0.0009999999992942321114242 

我的问题有两个：

1. 为什么会有如此巨大的差异？
2. 我们可以抵消差异吗？

谢谢！

解决方法

您需要运行“profile”命令来查看它把时间花在哪里。您正在调用计算到某个数字精度的“z1f.mpfr”函数。正如您在结果中看到的那样，这是非常昂贵的。如果你真的认为你需要 20 位数字，而不是你在浮点数中得到的数字，那么这就是你付出的代价。您无法优化它，因为所有时间都花在您正在调用的函数上。您可能希望获取这些函数的源代码，以了解您正在执行多少代码，而不是仅使用浮点数。功能不是免费提供的，如果您使用了配置文件，您会看到第二个功能运行得如此之快，以至于几乎没有记录太多 CPU 使用率。

,

（我是包维护者：@Data Munger 是正确的：如果您不使用非常快的 CPU 内置的浮点算法而是基于 MPFR 的算法，那么您将付出相当大的代价。还有更多——我的 Rmpfr 包中的额外低效率原则上可以消除。为此，我需要精通 R 和 C 编程的合作者......）。

第二：你真的应该在使用它们之前阅读重要函数的帮助页面......在任何 R 包中。
在这种情况下，您已经错误地开始了：precBits 是精度位，位是二进制数字，不是十进制数字......所以您的 {{1} } 在你的函数中是..抱歉地说..无稽之谈。

要从位到数字或返回，您必须用 / 乘以或除以 log10(2) ~= 0.3，即

precBits = DecimalDigits

例如

bits <- ceiling(decimals / log10(2))

...是的，请使用 > (digits <- 10*(1:4)) [1] 10 20 30 40 > (bits <- ceiling(digits / log10(2))) [1] 34 67 100 133 > ：它对您所有的代码阅读器都更加友好，包括您未来的自己 ;-)