<jupyter_start><jupyter_text>![](img/572_banner.png)  Appendix C: Computing Derivatives<jupyter_code>import sklearn
import sklearn.linear_model
import scipy.optimize
import autograd # pip install autograd
from autograd import grad
import autograd.numpy as anp
from scipy.optimize import minimize
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams.update({'font.size': 16, 'axes.labelweight': 'bold', 'figure.figsize': (8,6)})<jupyter_output><empty_output><jupyter_text>1. Gradients Introduction What if we are optimizing a function and we don't know its derivative/gradient. This could happen because:- You're optimizing the parameters of a computer simulation  - a robot  - weather / atmospheric science  - financial simulation  - etc.- You're optimizing the hyperparameters of a machine learning algorithm- The derivative is just too much effort to compute and it's slowing down your prototyping effortsWhat to do? Well, there are 3 main approaches to computing derivatives with a computer:1. Symbolic differentiation2. Numerical differentiation3. Automatic differentiationWe'll explore these different options in this short appendix. Note that this material has been modified after material originally created by [<NAME>](https://www.mikegelbart.com/).As a running example, consider $f(x)=\sin(x)$. We know that $f'(x)\equiv\frac{df(x)}{dx}=\cos(x)$ and thus $f'(1)=\cos(1)$  2. Symbolic Differentiation Symbolic differentiation is basically what you learned in calculus. We interpret the symbols in an equation and apply the rules of differentiation. Wolfram Alpha is a good tool for this, e.g., .Symbolic differentiation is useful but less so when dealing with data. We will not say more about this now. You should be aware of its existence.  3. Numerical Differentiation (Finite Differences) A derivative is defined as:$$\frac{df(x)}{dx} \equiv \lim_{h\rightarrow 0} \frac{f(x+h)-f(x)}{h} $$Numerical differentiation simply approximates the above using a very small $h$:$$\frac{df(x)}{dx} \approx \frac{f(x+h)-f(x)}{h}$$ for small $h$.This approach is called "finite differences". I like to think of it as an estimate of as the "rise over run" estimate of slope. Let's give it a try!<jupyter_code>h = 0.001  # small h
f = np.sin
x0 = 1
dfdx = (f(x0 + h) - f(x0)) / h
dfdx
np.cos(x0)<jupyter_output><empty_output><jupyter_text>Not bad!  3.1. Error of Derivative Estimates In [Lecture 1](lecture1_floating-point.ipynb) we discussed numerical issues due to the **floating point** representation. There's another entirely different type of numerical error that arises from the *discretization* inherent finite difference approximations. Let's explore this a bit. What do you think will happen to the absolute error of the estimate of the derivative if I change $h$?<jupyter_code>f = lambda x: np.sin(x)
dfdx_estimate = lambda h: (f(x0 + h) - f(x0)) / h
dfdx_true = np.cos(x0)
abs_error = lambda h: np.abs(dfdx_estimate(h) - dfdx_true)

h = 10 ** (np.linspace(-8, 0))
plt.loglog(h, abs_error(h))
plt.xlabel("h")
plt.ylabel("Error of derivative estimate");<jupyter_output><empty_output><jupyter_text>The above plot says that the error between our estimate and the true derivative grows as $h$ gets bigger. Makes sense, right? But the plot thickens when we look at very small values of $h$...<jupyter_code>h = 10 ** (np.linspace(-16, 0))
plt.loglog(h, abs_error(h))
plt.xlabel("h")
plt.ylabel("Error of derivative estimate");<jupyter_output><empty_output><jupyter_text>Wait... what?When $h\ll 1$ we have that $\left| f(x+h)-f(x) \right| \ll 1$. This leads to _roundoff errors_. Is it underflow?<jupyter_code>h = 10 ** (np.linspace(-25, 0))
plt.loglog(h, abs_error(h));<jupyter_output><empty_output><jupyter_text>Well that flat part is underflow (estimate of the derivative is exactly zero). But what's up with the part in between, from $h=10^{-16}$ to $10^{-7}$?The problem is regular old roundoff error (go back to [Lecture 1](lecture1_floating-point.ipynb) if you've forgottn what this is). Roughly speaking, the magnitude of these roundoff errors is independent of $h$ (e.g. when computing $f(x)$ it has nothing to do with $h$), and so the $h$ in the denominator causes $\frac{1}{h}$-like behavior. We can test this with a line, using our friend $10^{-16}$ from [Lecture 1](lecture1_floating-point.ipynb).<jupyter_code>h = 10 ** (np.linspace(-16, 0))
plt.loglog(h, abs_error(h))
plt.loglog(h, 1e-16 * dfdx_true / h);<jupyter_output><empty_output><jupyter_text>We've seen a bunch of plots, but what have we actually learned? Well:- the _numerical approximation_ to the derivative initially increases as $h$ gets smaller (this effect is linear in $h$)- however, when $h$ is too small then roundoff error kicks in and the error starts to increase (this effect is linear in $1/h$)What you really need to know: picking $h$ is tricky: it can't be too big or too small.  3.2. Gradients So, what about gradients?$$\nabla f(x) = \begin{bmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \\ \vdots \\ \frac{\partial f}{\partial x_d} \end{bmatrix}$$ Well, we can approximate each partial derivative as:$$\frac{\partial f(x)}{\partial x_1} \approx \frac{f(x+h e_1)-f(x)}{h}$$ where $e_1 = \begin{bmatrix}1 \\ 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}$. In general, $$\frac{\partial f(x)}{\partial x_j} \approx \frac{f(x+h e_j)-f(x)}{h}$$ where $e_j$ is $1$ at position $j$ and zero elsewhere.Put another way (if this helps),$$\frac{\partial f(x)}{\partial x_1} \approx \frac{f\left(\begin{bmatrix}x_1 + h \\ x_2 \\ x_3 \\ \vdots \\ x_d \end{bmatrix} \right)-f(x)}{h}$$ $$\frac{\partial f(x)}{\partial x_2} \approx \frac{f\left(\begin{bmatrix}x_1 \\ x_2 + h \\ x_3 \\ \vdots \\ x_d \end{bmatrix} \right)-f(x)}{h}$$ And so on and so forth.  3.3. Cost of Estimating a Gradient As you can see above, estimating a gradient requires $d+1$ function evaluations. That can be expensive - perhaps $d$ is $1000$ or $10000$. If you can compute the gradient symbolically, that is much better.  3.4. (Optional) Other Finite Difference Formulas There are better finite differences formulas than the one we are using, such as the centered difference formula:$$\frac{df(x)}{dx} \approx \frac{f(x+h)-f(x-h)}{2h}$$In this case the error goes down as $h^2$, which is great, but most of the benefits come in applications that don't concern us here (e.g., solving differential equations). However, in $d$ dimensions, this requires $2d$ evaluations of the function instead of $d+1$, because with the forward difference formula (the original one above) we reuse the point $f(x)$ for all the partial derivatives. Really, which formula you use depends on what you're doing.If you're really bored, you can also check out the [complex step method](https://timvieira.github.io/blog/post/2014/08/07/complex-step-derivative/).  3.5. Summary If we can't figure out the derivative of a function, but have code that computes the function, we can _approximate_ the derivative using finite differences. This works pretty well but there are some pitfalls. In particular you need to choose $h$ carefully and it costs $O(d)$ function evaluations to estimate a gradient. It seems better to actually know the derivative if possible, because the computations will be faster and more accurate.  4. Automatic Differentiation Automatic differentiation (AD) is another way to take the derivative given some code that computes $f$. It is advertised as the "best of both worlds" between symbolic and numerical differentiation:  - Pro: takes the code as input, not math.  - Pro: does not have the problem of choosing $h$, or the extra $O(d)$ cost.  - Con: some AD software packages have limitations, hurt code readability.  - Con: large space complexity.  There are two main types of AD, _forward mode_ and _reverse mode_ AD. We'll be discussing reverse mode here. Basically, AD keeps track of the gradient of primitive operations and uses the chain rule to link them together. We'll use a package called `AutoGrad` to demo AD. Demonstration 1: $\sin(x)$<jupyter_code>anp.cos(1.0)
grad(anp.sin)(1.0)<jupyter_output><empty_output><jupyter_text>What happened above is that AutoGrad's `grad` takes in a Python function, and returns another Python function that computes the gradient.Some more examples:<jupyter_code>x = anp.random.rand(4)

def foo(x):
    return anp.sum(x)

def foo_grad(x):
    return anp.ones(len(x)) ## SOLUTION

foo_grad(x)
foo_grad_AD = grad(foo)
foo_grad_AD(x)
def pin(x):
    return anp.sin(x[1])

def pin_grad(x):
    grad = anp.zeros(len(x)) ## SOLUTION
    grad[1] = anp.cos(x[1])
    return grad

pin_grad(x)
pin_grad_AD = grad(pin)
pin_grad_AD(x)<jupyter_output><empty_output><jupyter_text>As you can see above, Autograd knows how to deal with things like `np.sum`, `np.abs`, subtraction, addition, etc. But let's make sure it wasn't cheating and just doing numerical differentiation!<jupyter_code>def pin(x):
    print("Evaluated function!") # let us know what it's doing
    return anp.sin(x[1])

def pin_grad(x):
    grad = anp.zeros(len(x)) ## SOLUTION
    grad[1] = anp.cos(x[1])
    return grad<jupyter_output><empty_output><jupyter_text>If it were doing numerical differentiation, it would need to evaluation the function $d+1$ times (5 times in this case), like this:<jupyter_code>scipy.optimize.approx_fprime(x, pin, 1e-6) # do numerical differentiation
pin_grad_AD = grad(pin) # automatic differentiation
pin_grad_AD(x)<jupyter_output>Evaluated function!<jupyter_text>Looks like it's not doing numerical differentiation! Bam!<jupyter_code>def baz(x):
    result = 0
    for i in range(len(x)):
        result = result + x[i]**i
    return result

def baz_grad(x):
    result = [0] ## SOLUTION
    for i in range(1,len(x)):
        result.append(i*x[i]**(i-1))
    return result

baz_grad(x)
baz_grad_AD = grad(baz) # differentiate through `for` loops!
baz_grad_AD(x)
def bar(x):
    if anp.abs(x[1]) > 2:
        return 0
    else:
        return -(x[0] * x[0] + 1) * anp.cos(x[1] - 1)


def bar_grad(x):
    if anp.abs(x[1]) > 2:  ## SOLUTION
        return anp.zeros(len(x))
    else:
        result = 0 * x
        result[0] = -2 * x[0] * anp.cos(x[1] - 1)
        result[1] = (x[0] * x[0] + 1) * anp.sin(x[1] - 1)
        return result


bar_grad(x)
bar_grad_AD = grad(bar) # differentiate through `if` statements!
bar_grad_AD(x)<jupyter_output><empty_output><jupyter_text>As an additional demo, let's try robust regression with the Huber loss (whose derivative I don't really feel like dealing with):<jupyter_code>d = 10
n = 1000

# generate random data
X = anp.random.randn(n, d)
w_true = anp.random.randn(d)
y = X @ w_true
# add random outliers
Noutliers = 50
y[:Noutliers] += 100 * anp.random.randn(Noutliers)
w_true
from sklearn.linear_model import HuberRegressor

hr = HuberRegressor(fit_intercept=False, alpha=0)
hr.fit(X, y)
hr.coef_
huber = lambda z: 0.5 * z ** 2 * (anp.abs(z) <= 1) + (anp.abs(z) - 0.5) * (anp.abs(z) > 1)
f = lambda w: anp.sum(huber(X @ w - y))<jupyter_output><empty_output><jupyter_text>Okay here we go:<jupyter_code>df_dw = grad(f) # differentiate through matrix multiplications, etc.
w = np.zeros(d)

alpha = 0.001
while anp.linalg.norm(df_dw(w)) > 0.0001:
    w -= alpha * df_dw(w)

w<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12198153409090909, 'fraction_numerical': 0.017134232954545456, 'mean_word_length': 4.610059760956175, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 57, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '151542', 'license': 'MIT', 'n_tokens_mistral': 3899, 'n_tokens_neox': 3644, 'n_words': 1352, 'path': 'lectures/appendixC_computing-derivatives.ipynb', 'repo_name': 'UBC-MDS/DSCI_572_sup-learn-2', 'size': '11270'}