CODEWITHSUNDEEP
matrixjulialinear-algebramatrix-multiplication
Magdalen Ruth
Magdalen Ruth

Reputation: 45

Raising a matrix to a power in Julia

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?

Upvotes: 3

Views: 4297

Answers (2)

Oscar Smith
Oscar Smith

Reputation: 6378

One thing worth pointing out is that if your matrices are square the powers are very large (>50), you will probably save time by transforming to Jordan Normal Form, raising, and transforming back. This link gives the basics from the math side https://math.stackexchange.com/questions/354277/square-matrix-multiplication-when-raised-to-a-power

Upvotes: 3

mbauman
mbauman

Reputation: 31342

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

Upvotes: 4

Related Questions