Return to Answer

I didn't know Roman's trick with SparseArray. Here's another trick: use Method->"Arnoldi", asking for only the eigenvector corresponding to the smallest magnitude eigenvalue (approx. zero) 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-IdentityMatrix[4000,SparseArray],-1,Method->{"Arnoldi"}][[1]];//AbsoluteTiming
(* 0.092582 *)

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.

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.

I didn't know Roman's trick with SparseArray. Here's another trick: Method->"Arnoldi", asking for only the eigenvector corresponding to the smallest magnitude eigenvalue (approx. zero) 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-IdentityMatrix[4000,SparseArray],-1,Method->{"Arnoldi"}][[1]];//AbsoluteTiming
(* 0.092582 *)

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.

I didn't know Roman's trick with SparseArray. Here's another trick: use Method->"Arnoldi", asking for only the eigenvector corresponding to the smallest magnitude eigenvalue (approx. zero) 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-IdentityMatrix[4000,SparseArray],-1,Method->{"Arnoldi"}][[1]];//AbsoluteTiming
(* 0.092582 *)

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.

I didn't know Roman's trick with SparseArray. Here's another trick: Method->"Arnoldi" and ask for only one eigenvector with -1.

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

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

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.

I didn't know Roman's trick with SparseArray. Here's another trick: Method->"Arnoldi", asking for only the eigenvector corresponding to the smallest magnitude eigenvalue (approx. zero) 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-IdentityMatrix[4000,SparseArray],-1,Method->{"Arnoldi"}][[1]];//AbsoluteTiming
(* 0.092582 *)

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.