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: [EMAIL PROTECTED] Webpage: http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html ---------------------------------------------------------------------------- -------- -----Original Message----- From: Martin Maechler [mailto:[EMAIL PROTECTED] Sent: Tuesday, October 17, 2006 3:31 PM To: Ravi Varadhan Cc: R-help@stat.math.ethz.ch; 'Florent Bresson' Subject: Re: [R] Re : Generate a random bistochastic matrix Thank you, Ravi, >>>>> "Ravi" == Ravi Varadhan <[EMAIL PROTECTED]> >>>>> 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. ______________________________________________ R-help@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.