Derivatives and the tedious reality of optimization
The standard approach most people learn in calculus is straightforward enough: take the derivative, set it equal to zero, solve for x, then check whether that point is a maximum or minimum using the second derivative test or by comparing endpoint values. It works for clean polynomial and trigonometric functions, and you can usually grind through it in fifteen to twenty minutes if you are careful. The problem is that real functions rarely behave cleanly, and the method breaks down fast once you leave textbook problems behind. Start by confirming the domain you are working with. I cannot count the number of times a colleague handed me a function like f(x) = x*sqrt(4-x) and expected me to just differentiate and call it a day, when the domain is actually restricted to [-infinity, 4] because of the square root. Miss that restriction and your critical points include values where the function simply does not exist. Plot it first, even roughly. A quick sketch on graph paper or a five-second plot in Desmos saves you from chasing phantom extrema. Once you have your critical points, you need to determine which ones are maxima. The second derivative test is reliable but not foolproof. If f''(c) = 0, the test fails and you are back to first principles, checking sign changes of f' around c. I spent an afternoon once on a rational function where both the first and second derivatives vanished at the same critical point, and I had to do a sign chart analysis across seven intervals before I was confident about what was happening. That function was f(x) = (x^3 - 3x)/(x^2 + 1), and the point at x = 0 was a saddle-like inflection that looked like a local extremum on a coarse plot.
For functions of multiple variables, the process expands. You compute the gradient vector, set all partial derivatives to zero simultaneously, and then use the Hessian matrix to classify each critical point. The Hessian test tells you whether a point is a local maximum if the Hessian is negative definite at that point, which means all eigenvalues are negative. In two variables, that translates to checking that f_xx
0 and the determinant of the Hessian is positive. But here is the thing nobody warns you about: the Hessian test only confirms local behavior. A function can have a perfectly valid local maximum that is nowhere near the global maximum, and finding the global one requires checking boundary conditions, comparing all critical points, and sometimes accepting that you cannot analytically solve the system. I encountered a constrained optimization problem last year involving a production model where the objective was to maximize output subject to a cost constraint. The Lagrange multiplier method gave me three candidate solutions, and one of them was a saddle point that looked plausible at first glance. The real issue was that the constraint boundary itself contained the global maximum, not any interior critical point. I had to evaluate the objective function along the entire constraint curve using a parametrization, which turned a five-minute derivation into a numerical sweep. The workaround was straightforward: I substituted the constraint into the objective to reduce it to a single variable, then used a numerical optimizer to scan the resulting one-dimensional function. That cut my search time from roughly an hour of manual algebra to about twelve minutes of computation. When analytical methods fail, which they often do for higher-degree polynomials or transcendental combinations, numerical approaches become necessary. Gradient ascent is the workhorse here. You pick an initial point, compute the gradient, take a step in the direction of steepest ascent scaled by a learning rate, and iterate. The catch is that gradient ascent finds local maxima, not global ones, and the result depends heavily on where you start. Running the algorithm from five or six different starting points and keeping the best result is standard practice, though it does not guarantee you have found the global maximum.
Another numerical option is the Nelder-Mead simplex method, which does not require derivative information at all. It works by maintaining a simplex of points and iteratively replacing the worst vertex with a better one. It is robust for low-dimensional problems and handles non-smooth functions gracefully, but it can converge slowly and may stall near flat regions where the gradient provides almost no directional signal. For functions with sharp ridges or narrow valleys, you will often need a hybrid approach: use gradient-based methods for fast convergence in smooth regions, then switch to a direct search method near suspected optima to avoid getting trapped. There is also the question of computational cost. Analytical differentiation of a complicated expression can produce derivatives that are hundreds of times larger than the original function, which creates numerical instability when you evaluate them. Automatic differentiation tools like those available in PyTorch or JAX sidestep this by computing derivatives exactly up to machine precision without symbolic expansion. If you are working with a function defined by code rather than a closed-form expression, autodiff is usually the fastest path to a reliable solution, and it typically reduces the time spent debugging derivative calculations from hours to under ten minutes per function. The limitations are worth stating plainly. Analytical methods only work for functions where you can express the derivative in closed form and solve the resulting equations, which excludes most real-world objective functions. Numerical methods introduce approximation error and may miss narrow peaks if your search grid is too coarse. Global optimization is NP-hard in the general case, meaning there is no guaranteed efficient algorithm for arbitrary functions. You can use branch-and-bound methods or genetic algorithms for better global coverage, but these trade accuracy for speed and require tuning parameters that are not always obvious.
Get the Full Details
If your function is unimodal within a known interval, bracketing methods like the golden section search give you a maximum to arbitrary precision without computing derivatives, and they are considerably faster than grid search for one-dimensional problems. For multidimensional unconstrained problems with smooth objectives, L-BFGS is generally the default choice in optimization libraries because it approximates the Hessian efficiently and converges in far fewer iterations than basic gradient ascent. I usually reach for scipy.optimize.minimize with the method set to 'L-BFGS-B' when I need something that works without extensive configuration, and it handles bound constraints natively, which saves a separate step. The bottom line is that finding the maximum value of a function is not a single technique but a decision tree. Identify the domain. Check for analytical solutions if the function is simple enough. Fall back to numerical methods when it is not. Verify your result by evaluating the function at nearby points and checking boundary values. Most errors I see come from skipping the verification step and assuming the first critical point found is the answer.
