With Newton's method, you'll be solving Hx=g (H = hessian of f, g = gradient of f) at each iteration. For large number of variables N, building H is of order N^2 and solving Hx = g is of order N^3 with an usual solver. N^2 and N^3 are really large for an already large N. I believe the reason is as simple as that. It isn't that it is tedious and difficult to write down the formulas or the code to compute H. It's just too costly, computationally speaking. There is also an increased memory cost (possibly having to store H).

When people do have ways to go around this problem, they do use Newton's method for large scale problems.

Adding to sibling post from phao, Optimizers like ADAM essentially do a diagonal approximation to Newton's method, which gives some, though certainly not all, of the benefit of Newton.

There is the seminal paper "Deep Learning via Hessian-free Optimization" [1] that applies Newton's method (Hessian free because of a clever trick of solving Hx=-g iteratively without explicitly constructing the Hessian, and getting an approximation without full O(n^3)) - but n^2 or n^3 is just too large in practice when the diagonal approximations work just fine.

[1] https://www.cs.toronto.edu/~jmartens/docs/Deep_HessianFree.p...

Adding to the answers already given (complexity related to the computation of full Hessians). First order methods, i.e., methods that only require gradients, are the methods of choice when extreme solution accuracy is not required, which is typically the case in ML settings.

One of Newton's method major selling points is that once it's close to a local minimum, under the right smoothness assumptions it essentially converges faster and faster to an exact optimizer - "One more correct significant digit each iteration". It's called Newton's method region of quadratic convergence, see [0] Theorem 14.1, p.3

[0] https://www.stat.cmu.edu/~ryantibs/convexopt-F16/scribes/new...