在R中解决不定方程组

Question

我有一个方程组，希望使用数值方法来解决它。 我希望在给定初始种子的情况下得到一个接近的解。 让我解释一下。

我有一个常量向量X，其中包含以下值：

X <- (c(1,-2,3,4))

还有一个权重向量W：

W <- (c(0.25,0.25,0.25,0.25))

我希望W的所有元素之和等于1 (sum(W)=1)，并且X和W逐个元素相乘的结果之和为给定的数字N (sum(W*X)=N)。

在R中有简单的方法吗？我已经在Excel中使用了Solver，但我需要自动化它。

- arodrisa

Answer 1

这是您的常量和目标值：

x <- c(1, -2, 3, 4)
n <- 10

您需要一个最小化函数。第一行包含每个条件，第二行提供了如何将错误组合成单个分数的度量。您可能想要更改第二行。例如，您可以使用sum(c(1, 5) * errs ^ 2)使一个错误项比另一个错误项更加重要。

fn <- function(w)
{
  errs <- c(sum(w) - 1, sum(x * w) - n) 
  sum(errs ^ 2)
}

最简单的方法是从所有权重都相同的值开始。

init_w <- rep.int(1 / length(x), length(x))

使用 optim 进行优化。

optim(init_w, fn)
## $par
## [1]  0.1204827 -1.2438883  1.1023338  1.0212406
## 
## $value
## [1] 7.807847e-08
## 
## $counts
## function gradient 
##      111       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL

这个 par 元素包含你的权重。

- Richie Cotton

---

+1 很好的答案。我猜一个人也可以主观地加权考虑errs中每个值在最小化中的重要性？例如：errs <- c(sum(w) - 1, sum(x * w) - n) * c(5,1) - Marc in the box

哇，太棒了。终于我明白了优化器是如何正常工作的。非常感谢。非常好的答案。 - arodrisa

@Marcinthebox 或许在下一行将错误权重放置更语义化，这样你就可以将它们组合成一个单一的分数。我已经更新了答案来展示如何做到这一点。 - Richie Cotton

Answer 2

这个问题没有唯一的解。如果您尝试其他初始值来代替w，那么您很可能会从optim得到不同的结果。

这个问题可以被描述为解决一个欠定的线性方程组。

A <- matrix(c(rep(1,4),x), nrow=2,byrow=TRUE)
b <- matrix(c(1,n), nrow=2)

我们需要找到一个满足 A %*% w = b 的解决方案，但是应该选择哪个呢？最小范数解？还是其他的解？有无限多个解。可以使用矩阵A的伪逆来给出解决方案。请使用MASS包进行操作。

library(MASS)
Ag <- ginv(A)

最小范数解是：

wmnorm <- Ag %*% b

请查看线性方程组维基百科页面中的矩阵解法, 并使用A %*% wmnorm - b和fn(wmnorm)进行检查。

解为:

Az <- diag(nrow=nrow(Ag)) - Ag %*% A
w <- wmnorm + Az %*% z

其中z是一个任意向量，其包含了Az矩阵的ncol(Az)个元素。 现在生成一些解并进行检查。

xb <- wmnorm
z <- runif(4)
wsol.2 <- xb + Az %*% z
wsol.2
A %*% wsol.2 - b
fn(wsol.2)

z <- runif(4)
wsol.3 <- xb + Az %*% z
wsol.3
A %*% wsol.2 - b
fn(wsol.3)

当作为参数传递给 fn 时，您会发现这两个解决方案是有效的解决方案。它们与由 optim 找到的解决方案非常不同。您可以通过选择不同的起始点 init_w 来测试这一点，例如通过 init_w1 <- runif(4)/4。

- Bhas