[R] Generating a stochastic matrix with a specified second dominant eigenvalue
Albyn Jones
jones at reed.edu
Fri Oct 16 19:18:05 CEST 2009
ooops! I tested it on a 3x3, and a 4x4, and it worked both times.
Then I ran it once more and pasted it in without looking carefully...
Here's another half-baked idea. Generate a random graph with given
degree distribution, turn the incidence matrix into a stochastic
matrix (ie the transition probability is 1/d where d is the degree of
the vertex). Now all we need is the connection between the degree
distribution and the second eigenvalue :-) If you want lambda
close to one, perhaps you could generate two graphs, then
add vertices joining the two graphs until the second eigenvalue
is the right size.... of course, this will give a random lambda,
not exactly your chosen value.
albyn
On Fri, Oct 16, 2009 at 10:32:48AM -0400, Ravi Varadhan wrote:
> A valiant attempt, Albyn!
>
> Unfortunately, the matrix B is not guaranteed to be a stochastic matrix. In
> fact, it is not even guaranteed to be a real matrix. Your procedure can
> generate a B that contains negative elements or even complex elements.
>
> > M = matrix(runif(9),nrow=3)
> > M = M/apply(M,1,sum)
> > e=eigen(M)
> > e$values[2]= .7
> > Q = e$vectors
> > Qi = solve(Q)
> > B = Q %*% diag(e$values) %*% Qi
> >
> > eigen(B)$values
> [1] 1.00000000 0.70000000 -0.04436574
> > apply(B,1,sum)
> [1] 1 1 1
> > B
> [,1] [,2] [,3]
> [1,] 0.77737077 0.3340768 -0.1114476
> [2,] 0.20606226 0.2601840 0.5337537
> [3,] 0.08326022 0.2986603 0.6180794
> >
>
> Note that B[1,3] is negative.
>
> Another example:
>
> > M = matrix(runif(9),nrow=3)
> > M = M/apply(M,1,sum)
> > e=eigen(M)
> > e$values[2]= .95
> > Q = e$vectors
> > Qi = solve(Q)
> > B = Q %*% diag(e$values) %*% Qi
> >
> > eigen(B)$values
> [1] 1.00000000-0.00000000i 0.95000000+0.00000000i -0.09348883-0.02904173i
> > apply(B,1,sum)
> [1] 1+0i 1-0i 1+0i
> > B
> [,1] [,2] [,3]
> [1,] 0.6558652-0.550613i 0.2408879+0.2212234i 0.1032469+0.3293896i
> [2,] 0.1683119+1.594515i 0.6954317-0.7378503i 0.1362564-0.8566647i
> [3,] 0.2812210-2.462385i 0.2135648+1.2029636i 0.5052143+1.2594216i
> >
>
> Note that B has complex elements.
>
> So, I took your algorithm and embedded it in an iterative procedure to keep
> repeating your steps until it found a B matrix that is real and
> non-negative. Here is that function:
>
> e2stochMat <- function(N, e2, maxiter) {
> iter <- 0
>
> while (iter <= maxiter) {
> iter <- iter + 1
> M <- matrix(runif(N*N), nrow=N)
> M <- M / apply(M,1,sum)
> e <- eigen(M)
> e$values[2] <-e2
> Q <- e$vectors
> B <- Q %*% diag(e$values) %*% solve(Q)
> real <- all (abs(Im(B)) < 1.e-16)
> positive <- all (Re(B) > 0)
> if (real & positive) break
> }
> list(stochMat=B, iter=iter)
> }
>
> e2stochMat(N=3, e2=0.95, maxiter=10000) # works
> e2stochMat(N=5, e2=0.95, maxiter=10000) # fails
>
> This works for very small N, say, N <= 3, but it fails for larger N. The
> probability of success is a decreasing function of N and e2. So, the
> algorithm fails for large N and for values of e2 close to 1.
>
> Thanks for trying.
>
> 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_personal_pages/Varadhan.h
> tml
>
>
>
> ----------------------------------------------------------------------------
> --------
>
>
> -----Original Message-----
> From: Albyn Jones [mailto:jones at reed.edu]
> Sent: Thursday, October 15, 2009 6:56 PM
> To: Ravi Varadhan
> Cc: r-help at r-project.org
> Subject: Re: [R] Generating a stochastic matrix with a specified second
> dominant eigenvalue
>
> I just tried the following shot in the dark:
>
> generate an N by N stochastic matrix, M. I used
>
> M = matrix(runif(9),nrow=3)
> M = M/apply(M,1,sum)
> e=eigen(M)
> e$values[2]= .7 (pick your favorite lambda, you may need to fiddle
> with the others to guarantee this is second largest.)
> Q = e$vectors
> Qi = solve(Q)
> B = Q %*% diag(e$values) %*% Qi
>
> > eigen(B)$values
> [1] 1.00000000 0.70000000 -0.08518772
> > apply(B,1,sum)
> [1] 1 1 1
>
> I haven't proven that this must work, but it seems to. Since you can
> verify that it worked afterwards, perhaps the proof is in the pudding.
>
> albyn
>
>
>
> On Thu, Oct 15, 2009 at 06:24:20PM -0400, Ravi Varadhan wrote:
> > Hi,
> >
> >
> >
> > Given a positive integer N, and a real number \lambda such that 0 <
> \lambda
> > < 1, I would like to generate an N by N stochastic matrix (a matrix with
> > all the rows summing to 1), such that it has the second largest eigenvalue
> > equal to \lambda (Note: the dominant eigenvalue of a stochastic matrix is
> > 1).
> >
> >
> >
> > I don't care what the other eigenvalues are. The second eigenvalue is
> > important in that it governs the rate at which the random process given by
> > the stochastic matrix converges to its stationary distribution.
> >
> >
> >
> > Does anyone know of an algorithm to do this?
> >
> >
> >
> > Thanks for any help,
> >
> > 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_personal_pages/Varadhan.
> > html>
> >
> http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
> > tml
> >
> >
> >
> >
> ----------------------------------------------------------------------------
> > --------
> >
> >
> > [[alternative HTML version deleted]]
> >
> > ______________________________________________
> > R-help at r-project.org 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.
> >
>
>
More information about the R-help
mailing list