I have implemented quite a wide variety of non-linear solvers on the GPU, including LBFGS, Barzilai Borwein gradient descent and non-linear conjugate gradient.
For this, the non-linear conjugate gradient of Dai & Yuan has been the most efficient. In general, other version of the nonlinear conjugate gradient may be more efficient (such as CG-DESCENT), but can also be trickier to implement.
LBFGS is in general a very solid choice, and unless you're really strapped for memory it's probably the best place to start.
Both the conjugate gradient and BFGS require line searches though, which is where the fp32 becomes a problem. Rather than using the standard Wolfe conditions for the line search, I would suggest using the approximate Wolfe condition suggested here. The paper is a little involved, but the important stuff is equation 4.1. Essentially they explicitly introduce the precision with which you can calculate your function.
Considerations for the GPU:
You have a lot of small problems, which is slightly different from my use case of one large problem. Consider running 1 problem per GPU block (or warp, rather) if you can parallelize function and gradient evaluations to use all the threads in a block. That way it's not a problem if different problems require a different number of iterations.
If this is not an option, I would go with the LBFGS solver. If your function is well behaved, you might get away with simply using a step size of 1 (avoiding the line search) and just running all the problems for a fixed number of iterations.