constrOptim 函数：优化 R 代码

Question

我正在尝试解决优化问题。

下面是问题的数学解释和我使用的代码：

F = {f_1, f_2, ... f_n}

S = {s_1, s_2, ....s_m}

这里m总是大于n，sum(S)总是大于sum(F)

如果 ST = 转置(S)

求一个矩阵 P (n x m) = {p_ij}，使得：P %* % ST = F，其中 %* % 是矩阵乘法，关于以下约束：

1. p_ij>= 0，对于所有 i 和 j
2. sum (p_ij)

由于可能不存在确切的解决方案，我试图通过最小化 [ P %* % ST - F ].[ P %* % ST - F ] 来最小化错误，其中 .是点积

所以问题在于我使用以下代码的约束优化。

F = c(10,10,5)
S = c(8,8,9,8,4)

loss_fun <- function(P){
    T = matrix(S*P, nrow = n,ncol = m, byrow=T)
    F2 = rowSums(T) # Predicted values of F
    E = F - F2  # Error
    return(sum(E*E))    
}

n = length(F)
m = length(S)
P_init = c(rep(0.0001,n*m)) #Initial solution (theta)

# Creating constraint matrix
ui_1 = matrix(0,ncol = n*m, nrow= m) 
for (i in 1:m){
    for (j in 1:(n*m)) {
        if (i%%m==j%%m) ui_1[i,j] = -1
    }
}
ui_2 = diag(1,ncol = n*m, nrow = m*n)
my_ui <- rbind(ui_1,ui_2)

# Creating constraint vector
my_ci = c(rep(-1,m),rep(0,n*m))

z = constrOptim(P_init,loss_fun,NULL,ui=my_ui, ci=my_ci)

#result
P_final = matrix(z$par,nrow=n,byrow=T)

#verification of result
T = t(S*t(P_final)) #proportion matrix * S, transpose to ensure multiplication is by row.
F2 = rowSums(T) # Predicted values of F
E = F - F2  # Error
sum(E*E)

上面的代码运行良好，在我的机器上运行不到 0.5 秒，它有 i5 CPU、4 核、8 GB RAM、64 位 Windows 7 和 64 位 R 3.1.1。

但是，当我在实际问题中使用 F 和 S 时，它运行了大约 15 个小时而没有产生任何输出。 F 有 39 个元素，S 有 196 个。

F = c(212,359,186,396,460,449,206,180,383,264,294,179,256,294,173,415,363,323,389,219,298,338,287,434,195,450,120,460,164,395,198,108,72,345,54,450,420,488,262)
S = c(233.81,0,1.13,59.68,0,768.18,12.33,147.56,115.2,537.32,0,144.35,93.63,13.43,48.58,60,78.26,1280,369.62,8.11,0,342.96,452.99,521.72,4995.58,0,0,10.59,8.1,38.89,161.67,186.14,0,83.22,13.89,37.35,2370,0,0,8.61,4.95,6.31,0,1.53,3600,0,12.48,444.26,0,8490,615.25,27.11,402.95,393.46,1.26,0,44.36,728.85,37.61,159.06,103.63,145.38,0.51,0,0,18.6,3.24,44.5,17.46,210,128.03,19.48,340.79,54.79,54.42,48.48,0,44.76,0,0,0,43.19,102.03,0,0,470,0,101,0,9060,6.09,8.33,49.09,0,19.72,170,57.54,128.78,636.01,10.93,38.79,0,0,49.65,173.58,101.96,21.84,2.55,14.55,770,7419.13,216.21,238.15,582.95,57.93,26.97,71.88,4.63,0,31,103.37,570.58,45.79,540,348.9,151.82,207.41,29.56,51.73,92.25,0,0,51.39,25.14,0,0,95.21,298.94,5.77,154.29,280,1666.59,40.19,0,9.37,119.76,0,0,9.17,28.19,67.5,129.62,85.41,24.59,3607.98,0,130.28,99.57,0,0,0,36.23,1140,328.87,0,0,0,40,22.77,0,2.08,0,0,0,14.66,0,102.86,50.06,13.22,62.25,1410,860,930,646.15,0,0,0,0,890,0,0,12.61,86.4,95.35,19.31,87.74
)

rbind 本身需要 3 到 4 秒，但真正的问题是 constrOptim 所用的时间。

Answer 1

由于你的约束比较简单，当你使用一些可以将函数作为约束参数的包时，可以避免在约束部分进行大矩阵计算，例如alabama。

loss_fun <- function(P){
  T = matrix(S*P, nrow = n,ncol = m, byrow=T)
  F2 = rowSums(T) # Predicted values of F
  E = F - F2  # Error
  return(sum(E*E))    
}

n = length(F)
m = length(S)
P_init = c(rep(0.0001, n*m)) #Initial solution (theta)

# Creating inequality constraint function (this is much faster than my_ui %*% P - my_ci)
hin <- function(P){
  P_mat <- matrix(P, nrow = m)
  c(rowSums(P_mat) * -1 +1, P)
}

library(alabama)
aug_res <- auglag(P_init, loss_fun, hin = hin, control.outer = list(kkt2.check = FALSE))