【翻译】R 中的设计模式

R 中的设计模式

R 中的设计模式

本文翻译自 Design Patterns in R（By Sebastian Warnholz）。

本文的灵感来源于：

Stuart Sierra 的演讲，关于函数式编程中的设计模式；以及
我从 F# for fun an profit 想到的想法；以及
我在使用 R 的过程中用不同方法解决问题获得的反馈。

设计模式似乎是一个很大的词，特别是因为它在面向对象编程中的使用。但最终我认为它只不过是软件设计中的可重复策略。

不动点算法

下面，我使用 R 并且以计算正数平方根的不动点算法为例。算法定义如下：

\[x_{n+1} = f(x_n) \]

用于寻找平方根的不动点函数如下：

\[f(x \mid p) = \frac{p}{x} \]

要计算的正是 p 的平方根。用 R 代码描述算法：

fp <- function(f,
               x,
               converged,
               ...)
{
    value <- f(x, ...)
    if (converged(x, value)) value
    else Recall(f, value, converged, ...) 
}

x 是初始值或最后一次迭代得到的值。converged 是有 arguments 和 ... 两个参数的函数，... 用于 R 的“函数柯里化”^[1]。不动点函数如下：

fpsqrt <- function(x, p) p / x

并且

converged <- function(x, y) all(abs(x - y) < 0.001)

开始计算：

fp(fpsqrt, 2, converged, p = 2)

## Error: evaluation nested too deeply: infinite recursion / options(expressions=)?

第一次运行完全不起作用。在目前的代码实现中很难找出哪里出了问题，但我们会找到的。下面，我将应用不同的模式来修改上述框架以获得解决方案。

包装器模式

这个模式是我从 Stuart Sierra 的演讲中得到的。通过包装器模式我可以向一个函数增加新功能，却不改变原先的函数。我所要做的事就是给函数添加日志，或者记录函数的保留属性，R 中的许多函数并不会这样做。有时候要尝试调用一个函数并每两分钟重试一次，因为连接数据库失败或文件系统没有响应。

一个函数要有单一、明确的用途，日志和写入数据库是两件事。计算下一次迭代以及记录迭代次数也是两件事。这一个功能，以及用于记录是/否的额外参数，不会长期独立存在。

在我的上个例子中，遇到的问题是不动点函数在两个值之间振荡，而不是收敛于平方根。解决这个问题的一个技巧是使用平均阻尼。这就是说我们用 \(\frac{x_{n-1}+x_n}{2}\)，而不是 \(x_n\) 来计算 \(x_{n+1}\)。这实际上不是不动点函数逻辑的一部分，所以逻辑不应该被它污染：

averageDamp <- function(fun)
{
    function(x, ...) (x + fun(x, ...)) / 2
}

fp(averageDamp(fpsqrt), 2, converged, p = 2)

## [1] 1.414214

# and to compare:
sqrt(2)	
	
## [1] 1.414214

OK，看起来能跑通了！我需要一个额外的包装器来打印每一次迭代的值：

printValue <- function(fun)
{
    function(x,
             ...)
    {
        cat(x, "\n")
        fun(x, ...)
    }
}

fp(printValue(averageDamp(fpsqrt)), 2, converged, p = 2)

## 2 
## 1.5 
## 1.416667 
## 1.414216

## [1] 1.414214

现在的问题是，如果我们添加太多的包装器，计算就会变得复杂。事实上试图找出哪个包装器首先被调用，也许这对你来说并不明显。

包装器模式可以在原函数之前或之后（或前后同时）增加新功能。printValue 在原函数之前添加打印功能，averageDamp 在原函数之后作修正。如果看到 averageDamp 的另一种实现，模式将会显得更加清晰：

averageDamp <- function(fun)
{
    function(x, ...)
    {
        value <- fun(x, ...)
        (x + value) / 2
    }
}

接口模式

柯里化（Currying）

从我的角度来看，这项技术的价值在于你可以更轻松地构建接口（在语言能够充分支持的情况下）。例如，平方根的不动点函数需要两个参数。然而，该算法实际上只知道一个参数。在这种情况下柯里化只是意味着将双参数函数 fpsqrt 变成单参数函数。我们可以通过设置 p = 2 来实现这一点，这是我借助 ... 实现的。

在 R 中，有两个原生选项来模拟柯里化。通常看到的是使用点参数（...）来允许将额外的参数传递给该函数。然而，这给我的框架中的每个实现都带来了额外的负担，因为我需要让我定义的每个包装器函数使用点参数。另一种选择是使用匿名函数将原始版本封装在单参数函数中，如下所示：

fp(averageDamp(function(x) fpsqrt(x, p = 2)), 2, converged)
	
## [1] 1.414214

如果我依靠这个接口（单参数函数），我可以不是用点参数。然而，这种技术的语法支持在 R 中受到限制，这就是为什么会出现像 purrr 和 rlist 这样的软件包来试图改善这种情况；包 functional 和 pryr 提供专门的功能用于实现柯里化。

闭包（Closures）

R 中的每个函数都是一个闭包（除了原生函数）。闭包是一个具有与之相关环境的函数。例如，R 包中的函数可以访问包的名称空间，或在面向对象时类中的方法可以访问整个类。但是通常这个术语是在从其他函数中返回函数时使用的（每当你试图对闭包进行子集化时，除了 R 的错误消息外）。如果你对此还不了解，你可以阅读这篇文章或《Advanced R》的相关章节。

我的例子中，我使用闭包来为给定的 p 值重新定义平方根的不动点函数。我认为只有通过以下实施才能强调给定 p 值：

fpsqrt <- function(p)
{
    function(x) p / x
}

这实际上使算法的调用更加简洁：

fp(averageDamp(fpsqrt(2)), 2, converged)

## [1] 1.414214

缓存模式

在各种情况下，我都想缓存一些结果而不是重新计算它们。这是因为性能方面的考虑，因为无论使用哪个库，计算矩阵逆的时间关于样本大小都不是线性的。如果你有 10000 个观察值，并且要计算在蒙特卡罗模拟中得到的协方差矩阵（\(10000 \times 10000\)）的逆，你有得好等了。为了说明这一点，我设想计算一个线性估计量。尽管估计量可以通过解析方法来得到，但我使用了不动点算法。这种情况下，不动点函数中我使用 Newton-Raphson 算法，定义如下：

\[\beta_{n+1} = \beta_n - (f'' (\beta_n))^{-1} f'(\beta_n) \]

其中，

\[f'(\beta) = X^\top (y - X\beta) \\ f''(\beta) = -X^\top X \]

是正态分布情况下似然函数在 \(\beta\) 点的一阶和二阶导数。如果你对此没有什么概念，就把注意力完全放在这个模式上。在 R 中，我已经使用闭包应用了接口模式，请考虑以下实现：

nr <- function(X,
               y)
{
    function(beta) beta - solve(-crossprod(X)) %*% crossprod(X, y - X %*% beta)
}

# Some data to test the function:
set.seed(1)
X <- cbind(1, 1:10)
y <- X %*% c(1, 2) + rnorm(10)

# Average damping in this case will make the convergence a bit slower:
fp(printValue(averageDamp(nr(X, y))), c(1, 2), converged)

## 1 2 
## 0.9155882 2.027366 
## 0.8733823 2.041049 
## 0.8522794 2.047891 
## 0.8417279 2.051311 
## 0.8364522 2.053022 
## 0.8338143 2.053877 
## 0.8324954 2.054304

##           [,1]
## [1,] 0.8318359
## [2,] 2.0545183

# And to have a comparison:
stats::lm.fit(X, y)$coefficients

##        x1        x2 
## 0.8311764 2.0547321

结果看起来很好。我们应该选择一个不同的容忍度，以便更接近 R 的实现，目前看来这个容忍度选择得非常宽松，但这不是现在的重点。我现在想要做的是使 nr 的返回函数依赖于预先计算的值，以避免它们在每次迭代中重新计算。在这个例子中，我将它与局部函数的定义结合起来：

nr <- function(X, y)
{  
    # f1 relies on values in its scope:
    f1 <- function(beta) Xy - XX %*% beta

    Xy <- crossprod(X, y)
    XX <- crossprod(X)
    f2inv <- solve(-XX)

    function(beta) beta - f2inv %*% f1(beta)  
}

fp(averageDamp(nr(X, y)), c(1, 2), converged)

##           [,1]
## [1,] 0.8318359
## [2,] 2.0545183

一些评论：

在上面的例子中，像 f1 这样的局部函数定义是唯一依赖自由变量的地方。即 f1 依赖于封闭环境中定义的值 Xy 和 XX；这是我在上层函数中不惜一切代价要避免的。
我喜欢这种表示方法，因为我将不动点函数的逻辑保留在 nr；在给定数据的情况下，该函数知道如何计算下一次迭代的所有事。一种不同的方法是定义 nr，使其将 XX、Xy 和 f2inv 作为参数，这意味着我的代码的其他部分必须知道 nr 中的实现，并且我必须查看不同的地方以了解下一次迭代是如何进行的计算。

计数器模式

到目前为止，不动点框架不允许限制迭代次数。这当然是你总想要控制的东西。这次我使用闭包来模拟可变状态。考虑计数器的两种实现：

# Option 1:
counterConst <- function()
{
    # like a constructor function
    count <- 0
    function()
    {
        count <<- count + 1
        count
    }
}
counter <- counterConst()
counter()

## [1] 1

counter()

## [1] 2

# Option 2:
counter <- local({
    count <- 0
    function()
    {
        count <<- count + 1
        count
    }
})
counter()

## [1] 1

counter()

## [1] 2

我还记得闭包很难理解，因为在上面的例子中，当 R 中几乎所有东西都是不可变的时候，它们可以用来模拟可变状态。为从 Hadley （R 领域的知名专家）的一个例子中明白这一点，我可能用头撞了几个小时墙。

为了在不动点框架中实现最大迭代次数，我结合了包装器模式和计数器模式来修改收敛准则，以便算法在给定次数的迭代后终止：

addMaxIter <- function(converged,
                       maxIter)
{
    count <- 0
    function(...)
    {
        count <<- count + 1
        if (count >= maxIter) TRUE
        else converged(...)
    }
}

这使我们能够探索最初示例中发生的错误：

fp(printValue(fpsqrt(2)), 2, addMaxIter(converged, 4))

## 2 
## 1 
## 2 
## 1

## [1] 2

现在我们可以看到算法的初始版本在 1 到 2 之间振荡。您可能会说，在这种情况下，你还可以将迭代的次数看作算法的逻辑（作为 fp 的责任）。在这种情况下，我会争辩说，代码不再反映之前介绍的算法公式。但是让我们来比较一下不同的实现：

fpImp <- function(f,
                  x,
                  convCrit,
                  maxIter = 100)
{
    converged <- function()
    {
        convCrit(x, value) | count >= maxIter
    }
  
    count <- 0
    value <- NULL

    repeat {

        count <- count + 1 
        value <- f(x)

        if (converged()) break
        else
        {
            x <- value
            next
        }

    }

    list(result = value, iter = count)  
}

fpImp(averageDamp(fpsqrt(2)), 2, converged)

## $result
## [1] 1.414214
## 
## $iter
## [1] 4

现在让我们为 fp 的返回值添加迭代次数：

addIter <- function(fun)
{
    count <- 0
    function(x)
    {
        count <<- count + 1
        value <- fun(x)
        attr(value, "count") <- count
        value
    }
}

fp(addIter(averageDamp(fpsqrt(2))), 2, converged)

## [1] 1.414214
## attr(,"count")
## [1] 4

也许这个实现应该有一个自己的名字，但是你仍然可以看到包装器和计数器模式，它和 averageDamp 函数类似。与 fpImp 相比，围绕最大迭代次数的逻辑已经从算法的具体实现中分离出来。特别是如果我考虑添加更多功能，命令式的实现不可避免地必须应付越来越多的事情。相反，我可以在我的不动点框架中插入新功能。所以我认为对扩展开放对修改封闭，这不仅仅是一件好事，如果你喜欢面向对象的话。

计数器模式当然更一般化。它仅仅反映了模拟可变状态的一种策略。只有极少数情况下，我才真的需要这样做，计数是一个重复出现的主题。

函数柯里化（Currying）指的是把多个参数放进一个接受许多参数的函数，形成一个新的函数接受余下的参数，并返回结果 ↩︎

posted @ 2018-05-24 21:43 xuruilong100 阅读(676) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部