CODEWITHSUNDEEP
rloopsfor-looplogic
forecaster
forecaster

Reputation: 1159

matrix programming sum function

I have a square matrix "a" below as an example see below. matrix a, is nxn square matrix.

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)

I'm trying to write a function(x) in R as below in order to provide this to an optimization routine. I'm trying to minimize function(x), where x is unknown. x is vector.

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
}

Can you please help in programming this logic and function in R?

Much appreciated

Upvotes: 1

Views: 196

Answers (1)

Ferdinand.kraft
Ferdinand.kraft

Reputation: 12819

You can use this:

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),]
}

EDIT: I've phrased the conditions as i!=j & j!=k & i!=k & j>i to be more readable and include the condition you mentioned in comments.

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

Example:

sumx(1:5)
#[1] 3584.542

Note that I've replaced a[i,j] with a[cbind(i,j)] to allow vectorized access to elements of the matrix.

You can now put sumx to optimize, but it may be better to save comb3(length(x)) and the part that does not depend on x as a global object to reduce computation time, like this:

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))
}

For minimization, you can use optim. Note that I've found two different attractors:

> 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

Upvotes: 1

Related Questions