# 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:** 1\
**Showing post:** 6

<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

```

---

_[View the full topic](https://discourse.julialang.org/t/numerically-stable-computation-of-matrix-power/31584)._