如何获得R中矩阵的非负解？答案

【问题标题】：How to get non-negative solutions to matrix in R?如何获得R中矩阵的非负解？
【发布时间】：2016-07-27 20:20:00
【问题描述】：

我正在尝试在 R 中求解方程 AX = B。

我有两个矩阵，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 产生积极的结果，这就是我所追求的。

【问题讨论】：

我不明白。如果代数给你一个负值，那就是负数。没有办法强迫它是积极的。
@Roland 在许多情况下，线性方程组有无数个解。因此，有些解决方案具有负值，有些具有所有正值（在我的情况下更可取）。
我找到了使用nnls 包的解决方案。解是 (1900, 0, 3297, 0, 0, 4272, 10881, 0, 3782, 0, 2799, 0, 109, 2987, 7783, 0)。但是，我知道粗略地说X1 应该是~600 而X2 应该是~1300。无论如何在R中包含这个逻辑？（通过可能对X1（假设它必须>500）和X2（假设它必须>1200）施加限制？
@AhmedFasih 。不要在简历上交叉发帖。交叉发布强烈劝阻，并且是标记版主并在一个或两个站点上删除问题的理由。请参阅this link 进行讨论。
您对解决方案有什么更深入的了解吗？是否存在解决方案不太可能大于的最大值？

标签： r matrix linear-algebra

【解决方案1】：

您将需要一个误差函数来优化参数Xs。然后您可以使用许多包进行多变量优化。这里我使用了BB包中的两个函数，最小化平方误差和最大绝对误差。

错误定义为：

$ 错误 = 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

【讨论】：

【解决方案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 > 10^12) 解决方案，并且获得所有这些解决方案有点荒谬.这就是为什么在线性代数课程中，通常在每个变量上给出不等式的解决方案（例如330 <= x1 <= 738、1154 <= x2 < 1500 等（注意，这些不是实际的解决方案））。现在，我们可以通过明智地将某些值设置为零来轻松获得所需的解决方案。让我们尝试以下方法：

##   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

【讨论】：