[R] Re : Generate a random bistochastic matrix

[Ravi Varadhan]

Martin,

Sinkhorn's theorem excludes the possibility of obtaining a doubly stochastic
matrix of the form D%*%A%*%E, which is diagonally equivalent to a given
positive rectangular matrix A.  But it "doesn't" say that one can't obtain a
doubly stochastic matrix B from A by some other set of operations, other
than multiplying by diagonal matrices. This throws up a number of issues:
does a B always exist, is it unique in some sense, and if so, what is its
relationship to A?  This seems like a really hard problem.  If this problem
can be set up as an optimization problem, perhaps, then we could establish
conditions under which a solution would exist.

In the iterative proportional fitting for contingency tables, we have the
row sums = column sums = grand total, so there is no problem.  

Also, in the case of infinity norm, the constraints are much looser so
convergence is easy.

I also wonder about the "physical reality" of this problem - i.e. is there a
physical problem that can give rise to a rectangular doubly stochastic
matrix?  In the Markov chain problems, one always gets a square matrix.  I
am not familiar with graph theory applications, where doubly stochastic
matrices play a useful role.

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: Martin Maechler [mailto:maechler at stat.math.ethz.ch] 
Sent: Tuesday, October 17, 2006 3:31 PM
To: Ravi Varadhan
Cc: R-help at stat.math.ethz.ch; 'Florent Bresson'
Subject: Re: [R] Re : Generate a random bistochastic matrix

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.