R中约束优化的非单调输出

Question

问题

constOptim函数是R给我一组参数估计.这些参数估计值是一年中12个不同点的消费值,应该是单调递减的.

我需要它们是单调的,并且每个参数之间的间隙都适合我想到的应用程序.出于此目的,支出值中的模式很重要,而不是绝对值.我想在优化术语中,这意味着我需要与参数估计值的差异相比较小的容差.

最小工作示例(具有简单的实用功能)

# Initial Parameters and Functions
Budget          = 1
NumberOfPeriods = 12
rho = 0.996
Utility_Function <- function(x){ x^0.5 }
Time_Array = seq(0,NumberOfPeriods-1)

# Value Function at start of time.
ValueFunctionAtTime1   = function(X){
  Frame                = data.frame(X, time = Time_Array)
  Frame$Util           = apply(Frame, 1, function(Frame) Utility_Function(Frame["X"]))
  Frame$DiscountedUtil = apply(Frame, 1, function(Frame) Frame["Util"] * rho^(Frame["time"])) 
  return(sum(Frame$DiscountedUtil))
}

# The sum of all spending in the year should be less than than the annual budget.
# This gives the ui and ci arguments
Sum_Of_Annual_Spends   = c(rep(-1,NumberOfPeriods))

# The starting values for optimisation is an equal expenditure in each period.
# The denominator is multiplied by 1.1 to avoid an initial values out of range error.
InitialGuesses         = rep(Budget/(NumberOfPeriods*1.1), NumberOfPeriods)

# Optimisation
Optimal_Spending    = constrOptim(InitialGuesses,
                                  function(X) -ValueFunctionAtTime1(X),
                                  NULL,
                                  ui = Sum_Of_Annual_Spends,
                                  ci = -Budget,
                                  outer.iterations = 100,
                                  outer.eps = 1e-10)$par

结果:

函数的输出不是单调的.

plot( Time_Array , Optimal_Spending)

我尝试修复它

我试过了:

• 增加容差(这在代码中是上面的outer.eps = 1e-10)
• 增加迭代次数(这在上面的代码中outer.iterations = 100)
• 提高初始参数值的质量.我用我的实际案例做了这个(相同但具有更复杂的实用功能),但没有解决问题.
• 通过增加预算或将效用函数乘以标量来扩展问题.

关于constOptim的其他问题

其他SO问题集中在编写constOptim约束的困难,例如:

1. 在constrOptim中设置约束
2. R中的约束优化问题

我没有发现任何检查容差或对输出不满意的事情.

Answer 1

这不完全是一个答案，但它比评论长，应该会有所帮助。

我认为你的问题有一个分析解决方案 - 如果你正在测试优化算法，那么知道它是很好的。

这是预算固定为 1.0 时的情况。

analytical.solution <- function(rho=0.9, T=10) {
    sapply(seq_len(T) - 1, function(t) (rho ^ (2*t)) * (1 - rho^2) / (1 - rho^(2*T)))
}
sum(analytical.solution())  # Should be 1.0, i.e. the budget

此处，消费者在 {0, 1, ..., T-1} 期间消费。该解确实随着时间索引单调递减。我通过设置拉格朗日并处理一阶条件得到了这个。

编辑：

我重写了您的代码，一切似乎都正常工作： constrOptim 给出了与我的分析解决方案一致的解决方案。预算固定为1。

analytical.solution <- function(rho=0.9, T=10) {
    sapply(seq_len(T) - 1, function(t) (rho ^ (2*t)) * (1 - rho^2) / (1 - rho^(2*T)))
}
candidate.solution <- analytical.solution()
sum(candidate.solution)  # Should be 1.0, i.e. the budget

objfn <- function(x, rho=0.9, T=10) {
    stopifnot(length(x) == T)
    sum(sqrt(x) * rho ^ (seq_len(T) - 1))
}
objfn.grad <- function(x, rho=0.9, T=10) {
    rho ^ (seq_len(T) - 1) * 0.5 * (1/sqrt(x))
}

## Sanity check the gradient
library(numDeriv)
all.equal(grad(objfn, candidate.solution), objfn.grad(candidate.solution))  # True

ui <- rbind(matrix(data=-1, nrow=1, ncol=10), diag(10))  # First row: budget constraint; other rows: x >= 0
ci <- c(-1, rep(10^-8, 10))
all(ui %*% candidate.solution - ci >= 0)  # True, the candidate solution is admissible
result1 <- constrOptim(theta=rep(0.01, 10), f=objfn, ui=ui, ci=ci, grad=objfn.grad, control=list(fnscale=-1))
round(abs(result1$par - candidate.solution), 4)  # Essentially zero
result2 <- constrOptim(theta=candidate.solution, f=objfn, ui=ui, ci=ci, grad=objfn.grad, control=list(fnscale=-1))
round(abs(result2$par - candidate.solution), 4)  # Essentially zero

关于渐变的后续：

即使 grad=NULL ，优化似乎也有效，这意味着您的代码中可能存在错误。看这个：

result3 <- constrOptim(theta=rep(0.01, 10), f=objfn, ui=ui, ci=ci, grad=NULL, control=list(fnscale=-1))
round(abs(result3$par - candidate.solution), 4)  # Still very close to zero
result4 <- constrOptim(theta=c(10^-6, 1-10*10^-6, rep(10^-6, 8)), f=objfn, ui=ui, ci=ci, grad=NULL, control=list(fnscale=-1))
round(abs(result4$par - candidate.solution), 4)  # Still very close to zero

关于rho=0.996情况的后续：

当 rho->1 时，解应该收敛到rep(1/T, T)——这解释了为什么即使 constrOptim 产生的小错误也会对输出是否单调递减产生显着影响。

当 rho=0.996 时，调整参数似乎会影响 constrOptim 的输出，足以改变单调性 - 见下文：

candidate.solution <- analytical.solution(rho=0.996)
candidate.solution  # Should be close to rep(1/10, 10) as discount factor is close to 1.0

result5 <- constrOptim(theta=c(10^-6, 1-10*10^-6, rep(10^-6, 8)), f=objfn,
                       ui=ui, ci=ci, grad=objfn.grad, control=list(fnscale=-1), rho=0.996)
round(abs(result5$par - candidate.solution), 4)
plot(result5$par)  # Looks nice when we used objfn.grad, as you pointed out

play.with.tuning.parameter <- function(mu) {
    result <- constrOptim(theta=c(10^-6, 1-10*10^-6, rep(10^-6, 8)), f=objfn,
                          mu=mu, outer.iterations=200, outer.eps = 1e-08,
                          ui=ui, ci=ci, grad=NULL, control=list(fnscale=-1), rho=0.996)
    return(mean(diff(result$par) < 0))
}

candidate.mus <- seq(0.01, 1, 0.01)
fraction.decreasing <- sapply(candidate.mus, play.with.tuning.parameter)
candidate.mus[fraction.decreasing == max(fraction.decreasing)]  # A few little clusters at 1.0
plot(candidate.mus, fraction.decreasing)  # ...but very noisy

result6 <- constrOptim(theta=c(10^-6, 1-10*10^-6, rep(10^-6, 8)), f=objfn,
                          mu=candidate.mus[which.max(fraction.decreasing)], outer.iterations=200, outer.eps = 1e-08,
                          ui=ui, ci=ci, grad=NULL, control=list(fnscale=-1), rho=0.996)
plot(result6$par)
round(abs(result6$par - candidate.solution), 4)

当您选择正确的调整参数时，即使没有梯度，您也会得到单调递减的结果。