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

and `S_i`

respectively. I need to calculate `nu_i = A*mu_i + b`

and `P_i = A*S_i*A' + R`

, where `A`

, `b`

and `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.