RAS Algorithm Refinement

In the post http://yetanothermathprogrammingconsultant.blogspot.com/2013/06/simple-data-set-to-show-equivalence-of.html a simple RAS algorithm was shown. A slightly better way is to use a convergence criterion.

I got some feedback that users were not really familiar with how to write this in GAMS. So I thought this would be a good thing to share.

$ontext    Example from    Using PROC IML to do Matrix Balancing    Carol Alderman, University of Kansas    Institute for Public Policy and Business Research    MidWest SAS Users Group MWSUG 1992    In this version we stop on convergence $offtext sets    p 'products' /pA*pI/    s 'salesmen' /s1*s10/ ; table A0(*,*) 'estimated matrix, known totals'             s1   s2   s3   s4   s5   s6   s7   s8   s9  s10   rowTotal       pA   230  375  375  100    0  685  215    0   50    0       2029       pB   330  405  419  175   90  504  515    0  240  105       2798       pC   268  225  242    0   30  790  301   44  100    0       1998       pD   595  380  638  275   30  685  605   88  100  160       3566       pE   340  360  440  200   30  755  475   44  150    0       2794       pF   132  190  200    0    0  432  130    0    0    0       1071       pG   309  330  350  125    0  612  474    0   50   50       2305       pH   365  400  330  150   50  575  600   44  150  110       2747       pI   210  250  308  125    0  720  256    0  100   50       2015 colTotal  2772 2910 3300 1150  240 5760 3526  220  950  495 ; A0(p,'sum') = sum(s,A0(p,s)); A0('sum',s) = sum(p,A0(p,s)); A0(p,'diff') = A0(p,'sum') - A0(p,'rowTotal'); A0('diff',s) = A0('sum',s) - A0('colTotal',s); option A0:0; display A0; alias (p,i); alias (s,j); parameters    A(*,*)  'new RASsed matrix'    u(i)    'row total'    v(j)    'column total'    rho(i)   'used inside RAS algorithm'    sigma(j) 'used inside RAS algorithm' ; A(i,j) = A0(i,j); u(i) = A0(i,'rowTotal'); v(j) = A0('colTotal',j); abort$(smin((i,j),A(i,j))<0) "Entries A(i,j) should be nonnegative",A; abort$(smin(i,u(i))<0) "Entries u(i) should be nonnegative",u; abort$(smin(j,v(j))<0) "Entries v(j) should be nonnegative",v; set iter 'max iterations of RAS algorithm' /iter1*iter100/; scalars    continue /1/    tol /1e-7/    diff    count ; loop(iter$continue, * * step 1 row scaling *    rho(i) = u(i)/sum(j,A(i,j));    A(i,j) = rho(i)*A(i,j); * * step 2 column scaling *    sigma(j) = v(j)/sum(i,A(i,j));    A(i,j) =  A(i,j)*sigma(j); * * check for convergence *    diff = max(smax(i,abs(rho(i)-1)),smax(j,abs(sigma(j)-1)));    continue$(diff<tol) = 0;    count = ord(iter); ); abort$continue "No convergence"; option count:0; display "Iterations needed for convergence:",count; A(i,'rowTotal') = u(i); A('colTotal',j) = v(j); A(p,'sum') = sum(s,A(p,s)); A('sum',s) = sum(p,A(p,s)); A(p,'diff') = A(p,'sum') - A(p,'rowTotal'); A('diff',s) = A('sum',s) - A('colTotal',s); display "-------------------------------------------------------------------",         "Results",         "-------------------------------------------------------------------",         A;

I often use a large set to loop over: “loop(iter,”. This will make sure we always terminate (no infinite loop). The $ condition on the loop can be used to terminate the loop when convergence is observed: “loop(iter$continue,”. Finally we easily can check whether the loop terminated with or without convergence and give an appropriate message.