I need help optimising a piece of code which involves the computation of many small vectors and matrices. I have a list of vectors and matrices, given by
S_i respectively. I need to calculate
nu_i = A*mu_i + b and
P_i = A*S_i*A' + R, where
R are all the correct size. This operation is typical in multi-target tracking, where a target’s spatial statistics are put through a motion (or measurement) model. This operation is completed at every time-step of simulation and it helps if it is fast.
In the code below I consider two different options: a for loop and a “vectorised” approach I previously used in Matlab. The “vectorised” approach is faster for large lists of vectors and matrices, but I’m not sure if I’m using Julia effectively. Is there anything I can do to speed up this code? The operation can be computed in parallel, but I’ve not had much success with Julia’s distributed for loop.
using LinearAlgebra # De-vectorised code function devectorised!(mu, S, A, ATranspose, R, simulationLength) n = size(mu, 2) for i = 1:simulationLength for j = 1:n mu[:, j] = A*mu[:, j] + u S[:, :, j] = A*S[:, :, j]*ATranspose + R end end end # Vectorised code function vectorised(nu, P, A, ATranspose, R, simulationLength) (d, n) = size(nu) mu = copy(nu) S = copy(P) for i = 1:simulationLength mu = A*mu .+ u D = reshape(permutedims(S, (1, 3, 2)), (d*n, d)) E = reshape(D*ATranspose, (d, n, d)) F = reshape(permutedims(E, (1, 3, 2)), (d, d*n)) S = reshape(A*F, (d, d, n)) .+ R end return mu, S end # Admin d = 6 n = 1000 simulationLength = 3000 # Motion model T = 20e-3 A = [Array(1.0I, 3, 3) T*Array(1.0I, 3, 3); zeros(3, 3) Array(1.0I, 3, 3)] ATranspose = A' u = zeros(d, 1) u[3, 1] = -0.5*(9.81)*(T^2) u[6, 1] = -(9.81)*(T) # Process noise R = 4*Array(1.0I, d, d) # Time the functions mu = rand(d, n) S = rand(d, d, n) nu = copy(mu) P = copy(S) @time devectorised!(mu, S, A, ATranspose, R, simulationLength) @time (nuUpdated, PUpdated) = vectorised(nu, P, A, ATranspose, R, simulationLength) # Check answers maximumMeanError = maximum(abs.(mu[:] - nuUpdated[:])) maximumCovarianceError = maximum(abs.(S[:] - PUpdated[:])) println("The maximum mean error is ", maximumMeanError) println("The maximum covariance error is ", maximumCovarianceError)
Furthermore, is Julia appropriate for real-time operation? I’ve previously used heavily vectorised Matlab code for real-time operation and it worked well enough; Julia should perform just as well.