R：等渗回归最小化答案

【问题标题】：R: Isotonic regression MinimisationR：等渗回归最小化
【发布时间】：2015-06-20 00:01:45
【问题描述】：

我想最小化以下等式：

F=SUM{u 1:20}sum{w 1:10}   Quw(ruw-yuw)

具有以下约束：

yuw >= yu,w+1
yuw >= yu-1,w
y20,0 >= 100
y0,10 >= 0

我有一个 20*10 ruw 和 20*10 quw 矩阵，我现在需要生成一个符合约束的 yuw 矩阵。我正在用 R 编码并且熟悉 lpsolve 和 optimx 包，但不知道如何将它们用于这个特定问题。

【问题讨论】：

在我的解决方案中，我表明最佳解决方案不依赖于 ruw 值。您是否在配方中遗漏了某些内容，例如目标中的绝对值？
最后两个约束我说错了，最后两个约束实际上应该是 y20,0 = 100 和 y0,10 = 0。
Josilber，(ruw-yuw) 的平方应该是：F=SUM{u 1:20}sum{w 1:10} Quw(ruw-yuw)^2。这将如何影响设置？仅提供有关设置的指导将不胜感激。
quadprog 软件包有帮助吗？
不幸的是，这个错字改变了整个问题，因为它现在是一个二次程序而不是一个线性程序。鉴于这种差异有多大，我建议您只提出一个新问题，确保这次您的表述正确。

标签： r regression mathematical-optimization linear-programming minimization

【解决方案1】：

因为Quw 和ruw 都是数据，所以所有约束和目标在yuw 决策变量中都是线性的。因此，这是一个线性规划问题，可以通过 lpSolve 包解决。

为了稍微抽象一下，假设R=20 和C=10 描述了输入矩阵的维度。然后是R*C决策变量，我们可以为它们分配顺序y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC，读取变量矩阵的列。

目标中每个yuw变量的系数为-Quw；请注意，求和中的Quw*ruw 项是常数（也不受我们为决策变量选择的值的影响），因此不会输入到线性规划求解器。有趣的是，这意味着ruw实际上对优化模型解没有任何影响。

第一个R*(C-1) 约束对应于yuw >= yu,w+1 约束，接下来的(R-1)*C 约束对应于yuw >= yu-1,w 约束。最后两个约束对应于y20,1 >= 100 和y1,10 >= 0 约束。

我们可以使用以下 R 命令将该模型输入到 lpsolve 包中，将一个简单的 Q 矩阵作为输入，其中每个条目都是 -1（得到的解决方案应该将所有决策变量设置为 0，除了左下角, 应该是 100):

# Sample data
Quw <- matrix(-1, nrow=20, ncol=10)
R <- nrow(Quw)
C <- ncol(Quw)

# Build constraint matrix
part1 <- matrix(0, nrow=R*(C-1), ncol=R*C)
part1[cbind(1:(R*C-R), 1:(R*C-R))] <- 1
part1[cbind(1:(R*C-R), (R+1):(R*C))] <- -1
part2 <- matrix(0, nrow=(R-1)*C, ncol=R*C)
pos2 <- as.vector(sapply(2:R, function(r) r+seq(0, R*(C-1), by=R)))
part2[cbind(1:nrow(part2), pos2)] <- 1
part2[cbind(1:nrow(part2), pos2-1)] <- -1
part3 <- rep(0, R*C)
part3[R] <- 1
part4 <- rep(0, R*C)
part4[(C-1)*R + 1] <- 1
const.mat <- rbind(part1, part2, part3, part4)

library(lpSolve)
mod <- lp(direction = "min",
          objective.in = as.vector(-Quw),
          const.mat = const.mat,
          const.dir = rep(">=", nrow(const.mat)),
          const.rhs = c(rep(0, nrow(const.mat)-2), 100, 0))

我们现在可以访问模型解决方案：

mod
# Success: the objective function is 100
matrix(mod$solution, nrow=R)
#       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#  [1,]    0    0    0    0    0    0    0    0    0     0
#  [2,]    0    0    0    0    0    0    0    0    0     0
#  [3,]    0    0    0    0    0    0    0    0    0     0
#  [4,]    0    0    0    0    0    0    0    0    0     0
#  [5,]    0    0    0    0    0    0    0    0    0     0
#  [6,]    0    0    0    0    0    0    0    0    0     0
#  [7,]    0    0    0    0    0    0    0    0    0     0
#  [8,]    0    0    0    0    0    0    0    0    0     0
#  [9,]    0    0    0    0    0    0    0    0    0     0
# [10,]    0    0    0    0    0    0    0    0    0     0
# [11,]    0    0    0    0    0    0    0    0    0     0
# [12,]    0    0    0    0    0    0    0    0    0     0
# [13,]    0    0    0    0    0    0    0    0    0     0
# [14,]    0    0    0    0    0    0    0    0    0     0
# [15,]    0    0    0    0    0    0    0    0    0     0
# [16,]    0    0    0    0    0    0    0    0    0     0
# [17,]    0    0    0    0    0    0    0    0    0     0
# [18,]    0    0    0    0    0    0    0    0    0     0
# [19,]    0    0    0    0    0    0    0    0    0     0
# [20,]  100    0    0    0    0    0    0    0    0     0

请注意，如果 Quw 更改（例如，如果我们用 1 而不是 -1 填充它），您的模型很容易变得不可行。在这些情况下，模型将以状态 3 退出（您可以通过运行 mod 或 mod$status 看到这一点）。

【讨论】：