Hi all. I’m quite new to Julia and this is my first post.
I want to optimize a functional F[\phi] = \int_0^1 \left[\frac{1}{2}\left(\partial_x \phi\right)^2 + V(x,\phi(x))\right]\text{d}x. Of course, this can be done by using variational methods and solving ODEs, but I want to treat it as an optimization problem. I know \phi is periodic in the internal [0,1], so I write it as an expansion \phi = \sum_n \phi_n e^{i 2\pi n x} (or equivalently in terms of sines and cosines) using the ApproxFun package. The functional now becomes a function of the series coefficients: F[\phi] \to F(\phi_1,\ldots ,\phi_n), which can be minimized using optimization methods.
I have a code which is working, but slowly. It is not so bad in the example provided, but will be much slower when the number of coefficients needed becomes very large.
I suspect the problem is to do with memory allocations. Defining a series expansion using ApproxFun is fine and so is differentiating the function, but when I do other operations like squaring the function, it does a lot of memory allocations. This propagates quite poorly through the numerical integration and optimization functions.
I would really appreciate if somebody could help me find a way to get this method to work without the huge memory allocation (which will hopefully speed things up and make it scale better when going to 2D). Minimal working example below using V(x,\phi(x)) = k \cos{(2\pi(x+\phi(x)))} as an example:
using ApproxFun # For defining series expansions
using QuadGK # for numerical integration
using Optim # for optimization
using Plots ; gr()
using LaTeXStrings
s=Fourier(0..1) # Specify the space.
# CosSpace / SinSpace for just cos/sin
# or exponentials with Laurent(0..1)
# Returns an expansion using entries of vec as coefficients
function fourier(space,vec::Vector{Float64})
return Fun(space, vcat([0.0],vec)) # fix const term to 0.
end;
function V(space,vec::Vector{Float64},k::Float64) # Just the potential
return x->k*cos.(2π*(x+fourier(space,vec).(x)))
end
# The integrand
function integrand(space, vec::Vector{Float64}, k::Float64)
return x->0.5*(fourier(space,vec)'.(x))^2 + k*cos(2π*(x+fourier(space,vec)(x)))
end
# Do the integral numerically for a given vector of coefficients.
function integral(space,vec::Vector{Float64}, k::Float64)
return quadgk(integrand(space,vec,k), 0, 1, atol=1e-12)
end
Now I can use optim to determine \phi(x) like this:
@time res=optimize(v->integral(s,v,10.0)[1], zeros(20));
sol = fourier(s,res.minimizer)
# setting the initial condition to zeros(20) also sets the length of v
37.462515 seconds (702.81 M allocations: 29.900 GiB, 10.96% gc time, 0.05% compilation time)
This gives the correct answer, which I have also obtained by solving the ODE. But the number of memory allocations is alarming, and I suspect this would become problematic in situations where I need a large number of coefficients.
I checked the timing / memory allocations of some functions:
v1 = [0.0,1.0]
@time fourier(s,v1) # sin(2πx)
@time fourier(s,v1)' # 2π*cos(2πx)
@time (fourier(s,v1)')^2 # 4π²cos²(2πx)
0.000005 seconds (3 allocations: 256 bytes)
0.000008 seconds (6 allocations: 528 bytes)
0.000238 seconds (378 allocations: 19.734 KiB)
So the problem comes when I do operations on the Fun() objects from ApproxFun, e.g. multiplying or taking cos(Fun()). Obviously this gets much worse when I go to do the numerical integration:
@time quadgk((fourier(s,v1)')^2, 0, 1, atol=1e-12)
@time quadgk(V(s,vector,1.0), 0, 1, atol=1e-12) # k = 1
@time integral(s,vector,1.0)
0.000676 seconds (9.06 k allocations: 392.062 KiB)
0.004499 seconds (85.10 k allocations: 3.697 MiB)
0.011349 seconds (190.62 k allocations: 8.405 MiB)
And worse again when I pass to the optimizer. So any help would be appreciated (as well as any other general advice on improving this code).
PS: I know that in this example, \phi(x) is odd, which cuts the number of coefficients in half and speeds up the program. But in the actual problem I want to solve it is not odd, so I don’t do that here.