具有线性限制的优化解决方案答案

【问题标题】：Optimizing solutions with linear restrictions具有线性限制的优化解决方案
【发布时间】：2015-07-28 21:57:35
【问题描述】：

假设我有一个带有 N 变量的线性系统，但我只有 N-1 方程（约束）。如何使用 R 获得 N 个变量中每一个的可行集（范围）？

示例：

A <- matrix(data=c(0,1,0,1,0,1,0,1,
                   0,0,1,1,0,0,1,1,
                   0,0,0,0,1,1,1,1,
                   0,0,0,1,0,0,0,1,
                   0,0,0,0,0,1,0,1,
                   0,0,0,0,0,0,1,1,
                   1,1,1,1,1,1,1,1
                   ),
            ncol=8, byrow=T)
b <- matrix(data=c(0.2,0.4,0.6,0.06,0.18,0.12,1),
            ncol=1)
> A
##     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
##[1,]    0    1    0    1    0    1    0    1
##[2,]    0    0    1    1    0    0    1    1
##[3,]    0    0    0    0    1    1    1    1
##[4,]    0    0    0    1    0    0    0    1
##[5,]    0    0    0    0    0    1    0    1
##[6,]    0    0    0    0    0    0    1    1
##[7,]    1    1    1    1    1    1    1    1
> b
##     [,1]
##[1,] 0.20
##[2,] 0.40
##[3,] 0.60
##[4,] 0.06
##[5,] 0.18
##[6,] 0.12
##[7,] 1.00

附加约束：每个值都必须是正的，并且应该在 0 和 1 之间（可以省略最后一个，因为最后一个等式加起来只有正值 1）

【问题讨论】：

这类问题可以使用正则化回归等统计工具来解决。其中之一就是著名的 LASSO。

标签： r mathematical-optimization

【解决方案1】：

您所描述的问题可能最好通过线性规划来解决，因为您有一组线性约束（Ax = b，x >= 0，x lpSolve 包解决它：

library(lpSolve)
for (idx in seq_len(ncol(A))) {
  for (type in c("min", "max")) {
    mod <- lp(direction = type,
              objective.in = as.numeric(seq_len(ncol(A)) == idx),
              const.mat = rbind(A, diag(ncol(A))),
              const.dir = c(rep("=", length(b)), rep("<=", ncol(A))),
              const.rhs = c(b, rep(1, ncol(A))))
    print(paste("Variable:", idx, "type:", type, "value", mod$objval))
  }
}
# [1] "Variable: 1 type: min value 0.1"
# [1] "Variable: 1 type: max value 0.12"
# [1] "Variable: 2 type: min value 0"
# [1] "Variable: 2 type: max value 0.02"
# [1] "Variable: 3 type: min value 0.26"
# [1] "Variable: 3 type: max value 0.28"
# [1] "Variable: 4 type: min value 0"
# [1] "Variable: 4 type: max value 0.02"
# [1] "Variable: 5 type: min value 0.34"
# [1] "Variable: 5 type: max value 0.36"
# [1] "Variable: 6 type: min value 0.12"
# [1] "Variable: 6 type: max value 0.14"
# [1] "Variable: 7 type: min value 0.06"
# [1] "Variable: 7 type: max value 0.08"
# [1] "Variable: 8 type: min value 0.04"
# [1] "Variable: 8 type: max value 0.06"

【讨论】：

【解决方案2】：

您需要计算 A 的“零空间”，例如使用包 MASS 或 pracma 中的函数 nullspace()。首先计算一个最小二乘解，然后在零空间中添加向量的线性组合：

library(pracma)
N <- nullspace(A)
# 0.3535534 * c(-1, 1, 1, -1, 1, -1, 1, -1)
x0 <- qr.solve(A, b)
# [1]  0.16 -0.04  0.22  0.06  0.30  0.18  0.12  0.00

和 x0 + x*N, x real，生成所有可能的解决方案。

本例中的 N 为 1-dim。因为 A 的秩为 1。方程越少，零空间的维度就越多。

【讨论】：

有没有办法添加所有解决方案都必须是正数的限制？
没有 x 使得 x0 + x*N 的所有分量都是正数，您可以通过从所有分量中推导出 x 的不等式来验证（如果我计算正确的话）。

【解决方案3】：

我建议使用正则化回归方法。 R 中有包，如lars 或glmnet。我建议阅读此answer。

基本上，你会在R 中运行

library(lars)
model=lars(A,b)
coef(model)
      [,1] [,2]  [,3] [,4]   [,5]   [,6]    [,7] [,8]
[1,] 0.000 0.00 0.000    0 0.0000 0.0000 0.00000    0
[2,] 0.000 0.00 0.000    0 0.0982 0.0000 0.00000    0
[3,] 0.131 0.00 0.000    0 0.2000 0.0000 0.00000    0
[4,] 0.138 0.00 0.185    0 0.3849 0.0000 0.00000    0
[5,] 0.109 0.00 0.260    0 0.4000 0.0600 0.00000    0
[6,] 0.109 0.00 0.260    0 0.3997 0.0603 0.00029    0
[7,] 0.100 0.02 0.280    0 0.3600 0.1200 0.06000    0

结果是一个 7x8 矩阵，其中从上到下的每一行将最重要的变量输入到模型中。例如，第一行尝试仅使用截距拟合模型。第二个引入了 8 个中最强的变量，即第 5 个，依此类推。从结果中您可以看到变量 4 和 8 的系数一直为零，这表明它们在其他所有变量中并不重要。

【讨论】：