【问题标题】:R: Quadratic programming/ Isotonic regressionR:二次规划/等渗回归
【发布时间】:2015-06-21 21:41:14
【问题描述】:

我想最小化以下等式:

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 包,因为这是一个二次规划问题。我不是在寻找一个完整的答案,我想要一些关于如何构建约束矩阵以及解决问题的最佳方法的指导。

【问题讨论】:

    标签: r mathematical-optimization minimization quadratic-programming


    【解决方案1】:

    鉴于此处的优化问题与your previous question 的相似之处,我将直接从我对该问题的回答中借用一些语言。但是它们有很大的不同(上一个问题是线性规划问题,这是一个二次规划问题,并且约束不同),所以它们不是重复的。

    扩展我们得到的优化目标Quw*ruw^2 - 2*Quw*ruw*yuw + Quw*yuw^2。我们看到这是决策变量yuw的二次函数,因此可以使用quadProg包的solve.QP方法来解决优化问题。

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

    ?solve.QP 中,我们了解到,对于决策变量b,目标采用-d'b + 0.5b'Db 的形式。 d对应决策变量yuw的元素有值2*Quw*ruwD是一个对角矩阵,对应决策变量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
    

    【讨论】:

    • 嗨 Josilber,一些 yuw 值与您沿行的方向相同,也与您沿列的方向相同。这是否符合 yuw >=yu,w+1 和 yu,w >=yu-1,w 约束?
    • 另外,第一个值 (20,1) 应该是 100 并且 y0,10 = 0 yu,10 = 0 约束不成立,请问为什么?
    • @Ankit 元素 (1, 1) 在左上角,所以yu,w &gt;= yu,w+1 表示当您从右向左移动时数字不会减少,yu,w &gt;= yu-1,w 表示数字不会减少当您从上到下移动时。我没有显示y20,0,因为它没有出现在目标值中。从y20,0 = 100yu,w &gt;= yu,w+1 我得出结论y20,1 &lt;= 100,这是我添加的约束。如您所见,y20,1 的最终选择值为 0.6298100。
    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2012-09-03
    • 2023-03-30
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多