I have a single variable function g(x) where x(t) is a periodic function. I am interested in obtaining the jacobian of the Fourier harmonics of g when differentiated with respect to the harmonics of x(t).
It is well known that, when working in complex variables, this jacobian has a Toeplitz structure. However, I am interested in obtaining the jacobian when working with the real Fourier coefficients. The following function exploits the structure of the jacobian, where only the first column is computed by finite differences and then the rest of the jacobian follows from it.
Nevertheless, the performance gain from this function is very minor when compared to computing the whole jacobian entirely by finite differences. What could I improve to increase the speed?
function toeplitz_jacobian(f, x)
n = length(x)
h = 10^-8
x_pert = x + vcat(h, zeros(n - 1))
f_x = f(x)
f_x_pert = f(x_pert)
# Compute first column of the jacobian
∂f∂x⁰ = (f_x_pert - f_x) / h
J = zeros(n, n)
# Exploit toeplitz structure to generate the rest of the jacobian
@views begin
J[:, 1] .= ∂f∂x⁰[1:n]
J[1, 2:end] .= ∂f∂x⁰[2:n] .* 2
∂f∂x⁰ = vcat(∂f∂x⁰[1], 0, ∂f∂x⁰[2:end])
end
@views for k in 2:2:(n - 1)
for j in k:2:(n - 1)
J[k, j] = ∂f∂x⁰[j - k + 1] + ∂f∂x⁰[j + k + 1]
J[k + 1, j + 1] = ∂f∂x⁰[j - k + 1] - ∂f∂x⁰[j + k + 1]
J[k, j + 1] = ∂f∂x⁰[j + k + 2] + ∂f∂x⁰[j - k + 2]
J[k + 1, j] = ∂f∂x⁰[j + k + 2] - ∂f∂x⁰[j - k + 2]
if j > k
J[j, k] = J[k, j]
J[j + 1, k + 1] = J[k + 1, j + 1]
J[j, k + 1] = J[k + 1, j]
J[j + 1, k] = J[k, j + 1]
end
end
end
return J
end