[R] Re : Generate a random bistochastic matrix

[Martin Maechler]

Thank you, Ravi,

>>>>> "Ravi" == Ravi Varadhan <rvaradhan at jhmi.edu>
>>>>>     on Mon, 16 Oct 2006 18:54:16 -0400 writes:

    Ravi> Martin, I don't think that a doubly stochastic matrix
    Ravi> can be obtained from an arbitrary positive rectangular
    Ravi> matrix.  There is a theorem by Sinkhorn (Am Math Month
    Ravi> 1967) on the diagonal equivalence of matrices with
    Ravi> prescribed row and column sums.  It shows that given a
    Ravi> positive matrix A(m x n), there is a unique matrix DAE
    Ravi> (where D and E are m x m and n x n diagonal matrices)
    Ravi> with rows, k*r_i (i = 1, ..., m), and column sums, c_j
    Ravi> (j=1,...,n) such that k = \sum_j c_j / \sum_i r_i.
    Ravi> Therefore, the alternative row and column
    Ravi> normalization algorithm (same as the iterative
    Ravi> proportional fitting algorithm for contingency tables)
    Ravi> will alternate between row and column sums being
    Ravi> unity, while the other sum alternates between k and
    Ravi> 1/k.

    Ravi> Here is a slight modification of your algorithm for
    Ravi> the rectangular case:


    Ravi> bistochMat.rect <- function(m,n, tol = 1e-7, maxit = 1000) {
    Ravi> ## Purpose: Random bistochastic *square* matrix (M_{ij}):
    Ravi> ##	        M_{ij} >= 0;  sum_i M_{ij} = sum_j M_{ij} = 1   (for all i,
    Ravi> j)
    Ravi> ##
    Ravi> ----------------------------------------------------------------------
    Ravi> ## Arguments: n: (n * n) matrix dimension;
    Ravi> ##
    Ravi> ----------------------------------------------------------------------
    Ravi> ## Author: Martin Maechler, Date: 16 Oct 2006, 14:47
    Ravi> stopifnot(maxit >= 1, tol >= 0)
    Ravi> M <- matrix(runif(m*n), m,n)
    Ravi> for(i in 1:maxit) {
    Ravi> rM <- rowSum(M)
    Ravi> M <- M / rep((rM),n)
    Ravi> cM <- colSum(M)
    Ravi> M <- M / rep((cM),each=m)
    Ravi> if(all(abs(c(rM,cM) - 1) < tol))
    Ravi> break
    Ravi> }
    Ravi> ## cat("needed", i, "iterations\n")
    Ravi> ## or rather
    Ravi> attr(M, "iter") <- i
    Ravi> M
    Ravi> }

    Ravi> Using this algorithm we get for an 8 x 4 matrix, for example, we get:

    >> M <- bistochMat.rect(8,4)
    >> apply(M,1,sum)
    Ravi> [1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
    >> apply(M,2,sum)
    Ravi> [1] 1 1 1 1

"Of course" I had tried that too, before I posted and limited
the problem to square matrices.

    Ravi> Clearly the algorithm didn't converge according to
    Ravi> your convergence criterion, but the row sums oscillate
    Ravi> between 1 and 0.5, and the column sums oscillate
    Ravi> between 2 and 1, respectively.

indeed, and I had tried similar examples.

The interesting thing is really the theorem you mention
a consequence of which seems to be that indeed, simple row and
column scaling iterations would not converge.

Intuitively, I'd still expect that relatively simple
modification of the algorithm should lead to convergence.

Your following statement seems to indicate so,
or do I misunderstand its consequences?

    Ravi> It is interesting to note that the above algorithm
    Ravi> converges if we use the infinity norm, instead of the
    Ravi> 1-norm, to scale the rows and columns, i.e. we divide
    Ravi> rows and columns by their maxima.