如何在R中得到非負解的矩陣？

Question

我試圖求解方程AX = B在R.

我有兩個矩陣A和B：

A = matrix(c(1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 
     0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0, 
     0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0, 
     0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0, 
     0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0, 
     0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0, 
     0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1, 
     0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0, 
     0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1, 
     1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0), byrow = T, nrow = 10, ncol = 16) 

B = matrix(c(1900,2799,3096,3297,3782,4272,7783,10881,7259,30551), nrow = 10, ncol = 1) 

我的問題是，如何解決AX = B和保證非負解？我正在解決的數值（X1, X2,...X15, X16）是人口數字，所以他們不能是負數。理想情況下，它們也是整數值，但一次只能做一件事。

有沒有一種簡單的方法來做到這一點在R？

我找到了一種方法來做到這一點here，但它並沒有產生積極的結果，所有X這是我所追求的。

Answer 1

您將需要一個錯誤函數來優化參數Xs。然後，您可以使用多個包進行多變量優化。這裏我使用BB包中的兩個函數，最小化平方誤差的總和以及最大絕對誤差。

錯誤定義爲：

$ error = A。 X - B $

您可以將lower和upper參數定義爲常量或向量來添加框約束。

set.seed(321) 
p0 <- abs(runif(16))*1000 #same starting values 
library(BB) 
se2 <- function(x){ # sum of squared errors 
    sum(1000 * (A %*% x - B)^2) 
} 
se2(p0) 

seab <- function(x){ # max error 
    max(1000 * abs(A %*% x - B)) 
} 
seab(p0) 

s1<-spg(par=p0, fn=se2,lower=0) 
s2<-BBoptim(par=p0, fn=se2,lower=0) 
s3<-spg(par=p0, fn=seab,lower=0) 
s4<-BBoptim(par=p0, fn=seab,lower=0) 

linf=c(500,1200,rep(0,length(p0)-2)) #prior inferior 
s5<-spg(par=p0, fn=se2,lower=linf) 
s6<-BBoptim(par=p0, fn=se2,lower=linf) 

round(cbind(s1$par,s2$par,s3$par,s4$par,s5$par,s6$par),2) # Xs 

round(c(se2(s1$par),se2(s2$par),seab(s3$par),seab(s4$par),se2(s5$par),se2(s6$par)),3) #functions 

結果：

> round(cbind(s1$par,s2$par,s3$par,s4$par,s5$par,s6$par),2) # Xs 
     [,1] [,2] [,3] [,4] [,5] [,6] 
[1,] 1900.00 1900.00 1334.61 1898.22 700.00 700.00 
[2,] 0.00 0.00 565.39 1.78 1200.00 1200.00 
[3,] 3166.80 3166.80 3211.66 3177.49 3297.00 3297.00 
[4,] 130.20 130.20 85.34 119.51 0.00 0.00 
[5,] 3687.42 3687.42 3839.29 3966.84 3954.88 3954.87 
[6,] 584.58 584.58 432.71 305.16 317.12 317.13 
[7,] 7048.48 7048.48 7200.57 6852.74 7315.90 7315.93 
[8,] 3832.52 3832.52 3680.43 4028.26 3565.10 3565.07 
[9,] 3239.80 3239.80 3356.15 3509.46 3507.25 3507.24 
[10,] 542.20 542.20 425.85 272.54 274.75 274.76 
[11,] 2799.00 2799.00 2788.24 2796.45 2799.00 2799.00 
[12,] 0.00 0.00 10.76 2.55 0.00 0.00 
[13,] 3096.00 3096.00 3054.95 2954.85 3096.00 3096.00 
[14,] 0.00 0.00 41.05 141.15 0.00 0.00 
[15,] 5613.50 5613.50 5765.53 5394.95 5880.97 5880.97 
[16,] 2169.50 2169.50 2017.47 2388.05 1902.03 1902.03 
> 
> round(c(se2(s1$par),se2(s2$par),seab(s3$par),seab(s4$par),se2(s5$par),se2(s6$par)),3) #functions 
[1] 0.000 0.000 1.161 0.000 0.000 0.000 

Answer 2

我不知道會帶來你想要的結果的任何包，所以我寫了下面的一個基礎R解決方案的矩陣降低到row echelon form。之後，我們需要做一個小代數來獲得結果。

reduceMatrixSO <- function(mat) { 
    n1 <- ncol(mat); n2 <- nrow(mat) 
    mymax <- 1L 
    for (i in 1:(n1-1L)) { 
     temp <- which(mat[,i] != 0L) 
     t <- which(temp >= mymax) 
     if (length(temp)>0L && length(t)>0L) { 
      MyMin <- min(temp[t]) 
      if (!(MyMin==mymax)) { 
       vec <- mat[MyMin,] 
       mat[MyMin,] <- mat[mymax,] 
       mat[mymax,] <- vec 
      } 
      Coef1 <- mat[mymax, i] 
      t <- t[-1]; temp <- temp[t] 
      for (j in temp) { 
       Coef2 <- mat[j,i] 
       mat[j,] <- mat[j,] - mat[mymax,]*Coef2/Coef1 
      } 
      mymax <- mymax+1L 
     } 
    } 

    if (mymax<n2) {simpMat <- mat[-(mymax:n2),]} else {simpMat <- mat} 
    lenSimp <- nrow(simpMat) 
    if (is.null(lenSimp)) {lenSimp <- 0L} 
    mycols <- 1:n1 

    if (lenSimp>1L) { 
     ## "Diagonalizing" Matrix 
     for (i in 1:lenSimp) { 
      if (all(simpMat[i,]==0L)) {simpMat <- simpMat[-i,]; next} 
      if (simpMat[i,i]==0L) { 
       t <- min(which(simpMat[i,] != 0L)) 
       vec <- simpMat[,i]; tempCol <- mycols[i] 
       simpMat[,i] <- simpMat[,t]; mycols[i] <- mycols[t] 
       simpMat[,t] <- vec; mycols[t] <- tempCol 
      } 
     } 

     colnames(simpMat) <- c(sapply(mycols[-n1], function(r) paste(c("x",r),collapse = "")),"B") 
     lenSimp <- nrow(simpMat) 
     MyFree <- sapply(mycols[which((1:(n1-1L))>lenSimp)], function(r) paste(c("x",r),collapse = "")) 

     for (i in 1:lenSimp) { 
      temp <- which(simpMat[,i] != 0L) 
      t <- which(temp != i) 
      if (length(temp)>0L && length(t)>0L) { 
       Coef1 <- simpMat[i,i] 
       temp <- temp[t] 
       for (j in temp) { 
        Coef2 <- simpMat[j,i] 
        simpMat[j,] <- simpMat[j,] - simpMat[i,]*Coef2/Coef1 
       } 
      } 
     } 
     list(ReducedMatrix = simpMat, FreeVariables = MyFree) 
    } else { 
     list(NULL,NULL) 
    } 
} 

上面的代碼是對眼睛有點粗糙，但它是很簡單的和非常快的（我已經測試了更大的矩陣，並立即返回正確的結果）。下面是它給出的結果：

NewM <- cbind(A,B) 

reduceMatrixSO(NewM) 
$ReducedMatrix 
    x1 x2 x3 x5 x7 x9 x11 x13 x15 x10 x4 x12 x8 x14 x6 x16  B 
[1,] 1 0 0 0 0 0 0 0 0 -1 -1 -1 -1 -1 -1 -1 -5359 
[2,] 0 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 7259 
[3,] 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 3297 
[4,] 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 4272 
[5,] 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 10881 
[6,] 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 3782 
[7,] 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 2799 
[8,] 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 3096 
[9,] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 7783 

$FreeVariables 
[1] "x10" "x4" "x12" "x8" "x14" "x6" "x16" 

這告訴我們，有6 free variables。由於OP正在尋找特定的非負整數解，我們可以對每個變量施加約束（即x1> 500和x2 < 1200）。這可以通過很多方式來完成，雖然一般的解決方案（有這樣的約束）不包含無限多的解決方案，但有很多（>> 10^12）解決方案，並且獲得所有這些解決方案將會有點無意義。這就是爲什麼在線性代數課程中，解決方案通常在每個變量（例如330 <= x1 <= 738,1154 <= x2 < 1500等（N.B.這些不是實際的解決方案））上有不等式。現在，我們可以通過明智地將某些值設置爲零來輕鬆獲得所需的解決方案。讓我們嘗試以下操作：

## x10 = x4 = x12 = x8 = x14 = x6 = 0 
## x1 >= 500 && x2 <= 1200 

x1 <- -5359 + x16 ## ==>> -5359 + x16 >= 500 ==>> x16 >= 5859 
x2 <- 7259 - x16 ## ==>> 7259 - x16 >= 1200 ==>> x16 <= 6059 
x3 <- 3297 
x5 <- 4272 
x7 <- 10881 
x9 <- 3782 
x11 <- 2799 
x13 <- 3096 
x15 <- 7783 - x16 ## ==>> x16 <= 7783 

Ultimately 5859 <= x16 <= 6059 

讓我們嘗試了上述解決方案：

set.seed(5467) 
x16 <- sample(5859:6059, 1) 
## set variables above using the newly defined x16 
X <- c(x1,x2,x3,0,x5,0,x7,0,x9,0,x11,0,x13,0,x15,x16) 

all(A %*% X == B) 
[1] TRUE 

all(X >= 0) ## all non-negative 
[1] TRUE 

all(gmp::is.whole(X)) ## all integers 
[1] TRUE 

最後，我們確定的變量和打印解決方案：

names(X) <- sapply(1:16, function(r) paste(c("x",r), collapse = "")) 

X 
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 
590 1310 3297  0 4272  0 10881  0 3782  0 2799  0 3096  0 1834 5949