1

我在下面有一个方阵“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中编程这个逻辑和功能吗？

非常感激

4

1 回答 1

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更具可读性，并包括您在评论中提到的条件。

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
于 2013-09-01T16:57:41.793 回答