Why this thing happens with random matrix such tha

Question

I would like to generate a random matrix whose rows sum up to one. I found this question and answers How to create a random matrix such that all rows sum to 1 which do what I exactly want.

The problem is when I did it, sometimes the sum is not exactly one.

mat      = rand(2, 2);
rowsum   = sum(mat, 2);
mat      = bsxfun(@rdivide, mat, rowsum);

Sometimes I get something like this:

sum_test = sum(mat, 2)-1
1.0e-15 *
      0
-0.2220;

I do not know why?

Answer 1

Matlab uses double-precision numbers, which means there are only a finite set of numbers that can be represented. The gap between double-precision numbers can be found in Matlab using the eps function. For a value of 1, which is pretty much the value you get when you sum up each row, eps(1) = 2.220446049250313e-16, which is the same magnitude numerical error you're getting. The error is introduced when you divide each value in the matrix. More information about Matlab precision here.

Even doing something like subtracting the error from the first value in each row doesn't solve your problem completely (although it does help). For example,

mat = rand(10,10);
rowsum = sum(mat,2);
mat = bsxfun(@rdivide,mat,rowsum);
sum_test = sum(mat,2)-1    %finite error
mat(:,1) = mat(:,1) - sum_test; %subtract the error from the first column in each row
sum_test2 = sum(mat,2)-1   %still get error

results in

sum_test = [
-1.110223024625157e-16
                     0
                     0
                     0
                     0
-2.220446049250313e-16
                     0
-2.220446049250313e-16
-1.110223024625157e-16
-2.220446049250313e-16]

and

sum_test2 = [
 2.220446049250313e-16
                     0
                     0
                     0
                     0
                     0
                     0
                     0
-1.110223024625157e-16
                     0]

(You can show that even (sum(mat,2)-0.5)-0.5 gives a different answer than sum(mat,2)-1)