Thread Subject:

Updating a matrix (RAS algoritham)

Dear all,
I wish to know if it is possible to update a matrix using the so
called RAS algorithm.

What we have are, (1) the matrix (A0)to be updated, and (2) the row
and column sums "r"and "s" repsecively of the matrix to be derived
 from the original matrix.The procedure is the following.

a.Row scaling (k=1)

A1=diag(r./r0)*A0;

%where r0 is the row sum of the original matrix (A0)

b.Column scaling (k=2)

A2=A1*diag(s./s1);,Dear all,
I wish to know if it is possible to update a matrix using the so
called RAS algorithm.

What we have are, (1) the matrix (A0) to be updated, and (2) the row
and column sums "r" and "s" respectively of the matrix to be derived
 from the original matrix. The procedure is the following.

a. Row scaling (k=1)

A1=diag(r./r0)*A0;

%where r0 is the row sum of the original matrix (A0)

b.Column scaling (k=2)

A2=A1*diag(s./s1);

%where s1 is the column sum of the matrix (A1) derived in the previous
step.

 This process goes on until the row and column sums of the new matrix
(A*) equals that of the row and column sums "r"and "s" respectively.

My intuition is that it is possible to use some loop, but have not
been able to do it.


Any idea would be sincerely appreciated.


Kind Regards,

Jojo Jacob

%where s1 is the column sum of the matrix (A1) derived in the previous
step.

 This process goen on until the row and column sums of the new matrix
(A*) equals that of the row and column sums "r"and "s" repsecively.

My intution is that it is possible to use some loop, but have not been
able to do it.


Any idea would be sincerely appreciated.


Kind Regards,

Jojo Jacob

In article <xduie5od72z5@legacy>, j.jacob@tm.tue.nl (Jojo Jacob) wrote:

> I wish to know if it is possible to update a matrix using the so
> called RAS algorithm.
>
> What we have are, (1) the matrix (A0)to be updated, and (2) the row
> and column sums "r"and "s" repsecively of the matrix to be derived
> from the original matrix.The procedure is the following.
>
> a.Row scaling (k=1)
>
> A1=diag(r./r0)*A0;

(snip)

> This process goes on until the row and column sums of the new matrix
> (A*) equals that of the row and column sums "r"and "s" respectively.

This is a simple loop. You say that you were not
able to make it work. There are several reasons
why it can fail, besides programming error.

1. If the sums of the row sums and column sums are
not the same, it will never converge to your goals.

2. If A0 has many zero elements, so there are fewer
non-zeros than the total of its rows and columns,
it may be unsolvable. For example a diagonal matrix
A0 will cause this algorithm to fail, unless the
aim row sums and column sums are equal to each
other.

3. If A0 has any all zero rows, or all zero columns
it will fail.

There are other, less easy to describe, circumstances
under which your problem is unsolvable if we are
required to keep the sparsity structure of A fixed.

Below is a simple code which iterates your algorithm.

We might consider other variations of your problem.
For example, find a matrix B, (the same size and shape
as A0) such that A = A0 + B, A has the desired row
and column sums, and the perturbation matrix B is
minimized under some norm. This variation is solvable
using linear algebra, and as such, does not require
iterations to converge. It also always has a solution
as long as the aim row and column sums are consistent.

HTH,
John D'Errico


% ==========================================================
function A = ras(A0,rowsum,colsum)
tol=1.e-12;
rowsum=rowsum(:);
colsum=colsum(:)';

if (abs(sum(rowsum)-sum(colsum))/abs(sum(rowsum)))>tol
  error 'rowsum and colsum must both sum to the same amount'
end
  
A=A0;
flag=1;
iter=0;
while flag
  % row sums
  rs=sum(A,2);
  A=diag(rowsum./rs)*A;
  
  % column sums
  cs=sum(A,1);
  A=A*diag(colsum./cs);
  
  rerr=max(abs(rs-rowsum))/sum(rowsum);
  cerr=max(abs(cs-colsum))/sum(rowsum);
  
  if (rerr<=tol) & (cerr<=tol)
    flag=0;
  end
  iter=iter+1;
  
  disp(['iter,rerr,cerr: ',num2str([iter,rerr,cerr])])
  
end

--

John D'Errico wrote:
>
>
> This is a simple loop. You say that you were not
> able to make it work. There are several reasons
> why it can fail, besides programming error.
>
> 1. If the sums of the row sums and column sums are
> not the same, it will never converge to your goals.
>
> 2. If A0 has many zero elements, so there are fewer
> non-zeros than the total of its rows and columns,
> it may be unsolvable. For example a diagonal matrix
> A0 will cause this algorithm to fail, unless the
> aim row sums and column sums are equal to each
> other.
>
> 3. If A0 has any all zero rows, or all zero columns
> it will fail.
>
> There are other, less easy to describe, circumstances
> under which your problem is unsolvable if we are
> required to keep the sparsity structure of A fixed.
>
> Below is a simple code which iterates your algorithm.
>
> We might consider other variations of your problem.
> For example, find a matrix B, (the same size and shape
> as A0) such that A = A0 + B, A has the desired row
> and column sums, and the perturbation matrix B is
> minimized under some norm. This variation is solvable
> using linear algebra, and as such, does not require
> iterations to converge. It also always has a solution
> as long as the aim row and column sums are consistent.
>
> HTH,
> John D'Errico
>
>
> % ==========================================================
> function A = ras(A0,rowsum,colsum)
> tol=1.e-12;
> rowsum=rowsum(:);
> colsum=colsum(:)';
>
> if (abs(sum(rowsum)-sum(colsum))/abs(sum(rowsum)))>tol
> error 'rowsum and colsum must both sum to the same amount'
> end
>
> A=A0;
> flag=1;
> iter=0;
> while flag
> % row sums
> rs=sum(A,2);
> A=diag(rowsum./rs)*A;
>
> % column sums
> cs=sum(A,1);
> A=A*diag(colsum./cs);
>
> rerr=max(abs(rs-rowsum))/sum(rowsum);
> cerr=max(abs(cs-colsum))/sum(rowsum);
>
> if (rerr<=tol) & (cerr<=tol)
> flag=0;
> end
> iter=iter+1;
>
> disp(['iter,rerr,cerr: ',num2str([iter,rerr,cerr])])
>
> end
>
> --


Hi John


Thanks very much for the suggestions and the code. I tried it for the
following example.


A0=[20,40,80;20,20,40;20,20,0];% the original matrix
rowsum=[150;70;50];%row sum of the desired matrix
colsum=[70;70;130];%column sum of the desired matrix.


However, the following message comes.


??? Strings passed to EVAL cannot contain function declarations.


Excuse me if the question sounds silly as I am quite new to matlab.


Thanks again, Jojo

Hi John,


It worked. Many thanks,


Jojo

In article <eeaf2f6.1@WebX.raydaftYaTP>, "Jojo Jacob"
<j.jacob@tm.tue.nl> wrote:

> Thanks very much for the suggestions and the code. I tried it for the
> following example.
>
>
> A0=[20,40,80;20,20,40;20,20,0];% the original matrix
> rowsum=[150;70;50];%row sum of the desired matrix
> colsum=[70;70;130];%column sum of the desired matrix.
>
>
> However, the following message comes.
>
>
> ??? Strings passed to EVAL cannot contain function declarations.

I can only make a wild guess about what it is
you have done to get that error. Have you learned
how to create m-files and where to save them?
It appears that you were trying to execute (using
eval) the code I posted (which is a function.)

Here is what happens in matlab for me:


»A0=[20,40,80;20,20,40;20,20,0];% the original matrix
»rowsum=[150;70;50];%row sum of the desired matrix
»colsum=[70;70;130];%column sum of the desired matrix.
»ras(A0,rowsum,colsum)
iter,rerr,cerr: 1 0.037037 0.056878
iter,rerr,cerr: 2 0.0078652 0.0044432
iter,rerr,cerr: 3 0.0016275 0.0009853
iter,rerr,cerr: 4 0.0003676 0.00022309
iter,rerr,cerr: 5 8.3357e-05 5.0603e-05
iter,rerr,cerr: 6 1.8913e-05 1.1482e-05
iter,rerr,cerr: 7 4.2915e-06 2.6054e-06
iter,rerr,cerr: 8 9.7383e-07 5.9122e-07
iter,rerr,cerr: 9 2.2098e-07 1.3416e-07
iter,rerr,cerr: 10 5.0146e-08 3.0444e-08
iter,rerr,cerr: 11 1.13792e-08 6.9084e-09
iter,rerr,cerr: 12 2.5822e-09 1.56767e-09
iter,rerr,cerr: 13 5.85958e-10 3.55738e-10
iter,rerr,cerr: 14 1.32967e-10 8.07248e-11
iter,rerr,cerr: 15 3.01731e-11 1.83182e-11
iter,rerr,cerr: 16 6.84692e-12 4.15673e-12
iter,rerr,cerr: 17 1.55372e-12 9.43285e-13
iter,rerr,cerr: 18 3.52561e-13 2.14005e-13
ans =
       22.923 34.537 92.54
       18.559 13.981 37.46
       28.518 21.482 0


HTH,
John

--

In article <eeaf2f6.2@WebX.raydaftYaTP>, "Jojo Jacob"
<j.jacob@tm.tue.nl> wrote:

> Hi John,
>
>
> It worked. Many thanks,
>
>
> Jojo

I'm glad. Now you might want to add a few things
to that code.

A well written numerical code would have a limit
on the maximum iterations, to prevent an infinite
loop in the event of non-convergence. Perhaps
return an error, or at least a warning message if
the code fails in a reasonable number of iterations.

You might allow the user to specify the tolerance,
with a default value if they do not specify it.
You also might want to turn off the display after
every iteration. This will get boring after much
use.

Finally, I would look carefully at the tolerance
I chose, and the relative error computation I used
on the row and column sums. Determine whether they
are reasonable for both large and small matrices,
at least as large as you will expect to see.

Finally, add some documentation to the beginning
so you know what the arguments mean when you look
back at this code in a year.

John

--

Hi John,


Your comments are noted. I shall be trying to implement your RAS code
for different initial input-output matrix formats (eg. for national
IO tables, with and with out imports or value added etc.). I shall
keep you posted on the developments.


Thanks again.


Jojo