Optimize Like a Pro: Numerical Methods in SciPy#

The Golden-Section Search and Brent’s Method#

A common technique for single-variable functions is Golden-Section Search, which doesn’t require derivatives. The interval is repeatedly shrunk in a way inspired by the golden ratio until convergence.

Brent’s Method is another widely used solver that combines parabolic interpolation (faster convergence if the function is smooth) with the robust bracketing of Golden-Section Search. This hybrid approach often converges quickly.

Practical Usage with SciPy’s minimize_scalar#

Here’s a quick example of how to use minimize_scalar:

1
import numpy as np
2
from scipy.optimize import minimize_scalar
3

4
# Define a simple objective function
5
def objective(x):
6
    return (x - 2)**2 + 3
7

8
# Using the default Brent method
9
result = minimize_scalar(objective)
10
print("Minimum found at:", result.x)
11
print("Function value:", result.fun)

Output might look like:

1
Minimum found at: 2.0
2
Function value: 3.0

Unconstrained Multivariate Minimization Algorithm Overview#

1. Nelder-Mead: A gradient-free algorithm that uses a simplex and is quite robust but can be slow in higher-dimensional spaces.
2. BFGS: A quasi-Newton method that approximates the Hessian matrix. Often performs well if gradients are available or can be approximated.
3. CG (Conjugate Gradient): Suitable when your problem has many variables, but the function is well-behaved (e.g., no heavy constraints).
4. Newton-CG: Uses the Hessian or Hessian-vector products for faster convergence if the Hessian is reasonably easy to compute.

Gradient-Based Methods vs. Gradient-Free Methods#

• Gradient-Based:

  • Pros: Can converge quickly if the function is smooth and derivatives are easily computed.
  • Cons: Might not handle non-smooth or noisy functions well.

• Gradient-Free:

  • Pros: More robust for noisy or non-smooth functions.
  • Cons: Might be slower in high-dimensional problems because they lack direction from gradient information.

Examples With scipy.optimize.minimize#

Here’s an example using BFGS:

1
import numpy as np
2
from scipy.optimize import minimize
3

4
def rosenbrock(x):
5
    # Rosenbrock function, often used as a performance test problem
6
    return 100.0*(x[1] - x[0]**2)**2 + (1 - x[0])**2
7

8
x0 = np.array([-1.2, 1.0])  # Initial guess
9
res = minimize(rosenbrock, x0, method='BFGS')
10
print(res)

Output:

1
      fun: 1.7274e-12
2
 hess_inv: array([[ 0.9589783,  1.9162314],
3
                  [ 1.9162314,  3.8791316]])
4
      jac: array([  1.256e-06,  -5.332e-07])
5
  message: 'Optimization terminated successfully.'
6
     nfev: 65
7
      nit: 50
8
     njev: 65
9
   status: 0
10
  success: True
11
        x: array([1.0000034 , 1.0000068])

This shows the function value, approximate Hessian inverse, gradient, and final solution. The Rosenbrock function’s known global minimum is at [1,1].

Linear Constraints#

For linear constraints of the form:

• A·x �?b
• A_eq·x = b_eq

SciPy’s minimize function allows specifying dictionaries for constraints:

1
import numpy as np
2
from scipy.optimize import minimize
3

4
def objective(x):
5
    return x[0]**2 + x[1]**2
6

7
# Define constraints: x0 + x1 <= 1 (inequality)
8
constraints = ({
9
    'type': 'ineq',
10
    'fun': lambda x: 1 - (x[0] + x[1])
11
})
12

13
# Initial guess
14
x0 = np.array([0.5, 0.5])
15

16
res = minimize(objective, x0, method='SLSQP', constraints=constraints)
17
print(res)

The SLSQP (Sequential Least Squares Programming) method handles both equality and inequality constraints smoothly.

Nonlinear Constraints#

Nonlinear constraints are defined similarly, but the constraint function can be any callable returning a scalar or vector. For instance:

1
constraints = ({
2
    'type': 'eq',
3
    'fun': lambda x: np.array([x[0]**2 + x[1] - 4]),
4
})

Here, we want x₀² + x�?- 4 = 0.

Practical Example: Portfolio Optimization#

A practical use case is optimizing a portfolio for minimal risk subject to a target return:

• Objective function: a risk measure (e.g., variance).
• Constraints: sum of weights = 1 (equality), and each weight �?0 (inequality).

A minimal working example:

1
import numpy as np
2
from scipy.optimize import minimize
3

4
def portfolio_variance(weights, cov_matrix):
5
    return weights.T @ cov_matrix @ weights
6

7
cov_matrix = np.array([
8
    [0.1, 0.02, 0.01],
9
    [0.02, 0.08, 0.03],
10
    [0.01, 0.03, 0.15]
11
])
12

13
# Constraints: sum(weights) = 1, each weight >= 0
14
cons = (
15
    {'type': 'eq', 'fun': lambda w: np.sum(w) - 1}
16
)
17

18
bounds = [(0,1), (0,1), (0,1)]
19
initial_guess = [1/3, 1/3, 1/3]
20

21
res = minimize(portfolio_variance, initial_guess,
22
               args=(cov_matrix,),
23
               method='SLSQP',
24
               constraints=cons,
25
               bounds=bounds)
26

27
print(res)

Basics of Global vs. Local Optimization#

Local optimization methods converge to the nearest local optimum. In contrast, global optimization algorithms systematically search over the entire domain to find the absolute minimum (or maximum). Global algorithms are usually more computationally expensive.

basinhopping and differential_evolution#

Basin-hopping:

• Performs local optimization repeatedly but perturbs the solution to escape local minima.

Differential Evolution:

• A population-based stochastic algorithm that evolves solutions over generations.
• Belongs to the family of evolutionary algorithms.

1
import numpy as np
2
from scipy.optimize import basinhopping
3

4
def f(x):
5
    return np.sin(5*x[0]) + (x[0] - 2)**2
6

7
x0 = [0.0]
8
minimizer_kwargs = {"method": "BFGS"}
9
res = basinhopping(f, x0, minimizer_kwargs=minimizer_kwargs, niter=100)
10
print(res)

Case Study: Parameter Estimation#

Imagine you’re fitting a model parameter θ to data. The model might be complex, non-convex, and noisy. A global optimization method like differential evolution can be used for robust parameter estimation:

1
import numpy as np
2
from scipy.optimize import differential_evolution
3

4
def model_error(params):
5
    # Hypothetical residual sum of squares
6
    return np.sum((observed_data - complicated_model(params))**2)
7

8
bounds = [(0,10), (0,5), (1,100)]
9
result = differential_evolution(model_error, bounds)
10
print(result)

This approach is especially helpful when the solution surface is highly non-convex with multiple local minima.

Why Root Finding is Key to Optimization#

Root finding comes into play when your optimization problem can be turned into finding a zero of the gradient or derivative. SciPy’s optimize suite provides powerful root-finding algorithms for single and multivariate equations.

Popular Methods: Bisection, Newton, and Secant#

1. Bisection requires brackets [a, b] where f(a) and f(b) have opposite signs.
2. Newton uses derivative information for fast convergence but may diverge if the initial guess is poor.
3. Secant approximates the derivative from previous function evaluations.

Practical Example With SciPy Functions#

1
from scipy.optimize import bisect, newton, root_scalar
2

3
def f(x):
4
    return x**2 - 2
5

6
# Bisection
7
root_bisect = bisect(f, 0, 2)
8
print("Bisection root:", root_bisect)
9

10
# Newton
11
root_newton = newton(f, 1.0)
12
print("Newton root:", root_newton)
13

14
# root_scalar can handle multiple methods
15
res = root_scalar(f, bracket=[0, 2], method='brentq')
16
print("BrentQ root:", res.root)

All three methods should yield �? �?1.414213562.

Custom Objective Functions and Gradients#

When the default symbolic gradient is not sufficient or you have a complex function, you can provide a custom gradient (jacobian) and Hessian:

1
import numpy as np
2
from scipy.optimize import minimize
3

4
def objective_and_grad(x):
5
    # f = (x0^2 + x1 - 11)^2 + (x0 + x1^2 -7)^2
6
    # A well-known two-variable function for testing
7
    f = (x[0]**2 + x[1] - 11)**2 + (x[0] + x[1]**2 - 7)**2
8
    # Analytical gradients
9
    dfdx0 = 2*(x[0]**2 + x[1] - 11)*(2*x[0]) + 2*(x[0] + x[1]**2 - 7)*1
10
    dfdx1 = 2*(x[0]**2 + x[1] - 11)*1 + 2*(x[0] + x[1]**2 - 7)*(2*x[1])
11
    return f, np.array([dfdx0, dfdx1])
12

13
def objective(x):
14
    return objective_and_grad(x)[0]
15

16
def grad(x):
17
    return objective_and_grad(x)[1]
18

19
x0 = np.array([0, 0])
20
res = minimize(objective, x0, jac=grad, method='BFGS')
21
print(res)

Automatic Differentiation Integration#

If you’re using libraries like JAX, PyTorch, or TensorFlow, you can leverage automatic differentiation to compute gradients. Then, you can feed those gradients into SciPy’s routines, combining the best of both worlds.

Handling Large-Scale Nonlinear Problems#

For problems with thousands or millions of parameters, you might need to:

• Use specialized large-scale solvers (L-BFGS, truncated Newton).
• Exploit sparse matrix structures to handle large Hessians or Jacobians.
• Parallelize or distribute computations across multiple processors.

Sparse Matrices and Solver Tricks#

SciPy offers sparse matrices (scipy.sparse) that store only non-zero elements. Coupled with specialized linear algebra routines, you can handle huge systems more efficiently. For instance, if your Hessian is sparse, you can use the sparse.linalg subpackage to factor or solve large linear systems quickly.

Vectorization Strategies#

Where possible, vectorize your functions using NumPy arrays for faster evaluations. For instance:

1
import numpy as np
2

3
def vectorized_objective(x):
4
    # Suppose x is an array of parameters, we compute something in a single pass
5
    return np.sum(x**2)

By avoiding Python loops, the overhead is greatly reduced.

Caching and Memoization#

For expensive function evaluations, caching can be a lifesaver. Use Python’s functools.lru_cache (with care over mutable inputs) or implement your own caching logic. This is especially useful if your optimization algorithm repeatedly evaluates the same points (as in certain derivative-free methods).

1
from functools import lru_cache
2

3
@lru_cache(None)
4
def expensive_calc(x):
5
    # Convert x to a tuple if it’s not hashable
6
    # Perform an expensive calculation
7
    pass

Parallel Computing Considerations#

Some optimization approaches support parallel objective evaluations (like evolutionary algorithms). SciPy’s built-in routines may not fully exploit parallelization, but you can:

• Use the parallel argument in certain functions (e.g., differential_evolution in SciPy can accept a ‘workers�?parameter).
• Roll out your own logic with Python’s multiprocessing or joblib to distribute computations of the objective function.