How to get optim working with matrix multiplication inside the function to be maximized in R

Question

I am attempting to maximize a likelihood for a Matrix parameter of dimension 2x2. The Likelihood function needs to pass in a couple of fixed matrix parameters that the likelihood is also a function of. The data, denoted Y, and a covariance matrix, Sigma.star (which I am passing through as a lower triangular matrix), are necessary for the calculation but I would like to keep those fixed and run an optim function over this, in my code trying to optimize A

My issue is that optim seems to be erring from the fact it's optimizing something inside of an object I'm using for matrix algebra. Is there some way to make it work without programming every little calculation out?

The specific error is:

Error in diag(1, nrow = (m^2)) - A %x% A : non-conformable arrays

But A kronecker A should be an m^2 x m^2 matrix just like the identity…

Code:

library(MCMCpack)
library(mvtnorm)
set.seed(1000)


Likelihood.orig<-function(A, Y, Sigma.star){
Sigma<-xpnd(Sigma.star)
n<-nrow(Y)
if(is.vector(A)==TRUE){
  A<-as.matrix(A, nrow=nrow(Sigma), ncol=ncol(Sigma))
}
m<-nrow(A)
V<-matrix(solve(diag(1, nrow=(m^2))-A%x%A)%*%as.vector(Sigma), nrow=m, ncol=m)
temp1<- (-.5)*log(abs(det(V)))
temp2<- (-(n-1)/2)*log(abs(det(Sigma)))
temp3<- t(Y[,1, drop=FALSE]) %*% (solve(V)) %*% Y[,1, drop=FALSE]
terms<- numeric(n-1)
for(i in 2:n){
    terms[i-1]<- t(Y[,i, drop=FALSE] - A %*%Y[,i-1, drop=FALSE]) %*% (solve(Sigma)) %*% (Y[,i] - A %*%Y[,i-1])
}
return(temp1+temp2-.5*(temp3+sum(terms)))   
}




Generate.Y<-function(n, A, Sigma){
m<-nrow(A)
Y<-matrix(0, nrow=m, ncol=n)
V<-matrix(solve(diag(1, nrow=m^2)-A%x%A)%*%as.vector(Sigma), nrow=m, ncol=m)
Y[,1]<-rmvnorm(1, numeric(nrow(A)), V)
for(i in 2:n){
    Y[,i]<-A%*%Y[,i-1, drop=FALSE]+t(rmvnorm(1, mean = numeric(m), sigma = Sigma))
    }
return(Y)
}


n<-500
A.true<-matrix(c(.8, .3, 0, .5), nrow=2, ncol=2)
Sigma<-matrix(c(1, 0, 0, .5), nrow=2, ncol=2)
Y<-matrix(0, nrow=2, ncol=n)
Y<-Generate.Y(n, A.true, Sigma)
m=nrow(Y)
lower.Sigma<-vech(Sigma)



optim(par=c(1, 0, 0, 1), fn=Likelihood.orig, method="Nelder-Mead", 
  control=list(maxit=500, fnscale=-1), Sigma.star=lower.Sigma, Y=Y)

Answer 1

Your approach is correct, i.e., make optim optimize over a vector, and only turn that vector into a matrix inside the function you are trying to maximize.

However, you need to use matrix and not as.matrix to create that matrix. See the difference between:

as.matrix(1:4, nrow=2, ncol=2)  # wrong tool
#      [,1]
# [1,]    1
# [2,]    2
# [3,]    3
# [4,]    4

and

matrix(1:4, nrow=2, ncol=2)
#      [,1] [,2]
# [1,]    1    3
# [2,]    2    4

For problems of this type, I would highly recommend you learn the R debugging tools (browser, debug, debugonce, etc.). See General suggestions for debugging in R for examples.