Function for Generating bistochastic matrix ?

Bistochastic matrix P is matrix where the sum each column or each row is 1, and aslo for n>=0 P^n is also bistochastic,

Another definition : a stochastic matrix P is bistochastic if P' is also stochastic . the question is :

Is there any predefined or sophisticated function to generate such matrices ?

I already tried a method using "magic" function :

 >>H=magic(10); % say we want a Bistoch of dimension n
 >>N=sum(H(1,:)); %  to get the Unique SUM
 >>P=H/N;

Cordially

[Matt J]

Is there any predefined or sophisticated function to generate such functions ?

Generate them from what?

[Youssef Khmou]

like generating random numbers with "rand" function, given the dimensions

[Cedric]

The "magic" solution seems to be quite efficient already; did you ask because you need to be able to generate random bistochastic matrices?

[Youssef Khmou]

hi Cedric, the truth is that solution respects the definition of Bistochastic matrix but the repartition of the data is quite regular, a well built function can generate such quantity with randomness quality and why not having some parameters related to well known PDFs , to see the regularity try:

   H=magic(100); % say we want a Bistoch of dimension n
   N=sum(H(1,:)); %  to get the Unique SUM
   P=H/N;
   imagesc(P);

[Cedric]

Hi Youssef, I asked precisely because of the regularity. I have no clean solution, but if I had to find some solution quickly, I would certainly go for yours, using two successive RANDPERM to permute rows and columns. It would not be optimal, but ok for a temporary approach I guess.

[Youssef Khmou]

alright, using Random permutation is good point , thanks

Answer 1

[Matt J]

Here's one idea. It uses interpMatrix ( Available Here ) to create circulant matrices, but any other method of making them would do.

    function P=bistoch(N)
    %Randomly generates an NxN bistochastic matrix, P

         for ii=1:2

             x=rand(N,1);
             x=x/sum(x);
             P=full(interpMatrix(x,1,length(x),1,'circ'));
             A{ii}=P(randperm(N),randperm(N));

         end

         w=rand;
         P=w*A{1}+A{2}*(1-w);

[Youssef Khmou]

hi Mtt J,

Thanks for the answer, your function has a quality of generating randomness,

Just a small detail : The function interpMatrix.m seems to have an error :

   ??? Error: File: interpMatrix.m Line: 136 Column: 15
   Expression or statement is incorrect--possibly unbalanced (, {, or [.

So in the line 136 :

          case 'max'

            [~,origin]=max(kernel);

i made an alteration :

               case 'max'

            [origin]=max(kernel);

And its functionning well,

Thanks .

[Matt J]

Hi Youssef,

No, that is not the fix you want. You're obviously using a (very!) old MATLAB version which doesn't support tilde arguments. Substitute this instead,

 [discard,origin]=max(kernel);

[Youssef Khmou]

alright, that is good advice, thanks :