R:二次规划/等渗回归
R: Quadratic programming/ Isotonic regression
我想最小化以下等式:
F=SUM{u 1:20}sum{w 1:10} Quw(ruw-yuw)^2
具有以下限制条件:
yuw >= yu,w+1
yuw >= yu-1,w
y20,0 = 100
y0,10 = 0
yu,10 = 0
我有一个 20*10 ruw 和 20*10 quw 矩阵,我现在需要生成一个遵守约束的 yuw 矩阵。我在 R 中编码并且熟悉 lpsolve、optimx 和 quadprog 包,但不知道如何将它们用于这个特定问题。我知道我必须使用 quadprog 包,因为这是一个二次规划问题。我不是在寻找完整的答案,我想要一些关于如何构造约束矩阵和解决问题的最佳方法的指导。
考虑到这里的优化问题与 的相似之处,我将直接从我对该问题的回答中借用一些语言。不过还是有点不同(之前的问题是线性规划问题,这次是二次规划问题,约束不一样),所以不是重复的。
扩展优化 objective 我们得到 Quw*ruw^2 - 2*Quw*ruw*yuw + Quw*yuw^2
。我们看到这是决策变量yuw
的二次函数,因此可以使用quadProg
包的solve.QP
方法来解决优化问题。
为了抽象一点,我们假设 R=20
和 C=10
描述了输入矩阵的维度。然后有 R*C
个决策变量,我们可以为它们分配顺序 y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC
,读取变量矩阵的列。
从 ?solve.QP
中,我们了解到 objective 对于决策变量 b
采用 -d'b + 0.5b'Db
的形式。 d
对应决策变量yuw
的元素取值2*Quw*ruw
,D
是一个对角矩阵,决策变量yuw
对应的元素取值2*Quw
。请注意,solve.QP
函数要求 D
矩阵为正定矩阵,因此我们要求每个 u, w
对的 Quw > 0
。
前R*(C-1)
个约束对应yuw >= yu,w+1
个约束,接下来的(R-1)*C
个约束对应yuw >= yu-1,w
个约束。接下来的2*R
约束对应于yuC = 0
约束(输入为yuC >= 0
和-yuC >= 0
),最后一个约束为-yR1 >= -100
(逻辑上等价) yR0 = 100
).
我们可以使用以下 R 命令将此模型输入到 quadProg
包中,使用随机输入数据:
# Sample data
set.seed(144)
Quw <- matrix(runif(200), nrow=20)
ruw <- matrix(runif(200), nrow=20)
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 <- matrix(0, nrow=2*R, ncol=R*C)
part3[cbind(1:R, (R*C-R+1):(R*C))] <- 1
part3[cbind((R+1):(2*R), (R*C-R+1):(R*C))] <- -1
part4 <- rep(0, R*C)
part4[R] <- -1
const.mat <- rbind(part1, part2, part3, part4)
library(quadProg)
mod <- solve.QP(Dmat = 2*diag(as.vector(Quw)),
dvec = 2*as.vector(ruw)*as.vector(Quw),
Amat = t(const.mat),
bvec = c(rep(0, nrow(const.mat)-1), -100))
我们现在可以访问模型解决方案:
# Objective (including the constant term):
mod$value + sum(Quw*ruw^2)
# [1] 9.14478
matrix(mod$solution, nrow=R)
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.3215992 0.1818095 0.1818095 0.1818095 0.000000e+00
# [2,] 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [3,] 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 2.775558e-17
# [4,] 0.5728478 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [5,] 0.5728478 0.5111456 0.5111456 0.4699046 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [6,] 0.5728478 0.5111456 0.5111456 0.4699046 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [7,] 0.5728478 0.5111456 0.5111456 0.4699046 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [8,] 0.5728478 0.5111456 0.5111456 0.5111456 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [9,] 0.5728478 0.5111456 0.5111456 0.5111456 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 1.110223e-16
# [10,] 0.5728478 0.5111456 0.5111456 0.5111456 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [11,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [12,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [13,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [14,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [15,] 0.6298100 0.6009718 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [16,] 0.6298100 0.6009718 0.6009718 0.6009718 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [17,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [18,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.6009718 0.4346995 0.4346995 0.4346995 0.000000e+00
# [19,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.6009718 0.4346995 0.4346995 0.4346995 0.000000e+00
# [20,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.6009718 0.5643033 0.5643033 0.5643033 0.000000e+00
我想最小化以下等式:
F=SUM{u 1:20}sum{w 1:10} Quw(ruw-yuw)^2
具有以下限制条件:
yuw >= yu,w+1
yuw >= yu-1,w
y20,0 = 100
y0,10 = 0
yu,10 = 0
我有一个 20*10 ruw 和 20*10 quw 矩阵,我现在需要生成一个遵守约束的 yuw 矩阵。我在 R 中编码并且熟悉 lpsolve、optimx 和 quadprog 包,但不知道如何将它们用于这个特定问题。我知道我必须使用 quadprog 包,因为这是一个二次规划问题。我不是在寻找完整的答案,我想要一些关于如何构造约束矩阵和解决问题的最佳方法的指导。
考虑到这里的优化问题与
扩展优化 objective 我们得到 Quw*ruw^2 - 2*Quw*ruw*yuw + Quw*yuw^2
。我们看到这是决策变量yuw
的二次函数,因此可以使用quadProg
包的solve.QP
方法来解决优化问题。
为了抽象一点,我们假设 R=20
和 C=10
描述了输入矩阵的维度。然后有 R*C
个决策变量,我们可以为它们分配顺序 y11, y21, ... yR1, y12, y22, ... yR2, ..., y1C, y2C, ..., yRC
,读取变量矩阵的列。
从 ?solve.QP
中,我们了解到 objective 对于决策变量 b
采用 -d'b + 0.5b'Db
的形式。 d
对应决策变量yuw
的元素取值2*Quw*ruw
,D
是一个对角矩阵,决策变量yuw
对应的元素取值2*Quw
。请注意,solve.QP
函数要求 D
矩阵为正定矩阵,因此我们要求每个 u, w
对的 Quw > 0
。
前R*(C-1)
个约束对应yuw >= yu,w+1
个约束,接下来的(R-1)*C
个约束对应yuw >= yu-1,w
个约束。接下来的2*R
约束对应于yuC = 0
约束(输入为yuC >= 0
和-yuC >= 0
),最后一个约束为-yR1 >= -100
(逻辑上等价) yR0 = 100
).
我们可以使用以下 R 命令将此模型输入到 quadProg
包中,使用随机输入数据:
# Sample data
set.seed(144)
Quw <- matrix(runif(200), nrow=20)
ruw <- matrix(runif(200), nrow=20)
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 <- matrix(0, nrow=2*R, ncol=R*C)
part3[cbind(1:R, (R*C-R+1):(R*C))] <- 1
part3[cbind((R+1):(2*R), (R*C-R+1):(R*C))] <- -1
part4 <- rep(0, R*C)
part4[R] <- -1
const.mat <- rbind(part1, part2, part3, part4)
library(quadProg)
mod <- solve.QP(Dmat = 2*diag(as.vector(Quw)),
dvec = 2*as.vector(ruw)*as.vector(Quw),
Amat = t(const.mat),
bvec = c(rep(0, nrow(const.mat)-1), -100))
我们现在可以访问模型解决方案:
# Objective (including the constant term):
mod$value + sum(Quw*ruw^2)
# [1] 9.14478
matrix(mod$solution, nrow=R)
# [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
# [1,] 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.3215992 0.1818095 0.1818095 0.1818095 0.000000e+00
# [2,] 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [3,] 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 2.775558e-17
# [4,] 0.5728478 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [5,] 0.5728478 0.5111456 0.5111456 0.4699046 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [6,] 0.5728478 0.5111456 0.5111456 0.4699046 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [7,] 0.5728478 0.5111456 0.5111456 0.4699046 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [8,] 0.5728478 0.5111456 0.5111456 0.5111456 0.4346995 0.4346995 0.4346995 0.2882339 0.2882339 0.000000e+00
# [9,] 0.5728478 0.5111456 0.5111456 0.5111456 0.4346995 0.4346995 0.4346995 0.4346995 0.2882339 1.110223e-16
# [10,] 0.5728478 0.5111456 0.5111456 0.5111456 0.4346995 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [11,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [12,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [13,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [14,] 0.6298100 0.5111456 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [15,] 0.6298100 0.6009718 0.5111456 0.5111456 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [16,] 0.6298100 0.6009718 0.6009718 0.6009718 0.4518205 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [17,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.4346995 0.4346995 0.4346995 0.4346995 0.000000e+00
# [18,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.6009718 0.4346995 0.4346995 0.4346995 0.000000e+00
# [19,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.6009718 0.4346995 0.4346995 0.4346995 0.000000e+00
# [20,] 0.6298100 0.6009718 0.6009718 0.6009718 0.6009718 0.6009718 0.5643033 0.5643033 0.5643033 0.000000e+00