Difficulties using Stata's optimize() for a minimization problem with constraints

[Lexy Thompson]

Hi,

I'm trying to calculate weights for constituents in an index using Stata's optimize() function. I have some code so far, but I'm having some difficulties getting it to run.

The equation I am trying to minimize is W=min_w || (X^T*Omega^-1*X)*W - X^T*Omega^-1 ||² subject to an upper bound and a lower bound for each constituent. This was previously done in Matlab using the command lsqlin.

X is a 134 x 4 matrix with rows for each constituent and columns for a geographic region (E, M, S, W). The data in X is the percent of each constituent's assets in each geographic region. Therefore, each row in matrix X adds to 1.

X^T is the transpose of X.

Omega^-1 is a 134 x 134 diagonal matrix with a constituent's enterprise value along the diagonal.

I also have variables in my data for the upper and lower bounds for each constituent.

W is what I am trying to minimize.

Ideally, I would like to have the results be in matrix form/some form that I could export back to Stata (maybe using the st_store option?). I've pasted my code below. When I try and run the last line, I get error r3200 (conformability). I don't know if I've chosen the correct option for any of the commands, particularly the evaluatortype() and params() lines. I was also having a difficulty understanding how to add the upper and lower bound constraints for each constituent.

CODE:
use indexdata

/* start optimization process */
mata

X=st_data(., "XRegE XRegM XRegS XRegW") /* to create X matrix */
Xt=X' /* transpose */

omegadiag=st_data(., "omegadiag")
omega=diag(omegadiag) /* create Omega^-1 */

C=Xt*omega
D=C*X

LB=st_data(., "LB") /* lower bound constraint for each constituent 134 x 1 matrix */
UB=st_data(., "UB") /* upper bound constraint for each constituent 134 x 1 matrix */

void weights(todo, p, W, g, H)
{
W=(D*W - C)^2
}

S = optimize_init()
optimize_init_which(S,"min")
optimize_init_evaluator(S, &weights())
optimize_init_evaluatortype(S, "gf2")
optimize_init_params(S, 0)
W=optimize(S)

[Christophe Kolodziejczyk]

Some quick comments


1. According to your evaluator you try to minimize with respect to p, which is passed as an argument, but is not used in the evaluator. I guess the objective function should be

W= (D*p-C):^2

Note the colon operator to take the square element-wise
2. You use the gf2 evaluator which requires to compute the gradient an the hessian analytically. Use the gf0 evalatuor in order to compute the gradient and the hessian numerically.
3. D and C are unkwown to the evaluator since they are not passed to the function. So add these in the list of arguments. Use optimize_init_argument() accordingly.
4. Finally optimize() is a tool for optimizing a function without inequality constraints. I don't know if it's possible to trick optimize, i.e. modify the objective function, to incorporate these inequality constraints.
(5.) Not strictly necessary, but you can use the cross() function to compute cross-products like X'*Y and X'*diag(w)*X

...

C=cross(X,omega) // Xt*omega 
D=cross(X,omegadiag,X) // C*X = X'diag(omega)*X

LB=st_data(., "LB") /* lower bound constraint for each constituent 134 x 1 matrix */
UB=st_data(., "UB") /* upper bound constraint for each constituent 134 x 1 matrix */

void weights(todo, p,C,D, W, g, H) // <= add C and D as arguments
{
W=(D*p- C):^2
}

/* to be consistent with your notation
void weights(todo, W,C,D, f, g, H) // <= add C and D as arguments
{
f=(D*W- C):^2
}

*/

S = optimize_init()
optimize_init_which(S,"min")
optimize_init_evaluator(S,&weights())
optimize_init_evaluatortype(S, "gf0") // change evaluator
optimize_init_argument(S,1,C)   // Tell the optimizer what the two new arguments 
optimize_init_argument(S,2,D)
optimize_init_params(S, 0)
W=optimize(S)

[Lexy Thompson]

Hi Christophe Kolodziejczyk . Thank you for the comments.

When I ran through your suggestions I still received the conformability error (r3200). I've included the complete code I was using with your suggestions below.

In order to incorporate the LB and UB constraints, do you think there is a better Stata function to use?

Code:
use indexdata
mata

X=st_data(., "XRegE XRegM XRegS XRegW") /* to create X matrix */
Xt=X' /* transpose */
omegadiag=st_data(., "omegadiag")
omega=diag(omegadiag) /* create Omega-1 */
C=Xt*omega
D=C*X

void weights(todo, p,C,D, W, g, H)
{
W=(D*p-C):^2
}

S = optimize_init()
optimize_init_which(S,"min")
optimize_init_evaluator(S,&weights())
optimize_init_evaluatortype(S, "gf0")
optimize_init_argument(S,1,C)
optimize_init_argument(S,2,D)
optimize_init_params(S, 0)
W=optimize(S)

[Lexy Thompson]

I posted an incorrect update. It doesn't seem that I can delete my comment, so I am posting this update.

[Christophe Kolodziejczyk]

p has to be transposed

void weights(todo, p,C,D, W, g, H)
{
W=(D*p'-C):^2
}

Stata has no method to specifically optimize functions with inequality constraints. You could use a transformation of your parameters to ensure that they are within the bounds (with the logit() function for example), but you won't allow some of the parameters to be on the boundaries. Generally this kind of problems requires different methods than the line-search algorithms available in Stata/Mata.

[Lexy Thompson]

Christophe Kolodziejczyk

I appreciate the help. I keep getting the same error (3200) though. I'm trying to have Mata return a matrix, which I then can export back to Stata. I originally was completing this in MATLAB using the lsqlin function (https://www.mathworks.com/help/optim/ug/lsqlin.html). Do you know if there is another command in Stata or Mata that I should be using?

[Christophe Kolodziejczyk]

There seems to be a general problem with the comformability of your matrices in your objective function.
C is a matrix 4 x 134 and D a matrix 4 x 4
For the expression in the objective function to be comformable W should be a matrix 4 x 134.
If you're trying to optimize your function with optimize(), you expect W in your expression to be vector. It is also the case with lsqlin. Are you sure the objective function is correct?


There is another problem. The initial values you feed the problem are a scalar (0), it should be a rowvector instead

optimize_init_params(S, J(1, # of parameters of the problem, 0)

use st_matrix() to send the Mata matrix in a Stata matrix

st_matrix("st_W",W)

puts the Mata matrix W into a Stata matrix st_W

[Lexy Thompson]

Christophe Kolodziejczyk
I had also tried another version where D was a 536 x 536 matrix and C was a 536 x 1 matrix. So W (or do you mean p?) would be a column vector and it still didn't work. Maybe because of the optimize_init_params statement? For this version, would the optimize_init_params statement be:
optimize_init_params(S, J(1, 536, 0) ?

I was also wondering if the function should be

Code:
void weights(todo, p, C, D, W, g, H)
{
W= D*p-C
}

to correspond with lsqlin.