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: rvarad...@jhmi.edu
Webpage:
http://www.jhsph.edu/agingandhealth/People/Faculty_personal_pages/Varadhan.h
tml
----------------------------------------------------------------------------
--------
-----Original Message-----
From: Albyn Jones [mailto:jo...@reed.edu]
Sent: Thursday, October 15, 2009 6:56 PM
To: Ravi Varadhan
Cc: r-help@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: rvarad...@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@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.
>
______________________________________________
R-help@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.
