[R] Re : Generate a random bistochastic matrix

[Ravi Varadhan]

Martin,

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

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

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

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

> M <- bistochMat.rect(8,4)
> apply(M,1,sum)
[1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5
> apply(M,2,sum)
[1] 1 1 1 1
>
Clearly the algorithm didn't converge according to your convergence
criterion, but the row sums oscillate between 1 and 0.5, and the column sums
oscillate between 2 and 1, respectively.

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


Best,
Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology 

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: rvaradhan at jhmi.edu

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: r-help-bounces at stat.math.ethz.ch
[mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Martin Maechler
Sent: Monday, October 16, 2006 10:30 AM
To: Florent Bresson
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] Re : Generate a random bistochastic matrix

A simpler approach --- as in similar problems ---
is to use an iterative solution  which works quite fast for any 'n'.
Interestingly, the number of necessary iterations seems to
*decrease* for increasing 'n' :

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

attr(M <- bistochMat(70), "iter")
## typically:
## [1] 8

attr(M <- bistochMat(10), "iter")
## -> 13, 16, 15, 17, ...

set.seed(101) ; attr(M <- bistochMat(10), "iter")             ## 15
set.seed(101) ; attr(M <- bistochMat(10, tol = 1e-15), "iter")## 32

------------------------

Initially, I tried to solve the general [ n x m ] case.
It seems that needs a slightly smarter "bias correction"
instead of just 'sum(M) / n'.
If someone figures that out, please post your solution!

Regards,
Martin Maechler, ETH Zurich

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.