Numerically stable computation of matrix power

General Usage
6

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> 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> 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> 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> 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