Matrix-free methods will work fine on problems with thousands of variables. Specifically, a properly implemented Newton-CG or trust-region code does not need to form a dense Hessian in order to compute a Newton iteration, but rather requires only the action of the Hessian on a vector. These methods rely on a truncated Krylov solver such as truncated-CG or truncated-MINRES in order to partially solve the linear systems. In my experience, a matrix-free, trust-region algorithm using truncated-CG works the best. It tends to work better than Newton-CG, which is a line-search method, because the trust-region helps determine how many Krylov iterations are required per iteration. In Newton-CG, you have to guess and set this number a priori. For a short description of Newton-CG, see Nocedal and Wright's book Numerical Optimization on page 169. That book also contains a description of a trust-region algorithm on page 69. Use this with truncated-CG, which is described on page 171. In addition, the big blue Trust Region Methods book has the most complete description of the trust-region algorithm. That book is a little overwhelming. Start with the basic trust-region algorithm on page 116. In step 2, use truncated-CG, which is described on page 202.
Now, these methods are not perfect. In general, you'll get better performance than a first-order method for a variety of reasons, but a high-order convergence rate is only guaranteed when close to the actual solution and only if we solve the Newton system accurately enough. If we're using a truncated Krylov solver such as truncated-CG, this is dictated by the clustering of the eigenvalues of the Hessian, which we likely can't compute. However, like I said above, in general it works well. If you want me to be frank, I've never seen a situation where a quasi-Newton method worked better than a good matrix-free optimization algorithm. I'm sure it happens, but not that often.
Finally, in terms of a warm-start, the only real issue with warm-starts would come into play if you had inequality constraints present in the problem. In unconstrained optimization, I'm not aware of any algorithm that has a problem with warm-starts. Though, there's a bit of a numerical conundrum if you actually start at the optimal solution. In any case, inequality constrained algorithms like primal-dual interior-point methods can have trouble with warm starts. However, other inequality constrained algorithms such as the Coleman-Li reflective Newton algorithm have zero problems with this. If you're close to the solution, just ignore the reflective bits of the code and it works great. As a note, good equality constrained algorithms also have no issues with warm-starts. These are things like sequential quadratic programming (SQP), or better, composite step SQP methods.
Finally, all of these algorithms are already implemented and freely available in Optizelle. It's BSD licensed and has hooks for C++, Python, and MATLAB. If you want a limited memory quasi-Newton method, it also has BFGS and SR1. Certainly, you can implement these algorithms yourself. And, to be honest, the codes for unconstrained optimization aren't that bad. However, Optizelle already does it and has been used on problems with upwards of half a billion variables.