# Numerically stable computation of matrix power

**URL:** <https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584>\
**Category:** General Usage\
**Tags:** linearalgebra\
**Created:** [November 27, 2019, 4:20pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584 "2019-11-27T16:20:19Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 27, 2019, 4:20pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/1 "2019-11-27T16:20:20Z")

</div>

I admit the title is not great, but I couldn’t find a better one.

I have a matrix M whose eigenvalues lie in the unit disk and the only eigenvalue on the unit cycle is 1. Moreover the algebraic multiplicity of the eigenvalue 1 is equal to its geometric multiplicity, i.e. there is no Jordan block associated with it.

Let M=PJP^{-1} be the Jordan decomposition of M. The matrix J is M's Jordan form and the limit J\_\infty = \lim\_{n\to\infty}J^n is well defined. Finally I define M\_\infty=PJ\_\infty P^{-1}.

I want to compute M\_\infty numerically. The problem is that numerically computing the Jordan decomposition is not stable. Is there any numerically stable way to compute M\_\infty?

---

<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:** [November 27, 2019, 4:34pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/2 "2019-11-27T16:34:40Z")

</div>

> [@gmouts](#):
>
> I have a matrix M whose eigenvalues lie in the unit disk and the only eigenvalue on the unit cycle is 1. Moreover the algebraic multiplicity of the eigenvalue 1 is equal to its geometric multiplicity, i.e. there is no Jordan block associated with it. […] I want to compute M^\infty numerically.

~~If you know this property of the spectrum analytically, then you just need to find a basis Q for the nullspace of M-I, in which case M^\infty = QQ^\* is simply a projection onto this subspace. The function `Q = nullspace(M-I)` in the `LinearAlgebra` standard library will do this for you using the SVD.~~_Correction_: a non-orthogonal projection is needed, as described [below](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/7).

A problem still arises, of course, if your matrix has eigenvalues extremely close to `1`, such that they can’t be distinguished from `nullspace(M-I)` to machine precision.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 27, 2019, 4:40pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/3 "2019-11-27T16:40:33Z")

</div>

I think this is the solution. I suspect that I don’t have a problem with eigenvalues that are arbitrarily close to 1.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 27, 2019, 5:14pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/4 "2019-11-27T17:14:54Z")

</div>

Ah, I think that QQ^\* is not M\_\infty but J\_\infty.

---

<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:** [November 27, 2019, 5:34pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/5 "2019-11-27T17:34:13Z")

</div>

No. J\_\infty is diagonal

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 27, 2019, 6:01pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/6 "2019-11-27T18:01:25Z")

</div>

I have an example where this does not work. I wrote it in Julia because I’m not sure how to write matrices in latex here.

```julia
julia> using LinearAlgebra

julia> M = [1 0 1 0; 0 1 0 1/2; 0 0 0 1/2; 0 0 0 0]
4×4 Array{Float64,2}:
 1.0 0.0 1.0 0.0
 0.0 1.0 0.0 0.5
 0.0 0.0 0.0 0.5
 0.0 0.0 0.0 0.0

julia> A = M - I
4×4 Array{Float64,2}:
 0.0 0.0 1.0 0.0
 0.0 0.0 0.0 0.5
 0.0 0.0 -1.0 0.5
 0.0 0.0 0.0 -1.0

```

For this specific matrix M^n=M^2 for any n\>1:

```julia
julia> M*M
4×4 Array{Float64,2}:
 1.0 0.0 1.0 0.5
 0.0 1.0 0.0 0.5
 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0

julia> M*M*M
4×4 Array{Float64,2}:
 1.0 0.0 1.0 0.5
 0.0 1.0 0.0 0.5
 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0

```

The nullspace of A is straightforward:

```julia
julia> Q = nullspace(A)
4×2 Array{Float64,2}:
  0.0 1.0
 -1.0 0.0
  0.0 0.0
  0.0 0.0

```

However QQ^\* is not M^2:

```julia
julia> Q*transpose(Q)
4×4 Array{Float64,2}:
 1.0 0.0 0.0 0.0
 0.0 1.0 0.0 0.0
 0.0 0.0 0.0 0.0
 0.0 0.0 0.0 0.0

```

---

<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:** [November 27, 2019, 6:35pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/7 "2019-11-27T18:35:30Z")

</div>

Sorry, you’re quite right — the basic problem is that the nullspace of M-I is not orthogonal to the other Jordan vectors of M, so the correct projection is non-orthogonal.

I think what you need, instead, is to project using left and right nullspaces. I came up with the following a bit hastily, but it looks right to me and works for your example:

```julia
function Minf(M)
    R = nullspace(M-I)
    L = nullspace(M'-I)'
    return R / (L*R) * L
end

```

You might want to check it a bit more carefully.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 27, 2019, 6:42pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/8 "2019-11-27T18:42:54Z")

</div>

I think this is the solution. I will check it more, but in simple cases it works.

---

<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:** [November 27, 2019, 6:46pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/9 "2019-11-27T18:46:45Z")

</div>

Note that my `Minf` implementation is suboptimal (if you care about factors of 2), since really you can get the left and right nullspaces simultaneously from the left and right singular vectors of the SVD. Just look at the [source code for `nullspace`](https://github.com/JuliaLang/julia/blob/2fcb1b16fc734be41ebf348ddfe37809b081f027/stdlib/LinearAlgebra/src/dense.jl#L1400-L1409) and pull out the relevant bits.

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 28, 2019, 12:49pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/10 "2019-11-28T12:49:00Z")

</div>

Could you please explain briefly how you came up with this? Why does it work? I cannot figure it out. 😕

---

<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:** [November 28, 2019, 1:56pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/11 "2019-11-28T13:56:07Z")

</div>

> [@gmouts](#):
>
> Could you please explain briefly how you came up with this?

Let us analyze the case where M is diagonalizable, since the [defective](https://en.wikipedia.org/wiki/Defective_matrix) case of a non-diagonalizable matrix (Jordan blocks) can be obtained at the end as the limit of diagonalizable matrices.

The diagonalization is M = R\Lambda R^{-1} where R is a matrix whose columns are (right) eigenvectors. Then M^T = (R^{-1})^T \Lambda R^T, but this is also a diagonalization and so we see by inspection that (R^{-1})^T = L, a matrix whose columns are _left_ eigenvectors. (This is one proof of the general property that we can choose L^T R = I, i.e. left and right eigenvectors can be chosen orthonormal under the “unconjugated inner product” x^T y.)

That is, M = R \Lambda L^T with L^T R = I, and it follows that M^n = R \Lambda^n L^T.

In your case, all of the eigenvalues have |\lambda|\<1 except for eigenvalues with \lambda=1, so the limit \Lambda^\infty is a diagonal matrix that projects onto the rows/columns corresponding to the \lambda=1 eigenvalues. Hence M^\infty = R \Lambda^\infty L^T = R\_0 L\_0^T, where L\_0 and R\_0 are the left and right eigenvectors for \lambda=1, equivalent to left and right nullspaces \hat{L}\_0 and \hat{R}\_0 of M-I. Technically, the left nullspace is the nullspace of M^\* - I, but this conjugation is cancelled if we also use M^\infty = \hat{R}\_0 \hat{L}\_0^\* instead of the transpose.

However, we have to be careful, especially in the case where \lambda=1 has multiplicity \> 1, since there are multiple possible choices of the left and right eigenvectors/nullspaces and we have to choose them consistently to ensure that \hat{L}\_0^\* \hat{R}\_0 = I. If we choose them arbitrarily, we might in general need some change of basis \hat{R}\_0 \to \hat{R}\_0 B, giving M^\infty = \hat{R}\_0 B \hat{L}\_0^\* for some (invertible) B. We can easily determine B by the fact that M^\infty \hat{R}\_0 = \hat{R}\_0, which immediately gives B = (\hat{L}\_0^\* \hat{R}\_0)^{-1}.

So, in the end, for any diagonalizable matrix M with the spectrum you assumed, we have M^\infty = \hat{R}\_0 (\hat{L}\_0^\* \hat{R}\_0)^{-1} \hat{L}\_0^\*, in terms of left and right nullspace bases \hat{L}\_0 and \hat{R}\_0 of M-I (regardless of how they are obtained/normalized).

For a defective matrix M, as long as the \lambda=1 eigenvalue is not itself defective, if we take the defective matrix M as the limit of diagonalizable matrices (which is always possible since diagonalizable matrices are a dense subset of all square matrices), we get the same answer (the limit of the \lambda=1 eigenvectors still span the nullspaces of M-I).

PS. I guess you get this spectrum from some kind of Markov property?

---

<div class="post-metadata">

**Author:** ![gmouts](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gmouts/32/9896_2.png) [@gmouts](https://discourse.julialang.org/u/gmouts)\
**Post date:** [November 28, 2019, 6:11pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/12 "2019-11-28T18:11:04Z")

</div>

Thank you! That was very clear. Yes, it is the matrix of a random walk on a weakly connected graph where most of the vertices are transient.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [November 29, 2019, 1:57pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/13 "2019-11-29T13:57:32Z")

</div>

IIUC this would also work:

```julia
function Minf2(M, tol=nothing)
  if tol === nothing
    tol = 10*size(M,1)*eps(real(eltype(M)))
  end
  S = schur(M)
  select = abs.(S.values .- 1) .< tol
  S1 = ordschur(S, select)
  n1 = count(select)
  return S1.Z[:,1:n1] * (S1.Z[:,1:n1])'
end

```

You might want to verify that `S1.T[1:n1,1:n1]` really is close to the identity matrix, as a check on construction of `M` and choice of tolerance.

There are also interesting possibilities for iterative methods in case your graphs are large.

---

<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:** [November 29, 2019, 2:55pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/14 "2019-11-29T14:55:31Z")

</div>

> [@Ralph\_Smith](#):
>
> IIUC this would also work:

Nope. Try it on @gmout’s example [from above](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/6), for example.

---

<div class="post-metadata">

**Author:** ![Ralph\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ralph_smith/32/10344_2.png) [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)\
**Post date:** [November 29, 2019, 4:51pm UTC](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584/15 "2019-11-29T16:51:20Z")

</div>

You’re right, of course. I left out a “left” in my thinking; I’d have to do a second Schur decomposition and a linear solve like yours, so the suggested SVD approach is more efficient and just as stable AFAICT.