4
$\begingroup$

I am doing some simulations for simple random walks on directed random graphs. From a graph of n vertices, I get a n by n transition probability matrix transit, with 2 n non-zero entries. I would like to compute the stationary distribution pi (a vector of n non-negative entries which sum to one), such that pi * transit == pi.

Currently, I found two ways to do this. The first is to use NullSpace.

eigenVector=First[NullSpace[N[transit]-IdentityMatrix[n]]];
pi=eigenVector/Total[eigenVector]

The second is to define a Markov Process directly on the underline graph g and use StationaryDistribution as follows.

n=VertexCount[g];
mp=DiscreteMarkovProcess[1,g];
stationary=StationaryDistribution[mp];
pi=NProbability[x==#,x\[Distributed]stationary]&/@Range[n]

Both methods can handle about n = 2000 on my laptop. But I would really like to compute this for a bit higher n. Any suggestions?


Update: For example g and transit with 1000, 2000, 3000 and 4000 nodes, see this link.

$\endgroup$
6
  • $\begingroup$ I’ve got some ideas but can you supply examples of transit and g? $\endgroup$
    – Chris K
    Nov 13, 2020 at 14:29
  • $\begingroup$ @ChrisK Thanks. I have prepared a .wl file. Please see the edit. $\endgroup$
    – faceclean
    Nov 13, 2020 at 14:41
  • $\begingroup$ When you compute the null space, is the matrix comprised of exact of approximate numbers? At the size range in question, I would expect serious memory issues unless it is machine reals (or complexes, but that won't be the case in this particular scenario). $\endgroup$ Nov 13, 2020 at 16:24
  • $\begingroup$ @DanielLichtblau I actually turns transit to numeric values before calling NullSpace. I have updated the code. $\endgroup$
    – faceclean
    Nov 13, 2020 at 17:14
  • $\begingroup$ Could you also give a smaller example so that we can compare the results against your methods? $\endgroup$
    – Chris K
    Nov 14, 2020 at 0:53

2 Answers 2

7
$\begingroup$

You may run into memory problems because you are constructing a full matrix with IdentityMatrix[n] instead of a sparse one with IdentityMatrix[n, SparseArray]. By using a sparse identity matrix you can do 4000x4000 in about 8 seconds, without excessive memory use:

Dimensions[transit]
(*    {4000, 4000}    *)

eigenVector = First[NullSpace[
  N[transit - IdentityMatrix[Dimensions[transit], SparseArray]]]]; // AbsoluteTiming // First
(*    8.32465    *)

Just make sure you keep all matrices sparse and never construct a full ("normal") matrix.

$\endgroup$
3
  • $\begingroup$ Hi, thanks for pointing this out. With this improvement, I can now deal with 4000x4000 matrix, though with as much as 60 seconds time. I guess the big time difference is because there is a of problem with the examples I provided. I have updated them now. Can you try the 4000x4000 one again? $\endgroup$
    – faceclean
    Nov 14, 2020 at 12:42
  • 1
    $\begingroup$ Even after your update it still takes 8 seconds on my laptop. Maybe your computer is just slow. $\endgroup$
    – Roman
    Nov 14, 2020 at 14:43
  • $\begingroup$ Yeah, it could just be my laptop is bit slow. Thanks. I am happy with the current performance anyway. $\endgroup$
    – faceclean
    Nov 14, 2020 at 14:58
6
$\begingroup$

edit: it's even simpler since this is a discrete-time Markov chain

I didn't know Roman's trick with SparseArray. Here's another trick: use Method->"Arnoldi", asking for only the eigenvector corresponding to the largest eigenvalue (approx. one) using 1, which is the stationary distribution.

evNull=First[NullSpace[N[transit]-IdentityMatrix[4000,SparseArray]]];//AbsoluteTiming
(* 15.6608 -- I must need a faster computer! *)

evArnoldi=Eigenvectors[N@transit,1,Method->{"Arnoldi"}][[1]];//AbsoluteTiming
(* 0.024535 *)

So there's another couple of orders-of-magnitude speedup for you. I'm not patient enough to wait for StationaryDistribution to finish!

The eigenvectors look the same when plotted BTW.

$\endgroup$
6
  • 2
    $\begingroup$ In this case it would be advantageous to point the Arnoldi algorithm towards the "right" eigenvalue by using a Shift operation: see this answer. In the present example you got lucky to catch the right eigenvalue. $\endgroup$
    – Roman
    Nov 14, 2020 at 16:46
  • 1
    $\begingroup$ @Roman True in general, but in this case I don't think it's luck, because all eigenvalues should be nonpositive except the one associated with the stationary distribution, which should be zero (conservation of probability). See cims.nyu.edu/~holmes/teaching/asa19/handout_Lecture4_2019.pdf for that last point. At least I hope this is true, because we're using this technique in a paper now :) $\endgroup$
    – Chris K
    Nov 14, 2020 at 17:17
  • $\begingroup$ I should add, I'm a bit unclear on OP's actual problem, whether it is a continuous- or discrete-time process and the role of the -IdentityMatrix[n] part here. $\endgroup$
    – Chris K
    Nov 14, 2020 at 17:32
  • $\begingroup$ Agreed, you're right that if the matrix is constructed correctly, it will have the right spectral property. $\endgroup$
    – Roman
    Nov 14, 2020 at 17:35
  • 1
    $\begingroup$ @Roman Upon closer reading of the problem, it's clear that it's a discrete-time Markov chain, so the solution is even easier (and doesn't require IdentityMatrix at all). $\endgroup$
    – Chris K
    Nov 14, 2020 at 17:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.