SAS minimization problem - Proc Optmodel / proc IML

[Giraffe123]

Hi everyone,

I have the following problem which I need to solve with either PROC OPTMODEL or PROC IML. Unfortunately I'm not advanced in either.

I have a transition matrix = {0.1,0.4,0,5,

                                             0.6,0,4,0,

                                             0 , 0, 1};

I have periods = 12.

And then my problem statement is: Matrix_nperiod = transition_matrix**n_period,

and I have the error: abs(Matrix_nperiod - transition_matrix) which I would like to minimize. I have tried the following code to no avail: ( I haven't tried proc iml yet)

proc optmodel;

num transition_matrix {1..3,1..3} = [0.1,0.4,0,5,

                                                        0.6,0,4,0,

                                                            0 , 0, 1];

num n_period = 12;

impvar Matrix_nperiod = transition_matrix**n_period;

min error =Matrix_nperiod - transition_matrix;

solve;

print error;

quit;

Can somebody please help?

Kind Regards,

[Ksharp]

I think the best choice is PROC IML.

Calling @Rick_SAS  and @RobPratt 

[Rick_SAS]

I don't understand your problem. You say you want to minimize abs(Matrix_nperiod - transition_matrix), but what is the parameter in this problem. That is, what parameter are you allowed to change in order to find the optimal value?

Also, can you clarify the loss function? Do you mean the SUM of the difference between elements in the two matrices?

In your program, you set n_period = 12, and then compute

Matrix_nperiod = transition_matrix**n_period;

But Matrix_nperiod is a (known) constant, so the difference is a constant.

Maybe you want to vary the number of periods? If so, the answer is n_period=1, because then Matrix_nperiod = transition_matrix, and the difference is exactly 0.

[Giraffe123]

Sorry let me make a correction.

So actually you have:

Matrix_nperiod = One_period_matrix **n_period

and you error = abs(Matrix_nperiod - transition_matrix).

So the objective function is actually:

sum (error) in other words:

sum ((One_period_matrix **n_period) - transition_matrix)

In other words, One_period_matrix is your parameter

Sorry for the confusion 🙂

[Rick_SAS]

If I understand you correctly, you want to find a 3x3 matrix, A, such that A^12 is as close as possible to a given matrix, B? In other words, you seek the 12th-root of B?

What properties must A obey? Is it also a Markov transition matrix?

Are you sure you want the matrix that is closest according to sum(abs(diff))? Another choice is to use the L2 norm (Fobenius norm for matrices), which is sum(diff_{ij}^2).

[Rick_SAS]

Regarding your access to PROC IML, are you running SAS 9.4 or SAS Viya?  I ask because the SAS Viya version of SAS IML has some additional features that we could use to solve this problem.

[Giraffe123]

I'm using 9.4. And yes, it's a markov transition matrix.

[Rick_SAS · Accepted solution]

Here's a PROC IML solution. I replaced your loss function with a quadratic loss function because some nonlinear optimization methods perform better when the objective function is differentiable. 

proc iml;
start ObjFun(parms) global(periods, target);
   A = shape(parms, 3, 3);
   B = A**periods;
   error = sum(ssq(B - target));
   return error;
finish;

/* define the target matrix. We want to find transition matrix A
   such that A**12 is as close as possible to the target matrix. */
Target = {0.1 0.4 0.5,
          0.6 0.4   0,
          0     0   1};
periods = 12;

/* define linear and boundary constraints */
con = {  0 0 0  0 0 0  0 0 0 . .,  /* lower bound is 0 for all parameters */
         1 1 1  1 1 1  1 1 1 . .,  /* upper bound is 1 for all parameters */
         1 1 1  . . .  . . . 0 1,  /* sum of 1st row = 1 */
         . . .  1 1 1  . . . 0 1,  /* sum of 2nd row = 1 */
         . . .  . . .  1 1 1 0 1}; /* sum of 3rd row = 1 */

x = I(3);  /* start with the identity matrix */
x0 = rowvec(x);
opt = {0 2};
call nlpnra(rc, xres, "ObjFun", x0, opt) BLC=con;

A = shape(xres, 3, 3);
B = A**periods;
SSE = sum(ssq(B - target));
print A[L="Best Estimate"], B[L="Estimate**12"], SSE;

I should mention that this problem might not have a unique solution, and it could have local optima. In other words, if you use a different initial guess, you might obtain a different value for the local solution.

[RobPratt]

Here's a way to minimize the quadratic loss with PROC OPTMODEL:

proc optmodel;
   set ROWS = 1..3;
   set COLS = ROWS;
   num target {ROWS, COLS} = 
      [0.1 0.4 0.5,
       0.6 0.4   0,
       0     0   1];
   num numPeriods = 12;
   set PERIODS = 1..numPeriods;

   var X {i in ROWS, j in COLS} >= 0 <= 1;
   impvar XPower {i in ROWS, j in COLS, p in PERIODS} = 
      (if p = 1 then X[i,j] else sum {k in COLS} XPower[i,k,p-1] * XPower[k,j,p-1]);
   min SSE = sum {i in ROWS, j in COLS} (XPower[i,j,numPeriods] - target[i,j])^2;
   con RowSum {i in ROWS}:
      sum {j in COLS} X[i,j] = 1;

   solve with nlp / ms algorithm=as;
   print X;
   print {i in ROWS, j in COLS} XPower[i,j,numPeriods];
   print target;
quit;

As @Rick_SAS suggested might be the case, the problem has several local minima.

[Rick_SAS]

I wrote an article about this problem that has additional details:

"Estimate the pth root of a Markov transition matrix"