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.
