矩阵规划求和函数

Question

我在下面有一个方阵“a”作为示例，请参见下文。矩阵a，是nxn方阵。

a = matrix( 
  c(1, 5 , 3, 7 , 3,
    5, 1, 2, 2, 4,
    3, 2 , 1, 2,4,
    7, 2, 2,1,3,
    2, 4,4 ,3 , 1   
   ),ncol = 5,nrow =5) 

我正在尝试在 R 中编写如下所示的函数（x），以便将其提供给优化例程。我试图最小化函数（x），其中 x 是未知的。 x 是向量。

sumx <- function(x) {

sum(((a[i,j]*a[j,k])-(x[i]/x[j]))^2) for all i,j,k such that i not eq to j not eq to k
}

您能帮忙在 R 中编程这个逻辑和函数吗？

非常感谢

Answer 1

你可以用这个：

comb3 <- function(n){
    result <- expand.grid(i=1:n,j=1:n,k=1:n)
    result[with(result, i!=j & j!=k & i!=k & j>i),]
}

编辑：我将条件表述为i!=j & j!=k & i!=k & j>i 以便更具可读性并包含您在 cmets 中提到的条件。

sumx <- function(x) {
    sum(with(comb3(length(x)), ((a[cbind(i,j)]*a[cbind(j,k)])-(x[i]/x[j]))^2))
}

例子：

sumx(1:5)
#[1] 3584.542

请注意，我已将 a[i,j] 替换为 a[cbind(i,j)] 以允许对矩阵元素进行矢量化访问。

现在可以把sumx放到优化中，不过最好还是把comb3(length(x))和不依赖x的部分保存为全局对象来减少计算时间，像这样：

y <- within(comb3(nrow(a)), b <- a[cbind(i,j)]*a[cbind(j,k)])

sumx <- function(x) {
    sum(with(y, (b-(x[i]/x[j]))^2))
}

为了最小化，您可以使用optim。请注意，我发现了两种不同的吸引子：

> optim(rep(1,5), sumx)
$par
[1] 1.9739966 1.5882750 1.5626338 0.1592725 0.1521839

$value
[1] 1436.526

$counts
function gradient 
     502       NA 

$convergence
[1] 1

$message
NULL

> optim(1:5, sumx)
$par
[1] 5.4254668 4.3857303 4.3029354 0.4374246 0.4199909

$value
[1] 1436.503

$counts
function gradient 
     218       NA 

$convergence
[1] 0

$message
NULL