[R] Generating a stochastic matrix with a specified second dominant eigenvalue

[Ravi Varadhan]

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.
>