约束优化

Question

我正在处理模型的输出，其中存在可能不遵循先验预期的参数估计值。我想编写一个函数来强制这些效用估计值符合这些预期。为此，该函数应最小化起始值和新估计值之间的平方偏差之和。由于我们有先验期望，因此优化应该受到以下约束：

B0 < B1
B1 < B2
...
Bj < Bj+1 

例如，下面的原始参数估计值是针对 B2 和 B3 翻转的。 Delta 和 Delta^2 列显示了原始参数估计值与新系数之间的偏差。我正在尝试最小化列Delta^2。我已经在 Excel 中对此进行了编码，并展示了 Excel 的 Solver 如何通过提供一组约束来优化此问题：

Beta    BetaRaw    Delta    Delta^2    BetaNew
B0       1.2       0        0          1.2
B1       1.3       0        0          1.3
B2       1.6       -0.2     0.04       1.4
B3       1.4       0        0          1.4
B4       2.2       0        0          2.2

阅读?optim 和?constrOptim 后，我无法理解如何在 R 中进行设置。我确定我只是有点密集，但可以在右侧使用一些指针方向！

3/24/2012 - 添加赏金，因为我不够聪明，无法翻译第一个答案。

这里有一些 R 代码应该在正确的路径上。假设 beta 开始于：

betas <- c(1.2,1.3,1.6,1.4,2.2)

我想最小化以下函数，使得b0 <= b1 <= b2 <= b3 <= b4

f <- function(x) {
  x1 <- x[1]
  x2 <- x[2]
  x3 <- x[3]
  x4 <- x[4]
  x5 <- x[5]

 loss <- (x1 - betas[1]) ^ 2 + 
         (x2 - betas[2]) ^ 2 + 
         (x3 - betas[3]) ^ 2 + 
         (x4 - betas[4]) ^ 2 +
         (x5 - betas[5]) ^ 2    

  return(loss)
}

为了证明该函数有效，如果我们将原始 beta 传递进来，损失应该为零：

> f(betas)
[1] 0

而且相对较大，有一些随机输入：

> set.seed(42)
> f(rnorm(5))
[1] 8.849329

并以我能够在 Excel 中计算的值最小化：

> f(c(1.2,1.3,1.4,1.4,2.2))
[1] 0.04

Answer 1

1. 由于目标是二次的并且约束是线性的， 你可以使用solve.QP。

它找到最小化的b

(1/2) * t(b) %*% Dmat %*% b - t(dvec) %*% b 

在约束下

t(Amat) %*% b >= bvec. 

在这里，我们希望b 最小化

sum( (b-betas)^2 ) = sum(b^2) - 2 * sum(b*betas) + sum(beta^2)
                   = t(b) %*% t(b) - 2 * t(b) %*% betas + sum(beta^2).

由于最后一个词，sum(beta^2)，是不变的，我们可以去掉它， 我们可以设置

Dmat = diag(n)
dvec = betas.

约束是

b[1] <= b[2]
b[2] <= b[3]
...
b[n-1] <= b[n]

即，

-b[1] + b[2]                       >= 0
      - b[2] + b[3]                >= 0
               ...
                   - b[n-1] + b[n] >= 0

所以t(Amat) 是

[ -1  1                ]
[    -1  1             ]
[       -1  1          ]
[             ...      ]
[                -1  1 ]

而bvec 为零。

这导致以下代码。

# Sample data
betas <- c(1.2, 1.3, 1.6, 1.4, 2.2)

# Optimization
n <- length(betas)
Dmat <- diag(n)
dvec <- betas
Amat <- matrix(0,nr=n,nc=n-1)
Amat[cbind(1:(n-1), 1:(n-1))] <- -1
Amat[cbind(2:n,     1:(n-1))] <-  1
t(Amat)  # Check that it looks as it should
bvec <- rep(0,n-1)
library(quadprog)
r <- solve.QP(Dmat, dvec, Amat, bvec)

# Check the result, graphically
plot(betas)
points(r$solution, pch=16)

2. 您可以以相同的方式使用constrOptim（目标函数可以是任意的，但约束必须是线性的）。

3. 更一般地，如果您重新参数化问题，您可以使用optim 变成一个无约束的优化问题， 比如

b[1] = exp(x[1])
b[2] = b[1] + exp(x[2])
...
b[n] = b[n-1] + exp(x[n-1]).

有几个例子 here 或there。

Answer 2

好的，这已开始形成，但仍有一些错误。根据与@Joran 聊天的对话，我似乎可以包含一个条件，如果值不按顺序设置，则将损失函数设置为任意大的值。如果前两个系数之间出现差异，这似乎有效，但此后则不然。我很难解释为什么会这样。

最小化函数：

f <- function(x, x0) {
  x1 <- x[1]
  x2 <- x[2]
  x3 <- x[3]
  x4 <- x[4]
  x5 <- x[5]

 loss <- (x1 - x0[1]) ^ 2 + 
         (x2 - x0[2]) ^ 2 + 
         (x3 - x0[3]) ^ 2 + 
         (x4 - x0[4]) ^ 2 +
         (x5 - x0[5]) ^ 2    

  #Make sure the coefficients are in order
  if any(diff(c(x1,x2,x3,x4,x5)) > 0) loss = 10000000

  return(loss)
}

工作示例（有点，如果b0 = 1.24，似乎损失会最小化？）：

> betas <- c(1.22, 1.24, 1.18, 1.12, 1.10)
> optim(betas, f, x0 = betas)$par
[1] 1.282 1.240 1.180 1.120 1.100

非工作示例（注意第三个元素仍然大于第二个：

> betas <- c(1.20, 1.15, 1.18, 1.12, 1.10)
> optim(betas, f, x0 = betas)$par
[1] 1.20 1.15 1.18 1.12 1.10