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

问题

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约束的困难上，例如:

Setting constraints in constrOptim

Constrain Optimisation Problems in R

我没有发现任何检查公差或对输出不满意的东西。

最佳答案

这不是一个确切的答案，但是比评论长，应该会有所帮助。

我认为您的问题有一个解析解决方案-很高兴知道您是否正在测试优化算法。

这是将预算固定为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)

选择正确的调整参数后，即使没有渐变，您也会得到单调递减的结果。

关于r - R中约束优化的非单调输出，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/23740812/