Broadcast raising a matrix to a specific power

General Usage
1

Hi, new to Julia here. I’m trying to model the accumulation of mutations in a population using a Markov chain, I’ve got a transition matrix and an initial state vector, basically what I want to do is raise the transition matrix to the power of the number of generations and then multiply it by the initial state vector and plot the values on a graph.

My code is as follows:

using Plots
μ=0.00005
T=[μ^2-2μ+1 μ-μ^2 μ-μ^2 μ^2; 0 1-μ 0 μ; 0 0 1-μ μ; 0 0 0 1]
i=[1000 0 0 0]

x = range(0, 10, length=100)
p_every = i*(Ref(T).^x)
y = p_every[1]
plot(x, y;dpi=460)

When I run this though, I get the following error:
DimensionMismatch(“matrix A has dimensions (1,4), vector B has length 100”)

  1. generic_matvecmul!(::Array{Array{Float64,2},1}, ::Char, ::Array{Int64,2}, ::Array{Array{Float64,2},1}, ::LinearAlgebra.MulAddMul{true,true,Bool,Bool})@matmul.jl:639
  2. mul!@matmul.jl:81 [inlined]
  3. mul!@matmul.jl:208 [inlined]
  4. *****(::Array{Int64,2}, ::Array{Array{Float64,2},1})@matmul.jl:51
  5. top-level scope@*[Local: 3]

I’ve been having a hard time finding information about the broadcast function online, so I was hoping someone might be able to clarify what I’m doing wrong.

2

Broadcasting “peels” one layer of containers and broadcasts over it (kinda like implicit list/array/set/etc comprehension). You use it corectly to broadcast over the Ref container and the x vector, i.e. for each element in the Ref and for each element in the vector, perform the ^ operation. The result of the Ref * vector will be a vector. Then you ask for the * multiplication between the i 1xn matrix and the result of the previous operation, which happens to be a vector. There is no multiplication defined between 1xn matrix and vector. You probably want to broadcast put i in another container and broadcast over that as well as Ref(i) .* Ref(T).^x.

One source of mistakes here is that people usually think of .* as elementwise multiplication. You want to actually broadcast matrix multiplication over containers that themselves contain matrices.

On a side note, matrix-vector multiplication is much faster than matrix-matrix multiplication (quadratic vs qubic time complexity). You probably will get better results with something that repeatedly does matrix-vector multiplication:

p_every = [i] # maybe add a type to this
for iter in 1:iterations
    push!(p_every, T*p_every[end])
end

or if you need only the first component

y = zeros(iterations)
buffer = copy(i)
for iter in 1:iterations
    buffer = T*buffer # this can be done more efficiently with in-place ops
    y[iter] = buffer[1]
end

or you can use the reduce foldr/foldl function or a few other ways, depending on what is important for your use case.

Two more notes:

y=p_every[1] would just give you the initial i, not the first scalar in each vector (which is what I assume you want)

range(0,10,length=100) would include fractions, not only integers. Fractional power of a matrix would be quite a bit more expensive. Make sure you actually need it (and if you actually need it, it would be more efficient to do some factorization of the matrix first, to make the computation cheaper)

5 Likes
3

Thank you so much, this is super helpful!

Yeah the efficiency of this is not something I’m at all considering because 1. this is basically my first time using julia, and 2. I just want to get a quick sketch of this problem, but I really do appreciate the insight and it’s certainly something I’ll look more into in the future.

Again thanks!

1 Like