# R - Constrained optimization for a function taking a matrix

I would like to do constrained optimization for a function which takes a matrix as input. The matrix is sparse, representing a weighted adjacency matrix , and only the weights shall be subject to optimization.

An example matrix:

> AM
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
[1,] 0.00 0.00 0.00 0.00 0.05 0.00 0.14    0 0.00
[2,] 0.00 0.00 0.00 0.58 0.15 0.85 0.82    0 0.00
[3,] 0.00 0.00 0.00 0.00 0.71 0.92 0.60    0 0.00
[4,] 0.00 0.00 0.99 0.00 0.00 0.53 0.00    0 0.00
[5,] 0.00 0.00 0.00 0.00 0.00 0.89 0.43    0 0.00
[6,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00    0 0.98
[7,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00    0 0.02
[8,] 0.25 0.58 0.00 0.00 0.00 0.00 0.00    0 0.00
[9,] 0.00 0.00 0.00 0.00 0.00 0.00 0.00    0 0.00


My function

H <- function(AM) {
n <- dim(AM)[1]
Pin<-colSums(AM)[1:(n-2)]
Pout<-rowSums(AM)[1:(n-2)]
H<-0
for (i in 1:(n-2)) {
for (j in 1:(n-2)) {
if (AM[i,j]>0) {
H <- H - AM[i,j]/sum(AM)*log(AM[i,j]/sum(AM))
}
}
}
return(H)
}


The output for the given matrix:

> H(AM)
[1] 2.107764


What I would like to do is to maximize H(AM), given the position of non-zero entries in AM (zero-entries in AM shall not be changed) and a few other constraints, namely:

rowSums(AM)[n-1] == c
colSums(AM)[1:(n-2)] == rowSums(AM)[1:(n-2)]


The optimization process should ideally work for any kind of matrix input (so that I don't have to specify the position or even number of non-zero entries manually).

I hope this problem is specified well enough, otherwise please let me know! Any help is very much appreciated!!

Many thanks!

PS: I also asked this question on stack overflow (https://stackoverflow.com/questions/28394327/r-constrained-optimization-for-a-function-taking-a-matrix) but was referred to this platform.

• Is there any constraints on AM, e.g. should AM* be PSD? Feb 8, 2015 at 14:38
• AM should have only positive entries and is limited to describing acyclic graphs. Feb 8, 2015 at 15:42

Here's an outline of a solution. The trick is to make a function that that takes a vector with length equal to the number of non-zero entries in your matrix (AM) and returns a matrix with the appropriate format.

Here's an example of you you can create functions that will do just that given an initial value of your matrix AM:

# a factor function which returns functions that build
# matrices with non-zeros in the same entries as the
# argument AM
AM_builder_Factory <- function(AM,c){
force(c)
DIM <- dim(AM)
nonZeros <- which(as.vector(AM)>0)

# this inner function will be the return value of
# AM_builder_Factory
function(x){
out<-numeric(prod(DIM))
# Here I'm using the  exponential transformation so that the
# non-zero entries in the output matrix are all positive.
# Many other transformations are possible, and may be more
# efficient depending on the function H() being optimized.
out[nonZeros] <- exp(x)

# Here is where you'll have to implement these constraints
# rowSums(AM)[n-1] == c
# colSums(AM)[1:(n-2)] == rowSums(AM)[1:(n-2)]

dim(out) <- DIM
out
}
}

# the function 'AM_builder' takes a vector with
# numeric values and returns a matrix with the shape
# of the initial AM, and with zero and non-zero values
# in the same locations as AM
AM_builder <- AM_builder_Factory(AM)

# for example,
AM = matrix(c(1,0,0,
2,3,0,
0,4,5),nrow=3,byrow=T)
#>      [,1] [,2] [,3]
#> [1,]    1    0    0
#> [2,]    2    3    0
#> [3,]    0    4    5

builder <- AM_builder_Factory(AM)
builder(log(5:1))
#>      [,1] [,2] [,3]
#> [1,]    5    0    0
#> [2,]    4    3    0
#> [3,]    0    2    1


Now that you have a function builder that builds your matrices, you can use one of the optimization functions like nlm to optimize your function H()

# we need a set of starting parameters, which we'll take as the log
# of the non-zero values in AM.  [We're using log() of the non-zero
# values in AM because log() is he inverse of exp() which is used in
# the builder function to convert numeric values to positive values.]
start<- as.vector(AM)
start <- start[start>0]
val = nlm(f=function(x) - H(builder(x)), #  the function to be minimized
p=log(start)) # starting parameter values for the minimization

# the minimum value of H conditional on the format of AM
- val$minimum # the optimized matrix builder(val$estimate)


Note that the inner function inner function of AM_builder_Factory is not trivial b/c there are multiple ways to implement your constraints, and the constraints will necessarily reduce the number of independent (non-zero) entries in AM, so you'll have to adjust the definition the inner function to take a an argument (x) with length less than the number of non-zero entries in AM.

• Thank you so much Jthorpe, this was very very helpful! Actually I had been quite close already, but your answer helped me to finish up the task. Thank you so much for your quick, comprehensive and uttermost helpful post! Feb 9, 2015 at 15:51
• Very helpful indeed! I like the factory function approach! Mar 13, 2020 at 19:46