Thursday, June 6, 2013

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.