# Stationary distribution with sparse transition matrix

**URL:** <https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301>\
**Category:** Finance and Economics\
**Created:** [May 27, 2020, 11:22pm UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301 "2020-05-27T23:22:13Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [May 27, 2020, 11:22pm UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/1 "2020-05-27T23:22:13Z")

</div>

Dear All,

I would like to ask for help/advice with following problem. I am solving the Hugget type of model (stationary distribution/transition) using the transition matrix approach (one AR(1) and one IID shock). In the first version, I solved for stationary distribution just by brute-force (rising the matrix to power 3000). Unfortunately, this approach becomes very costly for finer discretizations of shock processes, when the size of matrix grows to (7000\*7000). Also classical eigendecomposition using LinearAlgebra.eigen is pretty slow (it solves for all eigenvalues).

I would like to ask which method of solving for stationary distribution (finding eigenvector of the transpose of the transition matrix associated with unit eigenvalue) would be most efficient, given following sparsity pattern of the transition matrix.

Any advice/suggestion would be welcome!

Best,  
Honza

 ![Sparsity](https://global.discourse-cdn.com/julialang/original/3X/f/0/f026dd1157f5e569f5a1d0858d8a79ca91426cf0.png)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 28, 2020, 12:48am UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/2 "2020-05-28T00:48:44Z")

</div>

Since you know the eigenvalue, solving for an eigenvector is just a linear solve. If you set one of the eigenvector components to 1, you get a sparse non singular system of smaller size (by 1) that you can solve with `\`

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 28, 2020, 12:59am UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/3 "2020-05-28T00:59:51Z")

</div>

For example:

```julia
using SparseArrays, LinearAlgebra

# construct a random 10000×10000 sparse Markov matrix M,
# whose columns sum to 1 (like the transpose of your transition matrix)
A = sprand(10000,10000, 1e-2)
M = A ./ sum(A, dims=1)

# find an eigenvector x corresponding to eigenvalue λ = 1,
# normalized so that the first component is 1
x = [1; (I - M[2:end,2:end]) \ Vector(M[2:end,1])]

```

You can easily check that this indeed produces an eigenvector:

```julia
julia> norm(M * x - x) / norm(x)
3.419954054123219e-16

```

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [May 28, 2020, 1:16am UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/4 "2020-05-28T01:16:54Z")

</div>

Isn’t the Arnoldi method more efficient in this case?

I get ~16 seconds for your `\` code, and 2.3 seconds for Arnoldi:

```julia
julia> using ArnoldiMethod

julia> A = sprand(10000,10000, 1e-2);

julia> M = A ./ sum(A, dims=1);

julia> @time decomp, history = partialschur(M, nev=1);
  2.044695 seconds (154 allocations: 3.144 MiB)

julia> decomp.eigenvalues
1-element Array{Complex{Float64},1}:
 1.0000000000000002 + 0.0im

julia> norm(M * decomp.Q - decomp.Q)
4.4199549104149377e-16

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 28, 2020, 1:53am UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/5 "2020-05-28T01:53:13Z")

</div>

> [@dpsanders](#):
>
> Isn’t the Arnoldi method more efficient in this case?

Depends.

- The `\` operation (sparse-direct solve) will be likely be faster for a sparsity pattern that is regular rather than random, like the actual one shown in the original post, and in general how fast it is depends on the sparsity pattern.

- The speed of Arnoldi depends what the spectrum of the actual matrix looks like — if there are many eigenvalues close to 1, then Arnoldi could take a long time to converge.

---

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [May 28, 2020, 9:26am UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/6 "2020-05-28T09:26:34Z")

</div>

Thank you very much!

So, I doesn’t have to declare some particular sparsity pattern to use sparse matrices rutines efficiently? My matrix unfortunately doesn’t fit some pattern exactly…

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 28, 2020, 12:41pm UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/7 "2020-05-28T12:41:55Z")

</div>

> [@Honza9723](#):
>
> So, I doesn’t have to declare some particular sparsity pattern to use sparse matrices rutines efficiently?

No, you don’t have to declare anything, and the matrices don’t have to have an _exactly_ regular pattern.

My point was simply that sparse-direct algorithms (for `\`) tend to be more efficient for non-random sparsity patterns, and in general their efficiency varies for matrices with different patterns of nonzeros even if the number of nonzeros is the same.

---

<div class="post-metadata">

**Author:** ![Honza9723](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/honza9723/32/7807_2.png) [@Honza9723](https://discourse.julialang.org/u/Honza9723)\
**Post date:** [May 28, 2020, 9:37pm UTC](https://discourse.julialang.org/t/stationary-distribution-with-sparse-transition-matrix/40301/8 "2020-05-28T21:37:54Z")

</div>

Thank you very much!

I get 170x speed up vis-á-vis my original “brute force” solution method.

Best,  
Honza