Raising a matrix to a power in Julia

Question

I'm writing code in Julia that involves raising a large matrix of integers to a high power and I want to make this code more efficient. I've been searching around on JuliaLang, but I'm not sure whether when I raise a matrix to a power in Julia, Julia will automatically use the fastest method available (binary exponentiation or something similar) or whether it will multiply the matrices sequentially, e.g. A^p = A* A * ... * A. Could I speed up my code by manually implementing binary exponentiation, or is Julia already doing this for me?

Answer 1

Julia provides all the introspection methods you need to figure this out. Since the base library is open source and almost entirely written in Julia, it's fairly easy to see. Just take a look at:

julia> A = rand(1:10, 4, 4); p = 3;

julia> @less A^p
function (^)(A::AbstractMatrix{T}, p::Integer) where T<:Integer
    # make sure that e.g. [1 1;1 0]^big(3)
    # gets promotes in a similar way as 2^big(3)
    TT = promote_op(^, T, typeof(p))
    return power_by_squaring(convert(AbstractMatrix{TT}, A), p)
end

So it's using the internal power_by_squaring function to do its work:

julia> @less Base.power_by_squaring(A, p)
          "(e.g., [2.0 1.0;1.0 0.0]^$p instead ",
          "of [2 1;1 0]^$p), or write float(x)^$p or Rational.(x)^$p")))
function power_by_squaring(x_, p::Integer)
    x = to_power_type(x_)
    # … skipping the obvious branches to find …
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end

Binary exponentiation! Now this doesn't re-use any temporaries, so you might be able to do a bit better with a judicious use of in-place mul!.