The generalization of matrix product states (MPS) to continuous systems, as proposed in the breakthrough Letter of Verstraete and Cirac [Phys. Rev. Lett. 104, 190405 (2010).PRLTAO0031-900710.1103/PhysRevLett.104.190405], provides a powerful variational ansatz for the ground state of strongly interacting quantum field theories in one spatial dimension. A continuous MPS (cMPS) approximation to the ground state can be obtained by simulating a Euclidean time evolution. In this Letter we propose a cMPS optimization algorithm based instead on energy minimization by gradient methods and demonstrate its performance by applying it to the Lieb-Liniger model (an integrable model of an interacting bosonic field) directly in the thermodynamic limit. We observe a very significant computational speed-up, of more than 2 orders of magnitude, with respect to simulating a Euclidean time evolution. As a result, a much larger cMPS bond dimension D can be reached (e.g., D=256 with moderate computational resources), thus helping unlock the full potential of the cMPS representation for ground state studies.