Table of Contents

Topic 27. Gradient Descent¶

Author: Jordan Brown (3jordanbrown@gmail.com)

Introduction¶

Gradient descent, or steepest descent, is a function optimization method which iteratively approaches a minimum-cost solution for a function. It accomplishes this goal by updating the function after each iteration in such a way that it moves in the "direction of motion," or towards a minimum (or at least a local minimum).

Gradient descent is important as it allows a simple, fairly resource non-intensive way to solve difficult problems for which a precise model may be unkown, or direct minima calculation is complex. It also allows for adaptive solutions to such problems where new information may become available, and the model needs to be updated to reflect the new information. This is an idea commonly used in machine learning, for example, to train a neural network model to best approximate some unknown model for which inputs and desired outputs are known (like in the Engineering Application section below).

Explanation of the theory¶

In my opinion, one of the easiest ways to introduce gradient descent is by visualizing a simplified example. For example, consider some equation for which the performance was evaluated by observing the error between the estimated model output, and some real desired output.

Now imagine that we didn't know the shape of the error curve. To move towards a minimum in the error, what we could do is differentiate the error function with regards to some parameter. Then the slope will tell us how to adjust that parameter so that the error function gets closer to its minimum. The adjustment is done by taking a "step" towards the minimum, or changing x by the derivative evaluated at the current parameter value, multiplied by some step size that we will call $\mu$. This is represented as

$$\begin{equation} x_{k+1}=x_{k}-\mu \frac{\delta E}{\delta x}, \end{equation}$$ where $x_{k+1}$ is the updated value, and $x_{k}$ is the previous. Let's now say we have a fairly large step size, and our updated parameter swings us to the other side of the minimum as in the next plot.

By looking at the parameter update equation, you can see that if the slope swings positive, the adjustment will now bring us backwards towards the minimum. We can continue this method, approaching the minimum closer and closer until we actually find a derivative of zero, get under some desired threshold of error, or reach some maximum number of iterations.

Now that we have this visualized example, we can dig into the details. To define the above idea in a mathematical way, the goal of steepest descent is to keep updating a function such that

$$\begin{equation} f(x_{k+1}) < f(x_{k}). \end{equation}$$

We know we have reached a minimum (or other extremum) when

$$\begin{equation} f(x_{k+1}) = f(x_{k}), \end{equation}$$

since the derivative of a function at a minimum or maximum is zero and our above-defined function-updating equation doesn't change if the slope is zero.

Overall, steepest descent is a simple idea which yields good results. Almost all of the complexity of this method lies simply in selecting the right value of $\mu$ for our steps. In our original example, with an $x^2$ shape, if we would have selected a step value of 0.5, we would have oscillated back and forth and never converged to the minimum. Values larger than 0.5 would even cause us to diverge away from the minimum. Thus we need to be wise in our selection of $\mu$ to ensure convergence to a minimum solution. There are some derived requirements for $\mu$ in example 14.5 of [1]. Beginning with an example function defined as

$$\begin{equation} f(x) = x^TRx-2b^Tx, \end{equation}$$

where x $\in R^n$ and $R$ is positive semi-definite. For this example the resultant step requirement was

$$\begin{equation} 0 < \mu < \frac{2}{\lambda_{max}}, \end{equation}$$

where $\lambda_{max}$ is the largest eigenvalue of R. Similar rules can be found mathematically for equations where practical, or with a trial-and-error approach for more complex solutions.

Simple Numerical Examples¶

We will work through a Least-Mean Squares adaptive filter as an example to emphasize the concepts presented above. The discrete filter can be defined as

$$\begin{equation} \mathbf{y}[t]=\mathbf{f}[t]^T\mathbf{h}, \end{equation}$$

where $\mathbf{f}$ is the input and $\mathbf{h}$ is the filter impulse response, and our cost function $J$ which we will seek to minimize is defined as

$$\begin{equation} J(\mathbf{h})=E\{(\mathbf{d}[t]-\mathbf{y}[t])^2\}. \end{equation}$$

Following the derivation in [1], equation 14.13 gives us an update rule (based on the derivative) which is

$$\begin{equation} \mathbf{h}[t+1]=\mathbf{h}[t]+\mu(E\{\mathbf{f}[t]\mathbf{d}[t]\}-E\{\mathbf{f}[t]\mathbf{f}[t]^T\}\mathbf{h}[t]). \end{equation}$$

Say our true system results in a portion of a cubic function, and we start with our adaptive filter impulse response as $h=\mathbf{0}$. From there, follow the plots and code as we progressively update $h$, and observe the inputs passed through the our system which is adapting its impulse response so that the outputs more closely match the desired response. We'll let the input be time. Feel free to adjust $\mu$ and observe how it affects the fit. (Note: in these plots, the green line is our adaptive filter output, and the blue line is our desired filter output)

An Engineering Application¶

Scatterometry is a field where microwave radar scattering off of the earth's surface is used to improve natural hazard and weather prediction, and study global climate patterns. Many, if not all of the relationships between physical phenomena and radar return echo are too complex, however, to be derived from known equations, especially in the presence of noise. This means methods like gradient descent can be very useful in deriving heuristic models based on real data. For example, the Advanced Scatterometer (ASCAT) was designed to measure ocean winds. It observes radar backscatter ($\sigma$) as seen from 3 different azimuth look angles. Current models, known as Geophysical Model Functions (GMF), are used to retrieve wind speed and direction from these backscatter values. GMFs, however, have inherent ambiguity in their solutions. Thus there is room to improve our models to better discriminate between similar solutions. Here we will explore the idea of creating a model using gradient descent to minimize our error and improve our transfer function. Although this won't be a great model, it will illustrate the idea.

To make things a simple, we'll create our model in this format:

$$\begin{equation} WSPD_{mod}=a\sigma_{av}+b\sigma_{av}^2+c\sigma_{av}^3, \end{equation}$$

where $\sigma_{av}$ is the average of the backscatter from the 3 azimuth look angles (fore, mid, and aft)

$$\begin{equation} \sigma_{av} = (\sigma_{f} + \sigma_{m} + \sigma_{a})/3. \end{equation}$$

If we define the error as

$$\begin{equation} err=(WSPD-WSPD_{mod})^2, \end{equation}$$

then taking the derivative of the error with regards to each of the three adjustable parameters ($a$, $b$, $c$) leads to our update steps, which can be seen in the code below.

The data is actual data from a region near Hawaii. All speeds are reported in meters/second, and all backscatter measurements are in decibels.

References¶

[1] Moon, Todd K., and Wynn C. Stirling. Mathematical Methods and Algorithms for Signal Processing. Prentice Hall, 2000.

[2] R. Beard, T. McLain, C. Peterson. Introduction to Feedback Control: Using Design Studies. 2020.

Homework Problem:¶

Author: Spencer Jensen

Gradient descent can be used for control of systems. For any system, a quadratic cost function on states and inputs can be respresented as

$$\begin{equation} J = (x-x_{desired})^{T}Q(x-x_{desired})+u^{T}Ru \end{equation}$$

where $J$ is the scalar cost, $x$ is a nx1 vector of weights, $u$ is a mx1 vector of inputs and $Q$ and $R$ are nxn and mxm matrices of the weights on states and inputs respectively. Control of this system can be stated as a gradient descent problem where

$$\begin{equation} u_{k+1}=u_{k}-\alpha\frac{\partial J}{\partial u}(u_{k}) \end{equation}$$

Here $u_{k+1}$ represents the input at the next timestep i.e. you will use a single step of gradient descent to find the next u to apply. One can estimate the gradient of $J$ numerically by

$$\begin{equation} \frac{\partial J}{\partial u}(u_{k}) = \frac{J(u_{k}+h)-J(u_{k})}{h} \end{equation}$$

In this problem, you will simulate the control of the mass-spring-damper system shown below.

Given an initial $x_{o}$, $x_{desired}$, $h$, $Q$, $R$, $\alpha$ and $u_{o}$ and a function $x_{k+1} = f(x_{k},u_{k})$ written in "mass_spring_damper.py", use gradient descent to find the inputs and states for 50 seconds (dynamics function returns the next state at .01 seconds). Show graphs for position, velocity and input force.

$$\begin{equation} \begin{aligned} u_{o} = 0 \qquad x_{o} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \qquad x_{desired} = \begin{pmatrix} 5 \\ 0 \end{pmatrix} \qquad h = .01 \qquad Q = \begin{pmatrix} 2000 & 0 \\ 0 & 100 \end{pmatrix} \qquad R = \begin{pmatrix} .00001 \end{pmatrix} \qquad \alpha = 5 \end{aligned} \end{equation}$$

Hints: You will probably want something like:

    from mass_spring_damper import Mass_spring_damper
    sys = Mass_spring_damper()
    newx = sys.get_next_xk(x,u)

    def J(x,x_des,u):
        Q = np.array([[2000,0],[0,100]])
        R = .00001
        cost = (x-x_des).T@Q@(x-x_des)+u*R*u
        return cost