turn on suggestions

Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

Find a Community

Topic Options

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

03-19-2014 01:36 PM

I am converting SAS OPTMODEL codes to SAS/IML, though I am not familiar with both. OPTMODEL is easy to understand but I know IML is much more powerful.

If in OPTMODEL, the variable is defined a matrix, it's easy to write constraint, such as x[i,j]-x[i,j+1]>=0.01.

How can I write the constraint in IML? The online example shows how to write constraint if the variable is a vector, I am not sure how to write when we are dealing with a matrix.

for example, how to convert x[i,j]-x[i,j-1]>=0.01?

Thanks

Accepted Solutions

Solution

03-19-2014
02:33 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-19-2014 02:33 PM

What an interesting question! I haven't done this before, but I think you have two options.

1) Put the contraint as part of the objective function. You can add a penalty to the objective function to penalize solutions that do not satisfy the constraint. For example, the following function returns 0 when x>a and is quadratic when x<a. You can then use the sum-of-squares of the penalized matrix as a penalty function:

proc iml;

start Penalty(m, a);

v = m-a;

return( choose(v>0, 0, v##2) );

finish;

x = j(4,6,1); call randseed(1); call randgen(x,"normal");

/* add lines like these to objective function to penalize elements when x[i,j]-x[i,j+1]>=0.01 */

d = x[,1:ncol(x)-1] - x[,2:ncol(x)];

p = ssq(penalty(d, 0.01));

print p;

2) The BLC matrix (which stands for "boundary and linear constraints") can be used to specify constraints on adjacent elements of x. Look at the segment of the IML documentation

"Nonlinear Optimization Examples"--> "Details"--> "Parameter Constraints"

You would treat the elements of the matrix as if they were a vector (in row-major order).

Personally, I like the constraint approach.

All Replies

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-19-2014 02:02 PM

Actually, I would re-consider changing from OPTMODEL to IML. OPTMODEL is specifically designed for optimization and offers the newest and fastest algorithms. You can combine the best features of both IML and OPTMODEL in PROC IML, by using a submit block to run OPTMODEL specifically for the optimization phase, and use IML for any other calculations more easily programmed in a matrix-based mathematical programming language.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Hutch_sas

03-19-2014 02:19 PM

Thank you so much.

Solution

03-19-2014
02:33 PM

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-19-2014 02:33 PM

What an interesting question! I haven't done this before, but I think you have two options.

1) Put the contraint as part of the objective function. You can add a penalty to the objective function to penalize solutions that do not satisfy the constraint. For example, the following function returns 0 when x>a and is quadratic when x<a. You can then use the sum-of-squares of the penalized matrix as a penalty function:

proc iml;

start Penalty(m, a);

v = m-a;

return( choose(v>0, 0, v##2) );

finish;

x = j(4,6,1); call randseed(1); call randgen(x,"normal");

/* add lines like these to objective function to penalize elements when x[i,j]-x[i,j+1]>=0.01 */

d = x[,1:ncol(x)-1] - x[,2:ncol(x)];

p = ssq(penalty(d, 0.01));

print p;

2) The BLC matrix (which stands for "boundary and linear constraints") can be used to specify constraints on adjacent elements of x. Look at the segment of the IML documentation

"Nonlinear Optimization Examples"--> "Details"--> "Parameter Constraints"

You would treat the elements of the matrix as if they were a vector (in row-major order).

Personally, I like the constraint approach.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

03-19-2014 03:08 PM

Thank you very much.

I will take a look at the option 1.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

03-20-2014 11:30 AM

Hi Rick,

Thanks for your reply. I watched your video on youtube awhile ago and I am following your blog now.

I prefer your option 1 and already did some initial work. I think it's easy to use penalty function.

I have a followup question:

in OPTMODEL, the constraint is x[i,j]-x[i,j+1]>=0.01, however i is from 1 to nrow(x), but j is from** i,** instead of **1.**

How can I do it in one statement?

Thanks

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-20-2014 02:16 PM

What I would do is to put missing values in the lower triangular elements of x (or a copy of x) and then proceed as before. An easy way to access the lower triangular elements is to use the ROW function and COL function, which were introduced in SAS/IML 12.3. If you are using an earlier version, use the definitions in the following blog post, which also shows how to define and use these function: Defining banded and triangular matrices: Another approach - The DO Loop

The revised program might look something like this:

/* add lines like these to objective function to penalize elements when x[i,j]-x[i,j+1]>=0.01, j>=i */

y = x; /* copy x */

y[ loc(row(x)>col(x)) ] = .; /* set lower triangular elements to missing */

d = y[,1:ncol(x)-1] - y[,2:ncol(x)];

p = ssq(penalty(d, 0.01));

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

03-26-2014 01:26 PM

Thanks Rick.

Now I finished converting all constraints from OPTMODEL to IML.

The structure is like this:

PROC IML;

START PENALTY FUNCTION;

FINISH;

START OBJECTIVE FUNCTION;

NEW OBJECTIVE FUNCTION=ORIGINAL FUNCTION + PENALTY

FINISH;

now i tried to set up the optimization part.

The initial guess is p=j(4,4,0.1) /*4 by 4 matrix*/

The option is opt={1,4} /*min function and print detail*/

However I don't know how to set up constraint here. Or I don't know how to tell SAS this 4 by 4 matrix does not need any constraint (constraints are in objective function).

Any suggestion?

Thanks

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-26-2014 06:41 PM

If the parameters are unconstrined, omit the constraint matrix. It is an optional argument.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

03-27-2014 11:24 AM

I tested my codes as followings:

START PENALTY FUNCTION;

FINISH;

START OBJECTIVE_FUNCTION(AAA);

NEW OBJECTIVE FUNCTION=ORIGINAL FUNCTION + PENALTY

RETURN(ZZZ);

FINISH;

AAA={ ............};

test=OBJECTIVE_FUNCTION(AAA);

print test;

and the result is what I expect. So this part is correct.

now I remove AAA and test part, then write:

x0=j(4,4,0.1)

opt={1,4};

call nlpnra(rc,result,'fun',x0,opt);

I got the error saying "NLPNRA call: Error in argument X0".

Either I should not use NLPNRA function nor the syntax of X0 is wrong.

Which one?

Thanks

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-27-2014 11:39 AM

Does the following work?

test=OBJECTIVE_FUNCTION(x0);

How about this?

call nlpnra(rc,result,'fun',AAA,opt);

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

03-27-2014 12:03 PM

tried test=OBJECTIVE_FUNCTION(x0);

no error.

tried call nlpnra(rc,result,'fun',AAA,opt);

no luck, same error message.

however, if I change x0 to x0=j(1,4,0.1);

well, I got this

ERROR: (execution) Invalid subscript or subscript out of range.

once again, thanks for your help

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-27-2014 12:58 PM

if you change x0 to x0=j(1,4,0.1); then x0 is a vector with four elements, which is why you get "Invalid subscript."

There is nothing more I can offer. You have not posted any code or data, so it is unlikely that I can guess the source of your problem.

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to Rick_SAS

03-27-2014 01:09 PM

Hi Rick, the following is the sample codes I have. (I only keep one penalty value)

proc iml;

start penalty1(m,a);

v=m-a;

return(choose(v>0,0,v##2));

finish;

start sumsqr(a) global(m1);

d1=a[1:nrow(a)-1,1:ncol(a)];

p1=ssq(penalty1(d1,0.01));

ccc=(a-m1)##2;

zzz=sum(ccc)+p1;

return (zzz);

finish;

aaa={0.8 0.15 0.04 0.01,

0.1 0.8 0.09 0.01,

0.05 0.1 0.8 0.05,

0 0 0 1};

m1={0.7 0.25 0.04 0.01,

0.2 0.7 0.09 0.01,

0.05 0.2 0.7 0.05,

0 0 0 1};

test=sumsqr(aaa);

print test;

x0=j(4,4,0.1);

opt={1,4};

call nlpnra(rc, result, "sumsqr", x0, opt);

quit;

- Mark as New
- Bookmark
- Subscribe
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

Posted in reply to neilxu

03-27-2014 04:58 PM

You are trying to minimize a function, right? So use

opt={0,4};

Also, I just realized that the NLP functions are expecting column vectors, not matrices.

So do this:

start sumsqr(_a) global(m1);

a = shape(_a, nrow(m1), ncol(m1));

...

finish;

and then in the main program use

x0 = colvec(x0);

before you call NLPNRA.