我正试图求解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中有简单的方法吗？
我在这里找到了一种方法，但它并不能为所有X产生积极的结果，这就是我所追求的。
您将需要一个错误函数来优化参数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
我不知道有任何包会带来你想要的结果，所以我在下面写了一个基本的R解决方案，将矩阵简化为行-梯队形式。之后，我们需要做一点代数来获得结果。
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个自由变量。由于OP正在寻找特定的非负整数解，我们可以对每个变量施加约束（即x1>500和x2<1200）。这可以通过多种方式实现，虽然一般解决方案（具有此类约束）不包含无限多的解决方案，但有许多（>>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