解释 Rprofile 输出：<Anonymous> 函数是什么？

Question

所以我有一个运行 MCMC 算法的大函数。我相信大多数昂贵的操作是大矩阵的乘法，但是这个 Rprof 输出相当令人困惑。

$by.self
                self.time self.pct total.time total.pct
"<Anonymous>"      328.90    81.84     329.34     81.95
"fprod"             46.16    11.49     376.02     93.57
"Dikin_Walk"         7.42     1.85     401.32     99.86
"as.vector"          5.98     1.49      57.56     14.32
".External"          2.54     0.63       2.54      0.63
"-"                  1.84     0.46       1.84      0.46
"H_x"                1.16     0.29     225.82     56.19
"fcrossprod"         1.14     0.28     226.12     56.27

编辑：这是我在大包装函数中定义的 3 个函数：

## first, augment A | b 
A_b <- cbind (b, A)

## H(x) is the hessian
H_x <- function(x) {
  D <- as.vector(1/(A_b[,1] - fprod(A_b[,-1], x)))
  D_squared <- fdiag(D^2)
  return(fcrossprod(A, fprod(D_squared, A)))  
}

## D(x) is the diagonalized matrix of the log-barrier function of Ax <= b
D_x <- function(x) {

  D <- as.vector(1/(A_b[,1] - fprod(A_b[,-1], x)))
  return(fdiag(D))
} 

## checks whether a point z is in Ellip(x)
ellipsoid <- function(z, x) {

  ## as.numeric converts the expression into an atom, so we get boolean
  return( as.numeric(fcrossprod(z-x, fprod(H_x(x), (z-x)))) <= r^2)

}

fdiag、fcrossprod 和 fprod 都是 R 对应版本的 RcppArmEigen 版本。我使用它们是因为它们快得多。

主要算法：

> for (i in 1:n) {
>   
>    zeta <- rnorm(length(b), 0, 1)
>    zeta <- r * zeta / sqrt(as.numeric(fcrossprod(zeta,zeta)))
>  
>    rhs <- fcrossprod(A, fprod(D_x(current.point), zeta))
>    
>    ## DONE
>    
>    y <- fprod(fsolve(H_x(current.point)), rhs) 
>    y <- y + current.point
>    
>
>    while(!ellipsoid(current.point, y)) {
>    zeta <- rnorm(length(b), 0, 1)
>      
>     ## normalise to be on the m- unit sphere
>     ## and then compute lhs as a m-vector
>     zeta <- r * zeta / sqrt(sum(zeta * zeta))
>     
>     
>     rhs <- fcrossprod(A, fprod(D_x(current.point), zeta))
>     
>     ## 
>     y <- fprod(fsolve(H_x(current.point)), rhs)  
>     y <- y + current.point
>     
>
>     if(ellipsoid(current.point, y)) {
>      
>       probability <- min(1, sqrt(fdet(fprod(fsolve(H_x(current.point)),H_x(y))      )))
>        
>        
>       bool <- sample(c(TRUE, FALSE), 1, prob = c(probability, 1-?>probability))
>       if(bool) {
>         break
>       } 
>     }
>   }

这里是 by.total 输出：

$by.total
                      total.time total.pct self.time self.pct
"Dikin_Walk"              401.32     99.86      7.42     1.85
"fprod"                   376.02     93.57     46.16    11.49
"<Anonymous>"             329.34     81.95    328.90    81.84
"cbind"                   268.58     66.83      0.04     0.01
"fcrossprod"              226.12     56.27      1.14     0.28
"H_x"                     225.82     56.19      1.16     0.29
"fsolve"                  203.82     50.72      0.14     0.03
"ellipsoid"               126.30     31.43      0.56     0.14
"fdet"                     64.84     16.13      0.02     0.00
"as.vector"                57.56     14.32      5.98     1.49
"fdiag"                    35.68      8.88      0.50     0.12

fprod 定义为：

prodCpp <- 'typedef Eigen::Map<Eigen::MatrixXd> MapMatd;
const MapMatd B(as<MapMatd>(BB));
const MapMatd C(as<MapMatd>(CC));
return wrap(B * C);'

fprod <- cxxfunction(signature(BB = "matrix", CC = "matrix"),
                      prodCpp, "RcppEigen")

Answer 1

<Anonymous> 指的是匿名（未命名）函数。如果你在运行这样的函数中循环，一般大部分时间都会花在这个函数上。

显然 A_b 是一个矩阵，x 是一个向量。使用矩阵代数代替循环：

A_b <- matrix(1:16, 4)
x <- 1:3
D <- apply(A_b, 1, function(row) {1 / (row[1] - sum(row[-1] * x))})

D1 <- as.vector(1/(A_b[,1] - A_b[,-1] %*% x))

identical (D, D1)
#[1] TRUE

编辑：

匿名函数在fprod的Rcpp魔法中：

B <- matrix(rnorm(1e6),1e3)
C <- matrix(rnorm(1e6),1e3)

Rprof()
for (i in 1:30) BC <- fprod(B, C)
Rprof(NULL)
summaryRprof()
#$by.self
#              self.time self.pct total.time total.pct
#"<Anonymous>"      4.24      100       4.24       100
#
#$by.total
#              total.time total.pct self.time self.pct
#"<Anonymous>"       4.24       100      4.24      100
#"fprod"             4.24       100      0.00        0
#
#$sample.interval
#[1] 0.02
#
#$sampling.time
#[1] 4.24

你的大部分时间都花在了矩阵乘法上。您可能会受益于优化的 BLAS，例如，您可以尝试 OpenBLAS。

Answer 2

首先，忽略 "self time"，因为 "total time" 包含了那个加上被调用者。如果您花费任何不需要的时间，您更有可能通过调用函数来完成它，而不是通过处理。**

其次，别看那个。 Rprofile 生成堆栈跟踪文件。看看其中的几个，随机选择。如果一个函数负责 80% 的时间，您将在 5 个堆栈跟踪中大约有 4 个看到它。更重要的是，您会看到谁在调用它，您会看到它在调用谁，以花费时间。简单的数字不会告诉你这一点。对堆栈跟踪进行排序也不会告诉您这一点。

如果它给出调用的行号就更好了，但它没有。即便如此，只是展示功能还是很有用的。

** 分析器只显示 "self time"，因为他们总是这样做，而且其他人都这样做，很少有人意识到这只是一种干扰。如果函数位于堆栈跟踪的终点，则它位于 "self time" 中。无论哪种方式，它都在 "inclusive time".

解释 Rprofile 输出：<Anonymous> 函数是什么？

Interpreting Rprofile Output: What is this <Anonymous> Function?

profiling

r

anonymous-function