Unconstrained Optimization

Teng-Jui Lin

Content adapted from UW AMATH 301, Beginning Scientific Computing, in Spring 2020.

• Unconstrained optimization methods

  • Derivative-free method

    • Golden section search

  • Derivative method

    • Gradient descent

• scipy implementations

  • Minimize scalar functions by scipy.optimize.fminbound()

  • Minimize multivariable functions by scipy.optimize.fmin()

Golden section search

Goal: find the minimum of a single-variable unimodal function $f(x)$.

• unconstrained

• derivative-free

In a given interval $[a,b]$, we find two points $x_1$ and $x_2$ where $a<x_1<x_2<b$. The two points can be compared to determine the new interval:

• If $f(x_1)\le f(x_2)$, then change the new interval to $[a,x_2]$

• If $f(x_1)>f(x_2)$, then change the new interval to $[x_1,b]$

The golden section search minimizes the number of calculation rather than iteration, so we want to

• have a constant reduction factor $c$

  • If $f(x_1)\le f(x_2)$, $c={\displaystyle \frac{x_2-a}{b-a}}⟹x_2=(1-c)a+cb$

  • If $f(x_1)>f(x_2)$, $c={\displaystyle \frac{b-x_1}{b-a}}⟹x_1=ca+(1-c)b$

• reuse calculated points

  • $x_2^{\text{new}}=x_1^{\text{old}}=ca+(1-c)b$

  • $x_2^{\text{new}}=(1-c)a+c{x}_2^{\text{old}}=(1-c)a+c[(1-c)a+cb]$

Equate the above two equations, we get

$$\begin{array}{r}{c}^2+c+1=0\\ c={\displaystyle \frac{\sqrt{5}-1}{2}}\approx 0.618\end{array}$$

Therefore, we can reuse the points during iteration by choosing:

• If $f(x_1)\le f(x_2)$

  • $x_1^{\text{new}}=x_2^{\text{old}}$

  • $x_2^{\text{new}}=(1-c)a+cb$

• If $f(x_1)>f(x_2)$

  • $x_2^{\text{new}}=x_1^{\text{old}}$

  • $x_1^{\text{new}}=ca+(1-c)b$

Implementation

Problem Statement. Find the minimum of the function

$$f(x)=-{\displaystyle \frac{1}{3}}{e}^{-x/2}-3x{e}^{-x/2}$$

using golden section search.

import numpy as np

def golden_section_search(f, a, b, tolerance=1e-6, max_iter=5000):
    '''
    Find the minimum of single-variable unimodal function using golden section search method.
    
    :param f: objective function
    :param a: lower bound of interval
    :param b: upper bound of interval
    :param tolerance: tolerance of stopping criteria
    :param max_iter: maximum iterations allowed for calculation
    :returns: minimum of objective function
    '''
    c = (np.sqrt(5) - 1)/2
    x1 = c*a + (1 - c)*b
    x2 = (1 - c)*a + c*b
    f1 = f(x1)
    f2 = f(x2)
    
    i = 0  # iteration counter
    
    # golden section search logic
    while (b-a) > tolerance and i < max_iter:
        if f1 <= f2:
            b = x2
            x2 = x1
            f2 = f1
            x1 = c*a + (1 - c)*b
            f1 = f(x1)
        else:
            a = x1
            x1 = x2
            f1 = f2
            x2 = (1 - c)*a + c*b
            f2 = f(x2)
    
    # maximum iteration warning
    if i == max_iter:
        import warnings
        warnings.warn(f'Maximum iteration reached. Current stopping criteria is {b-a}', UserWarning)
    
    result = (a+b)/2  # take avg of upper and lower bound
    
    return result

def test_func(x):
    return -(1/3*np.exp(-x/2) + 3*x*np.exp(-x/2))

a = 0
b = 10
golden_section_search(test_func, a, b)

1.888888807141935

Finding minimum using scipy.optimize.fminbound()

scipy.optimize.fminbound() find the minimum for scalar function in given interval.

import scipy
from scipy import optimize

scipy.optimize.fminbound(test_func, a, b)

1.8888887480521486

Gradient descent

Goal: find the minimum of a multi-variable function $f(\mathbf{x})$

• unconstrained

• derivative method

Gradient descent chooses the next point in the steepest downhill direction, where the function value is minimum. For each guess ${\mathbf{p}}_{\mathbf{n}}$, we calculate ${\mathbf{p}}_{\mathbf{n}\mathbf{+}\mathbf{1}}$ by

1. Calculate $\mathrm{\nabla }f({\mathbf{p}}_{\mathbf{n}})$

2. Define $\varphi (t)={\mathbf{p}}_{\mathbf{n}}-t\mathrm{\nabla }f({\mathbf{p}}_{\mathbf{n}})$

3. Define $f(\varphi (t))$

4. Find $t^{\ast }$ using fminbound() to minimize $f(\varphi (t))$

5. Define ${\mathbf{p}}_{\mathbf{n}\mathbf{+}\mathbf{1}}=\varphi (t^{\ast })$

Optimization is the speed limiting step.

Implementation

Problem Statement. Find the minimum of the function

$$f(x,y)=(x-2{)}^2+(y+1{)}^2+5\sin x\sin y+100$$

using gradient descent and compare with the result from scipy.optimize.fmin().

def gradient_descent(f, fgrad, x0, tolerance=1e-6, max_iter=10000):
    '''
    Find a minimum of multivariable function by gradient descent.
    
    :param f: objective function
    :param fgrad: gradient of objective function
    :param x0: initial guess (starting point)
    :param tolerance: tolerance of stopping criteria
    :param max_iter: maximum iterations allowed for calculation
    :returns: minimum of objective function
    '''
    import scipy
    from scipy import optimize
    
    i = 0  # iteration counter
    x = x0
    grad = fgrad(x)
    
    # gradient descent logic
    while np.linalg.norm(grad, np.inf) > tolerance and i < max_iter:
        grad = fgrad(x)
        phi = lambda t : x - t*grad
        f_of_phi = lambda t : f(phi(t))
        tmin = scipy.optimize.fminbound(f_of_phi, 0, 1)
        x = phi(tmin)  # use next point with minimum function value
        i += 1
    
    # maximum iteration warning
    if i == max_iter:
        import warnings
        warnings.warn(f'Maximum iteration reached. Current stopping criteria is {np.linalg.norm(grad, np.inf)}', UserWarning)
    
    return x

fxy = lambda x, y : (x - 2)**2 + (y + 1)**2 + 5*np.sin(x)*np.sin(y) + 100
f = lambda p : fxy(*p)

fgradx = lambda x, y : 2*(x - 2) + 5*np.cos(x)*np.sin(y) 
fgrady = lambda x, y : 2*(y + 1) + 5*np.sin(x)*np.cos(y)
fgrad =  lambda p : np.array([fgradx(*p), fgrady(*p)])

x0 = [6, 4]
gradient_descent(f, fgrad, x0)

array([ 1.69484631, -1.40628341])

# compare with scipy
scipy.optimize.fmin(f, x0)

Optimization terminated successfully.
         Current function value: 95.363597
         Iterations: 42
         Function evaluations: 82

array([ 1.69487388, -1.40630558])

Gradient descent with fixed steps

Goal: find the minimum of a multi-variable function $f(\mathbf{x})$

• unconstrained

• derivative method

Gradient descent chooses the next point in the steepest downhill direction, with a constant learning rate (fixed step) $t$. For each guess ${\mathbf{p}}_{\mathbf{n}}$, we calculate ${\mathbf{p}}_{\mathbf{n}\mathbf{+}\mathbf{1}}$ by

1. Calculate $\mathrm{\nabla }f({\mathbf{p}}_{\mathbf{n}})$

2. Define ${\mathbf{p}}_{\mathbf{n}\mathbf{+}\mathbf{1}}={\mathbf{p}}_{\mathbf{n}}-t\mathrm{\nabla }f({\mathbf{p}}_{\mathbf{n}})$

Calculating gradient is the speed limiting step. The learning rate need to be tuned.

Implementation

Problem Statement. Find the minimum of the function

$$f(x,y)=(x-2{)}^2+(y+1{)}^2+5\sin x\sin y+100$$

using gradient descent with fixed steps and compare with the result from scipy.optimize.fmin().

def gradient_descent_fixed_step(f, fgrad, x0, learning_rate, tolerance=1e-6, max_iter=10000):
    '''
    Find a minimum of multivariable function by gradient descent.
    
    :param f: objective function
    :param fgrad: gradient of objective function
    :param x0: initial guess (starting point)
    :param learning_rate: learning rate
    :param tolerance: tolerance of stopping criteria
    :param max_iter: maximum iterations allowed for calculation
    :returns: minimum of objective function
    '''
    import scipy
    from scipy import optimize
    
    i = 0  # iteration counter
    x = x0
    grad = fgrad(x)
    
    # gradient descent with fixed step logic
    while np.linalg.norm(grad, np.inf) > tolerance and i < max_iter:
        grad = fgrad(x)
        x = x - learning_rate*grad  # use next point with fixed step
        i += 1
    
    # maximum iteration warning
    if i == max_iter:
        import warnings
        warnings.warn(f'Maximum iteration reached. Current stopping criteria is {np.linalg.norm(grad, np.inf)}', UserWarning)
    
    return x

x0 = [6, 4]
gradient_descent_fixed_step(f, fgrad, x0, 0.2)

array([ 1.69484629, -1.40628343])

# compare with scipy
scipy.optimize.fmin(f, x0)

Optimization terminated successfully.
         Current function value: 95.363597
         Iterations: 42
         Function evaluations: 82

array([ 1.69487388, -1.40630558])

Finding minimum using scipy.optimize.fmin()

scipy.optimize.fmin() finds a minimum for multivariable function.

scipy.optimize.fmin(f, x0)

Optimization terminated successfully.
         Current function value: 95.363597
         Iterations: 42
         Function evaluations: 82

array([ 1.69487388, -1.40630558])