Как получить неотрицательные решения матрицы в R?

Question

Я пытаюсь решить уравнение AX = B в R.

У меня есть две матрицы А и В:

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, сводя к минимуму сумму квадратов ошибок и максимальную абсолютную ошибку.

Ошибка определяется как:

$ ошибок = 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. Поскольку ОП ищет конкретные неотрицательные целочисленные решения, мы можем установить ограничения на каждую переменную (т. Е. 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