SciPy - Optimize

SciPy - Optimization (scipy.optimize)

The scipy.optimize package provides a rich set of algorithms for function minimization (optimization), root finding, and curve fitting. Optimization is the process of finding the input values (parameters) that result in the minimum (or maximum) value of a given function.

This module includes tools for:

• Scalar function minimization: Finding the minimum of a function of one variable (minimize_scalar).
• Multivariate function minimization: Finding the minimum of a function of multiple variables, both unconstrained and constrained (minimize). Various algorithms are available (e.g., BFGS, Nelder-Mead, SLSQP).
• Global optimization: Routines designed to find the global minimum, potentially avoiding local minima (e.g., basinhopping, differential_evolution).
• Least-squares minimization: Minimizing the sum of squares of residuals, often used in curve fitting (least_squares, curve_fit).
• Root finding: Solving equations of the form $f(x) = 0$. For scalar functions (newton, brentq) and multivariate systems (root).

Minimizing Multivariate Scalar Functions with minimize()

The minimize() function offers a unified interface to various algorithms for minimizing scalar functions of one or more variables. Let’s consider the Rosenbrock function, a common benchmark for optimization algorithms. For N variables, it’s defined as:

$$f(\mathbf{x}) = \sum_{i=0}^{N-2} [100(x_{i+1} - x_i^2)^2 + (1 - x_i)^2]$$

The global minimum is $f(\mathbf{x}) = 0$, which occurs at $\mathbf{x} = (1, 1, \ldots, 1)$.

Nelder-Mead Simplex Algorithm

The Nelder-Mead algorithm is a direct search method that doesn’t require gradient information. It’s often robust for non-smooth functions but can be slower than gradient-based methods.

import numpy as np
from scipy.optimize import minimize

# Define the Rosenbrock function for N variables
def rosenbrock_func(x):
    """The Rosenbrock function."""
    return sum(100.0 * (x[i+1] - x[i]**2.0)**2.0 + (1 - x[i])**2.0
               for i in range(len(x)-1))

# Initial guess (e.g., for a 5-dimensional Rosenbrock function)
x_initial = np.array([1.5, 0.8, 1.2, 0.7, 1.1])

print(f"Initial guess x0: {x_initial}")
print(f"Rosenbrock function value at initial guess: {rosenbrock_func(x_initial):.4f}")

# Use minimize with the Nelder-Mead method
# `method='Nelder-Mead'`
result_nelder_mead = minimize(rosenbrock_func, x_initial, method='Nelder-Mead',
                              options={'xatol': 1e-8, 'disp': True})

print("\nResult from Nelder-Mead:")
print(f"Optimal x: {result_nelder_mead.x}")
print(f"Function value at optimum: {result_nelder_mead.fun}")
print(f"Success: {result_nelder_mead.success}")
# print(f"Message: {result_nelder_mead.message}")

The output will show the optimization process (if disp=True) and the final result:

Initial guess x0: [1.5 0.8 1.2 0.7 1.1]
Rosenbrock function value at initial guess: 145.1000
Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 378
         Function evaluations: 662

Result from Nelder-Mead:
Optimal x: [1.         1.         1.         1.         0.99999999]
Function value at optimum: 5.959159901745691e-17
Success: True

The simplex algorithm is simple and doesn’t need gradients, making it suitable for some problems. However, for smoother functions, gradient-based methods like BFGS (Broyden-Fletcher-Goldfarb-Shanno), available via method='BFGS', are often much faster and more accurate. To use BFGS, you might need to provide the Jacobian (gradient) of the function, or the algorithm can approximate it.

Least-Squares Minimization with least_squares

Least-squares optimization aims to find parameters $\mathbf{x}$ that minimize the sum of squares of $m$ residual functions $r_i(\mathbf{x})$: $Cost(\mathbf{x}) = \sum_{i=0}^{m-1} r_i(\mathbf{x})^2$. This is very common in data fitting (curve fitting).

Let’s again use the Rosenbrock function, but frame it as a least-squares problem. The Rosenbrock function can be written as $f(x_0, x_1) = r_0(x_0, x_1)^2 + r_1(x_0, x_1)^2$, where $r_0 = 10(x_1 - x_0^2)$ and $r_1 = (1 - x_0)$.

import numpy as np
from scipy.optimize import least_squares

# Define the residual functions for Rosenbrock (2D case for simplicity here)
def rosenbrock_residuals(x):
    # x is expected to be an array [x0, x1]
    r0 = 10 * (x[1] - x[0]**2)
    r1 = (1 - x[0])
    return np.array([r0, r1])

# Initial guess for x = [x0, x1]
x_initial_ls = np.array([2.0, 2.0])

print(f"Initial guess for least_squares: {x_initial_ls}")
print(f"Residuals at initial guess: {rosenbrock_residuals(x_initial_ls)}")

# Solve the least-squares problem
result_least_squares = least_squares(rosenbrock_residuals, x_initial_ls)

print("\nResult from least_squares:")
print(f"Optimal x: {result_least_squares.x}")
print(f"Value of residuals at optimum: {result_least_squares.fun}") # These are the r_i values
print(f"Cost function (sum of squares of residuals): {result_least_squares.cost}")
print(f"Success: {result_least_squares.success}")
# print(f"Message: {result_least_squares.message}")

The least_squares function requires a function that returns the vector of residuals. It then minimizes the sum of their squares. The exact minimum for the 2D Rosenbrock function is at $x = [1.0, 1.0]$.

The output shows convergence to the known minimum:

Initial guess for least_squares: [2. 2.]
Residuals at initial guess: [-20.  -1.]

Result from least_squares:
Optimal x: [1. 1.]
Value of residuals at optimum: [4.4408921e-15 0.0000000e+00]
Cost function (sum of squares of residuals): 9.866924291084687e-30
Success: True

For curve fitting specifically, scipy.optimize.curve_fit is a specialized and often more convenient function that uses least-squares to fit a user-defined function to data.

Root Finding

Root finding involves solving equations of the form $f(x) = 0$ (for scalar functions) or $\mathbf{f}(\mathbf{x}) = \mathbf{0}$ (for systems of equations).

Scalar Functions Root Finding

For single-variable equations, scipy.optimize offers several algorithms like brentq (often the best choice if an interval bracketing the root is known), newton (requires derivative or approximates it), and bisect.

Systems of Equations with root()

To find roots of a set of non-linear equations $\mathbf{f}(\mathbf{x}) = \mathbf{0}$, use the root() function. It supports various methods, including ‘hybr’ (Powell’s hybrid method, default) and ‘lm’ (Levenberg-Marquardt).

Consider solving the single-variable transcendental equation $2x + 2\cos(x) = 0$. We are looking for $x$ such that $f(x) = 2x + 2\cos(x)$ is zero.

import numpy as np
from scipy.optimize import root

# Define the function whose root we want to find
def func_for_root(x):
    # Note: x here is an array, even for a single variable problem with root()
    return 2 * x[0] + 2 * np.cos(x[0])

# Initial guess for the root
initial_guess_root = np.array([0.0]) # Must be an array for root()

solution = root(func_for_root, initial_guess_root)

print("Solution from root():")
print(f"Root x: {solution.x}")
print(f"Function value at root f(x): {solution.fun}") # Should be close to 0
print(f"Success: {solution.success}")
# print(f"Message: {solution.message}")

The output will show the found root:

Solution from root():
Root x: [-0.73908513]
Function value at root f(x): [0.]
Success: True

This means $x approx -0.739$ is a root of the equation $2x + 2\cos(x) = 0$.

Practical Applications of Optimization

Optimization techniques are fundamental in:

• Machine Learning: Training models by minimizing loss functions.
• Engineering: Designing structures, processes, or systems for maximum efficiency or minimum cost.
• Finance: Portfolio optimization, risk management.
• Logistics: Route planning, resource allocation.
• Scientific Research: Parameter estimation in models to fit experimental data.

Choosing the right optimization algorithm depends on the problem’s characteristics: whether it’s constrained or unconstrained, if derivatives are available, the function’s smoothness, and whether a global or local minimum is sought. Always consult the SciPy documentation for detailed guidance on each optimizer’s strengths and weaknesses.