<jupyter_start><jupyter_text>Salad DocumentationGroup 22: <NAME>, <NAME>, <NAME>, <NAME>, 2021---  Table of Contents* [Introduction](c1)* [Background](c2)    * [Computational Graph](c21)    * [Chain Rule](c22)    * [Forward Mode](c23)    * [Gradient & Jacobian](c24)    * [Dual Numbers](c25)* [Software Organization](c4)    * [Directory Structure](c41)    * [Modules](c42)    * [Test Suite](c43)    * [Packaging and Distribution](c44)* [How to Use Salad](c3)    * [Installation](c31)    * [Methods](c32)        * [1 input to 1 output](c321)        * [m inputs to 1 output](c322)        * [m inputs to n outputs](s1)* [Implementation](c5)    * [Data Structure](c51)    * [Classes and Methods](c52)        * [Variable](c521)        * [Forward](c522)* [Extensions](c6)    * [Gradient Descent (GD)](c61)    * [Broyden–Fletcher–Goldfarb–Shanno (BFGS)](c62)    * [Stochastic Gradient Descent (SGD)](c63)* [Broader Impact](c7)* [Software Inclusivity](c8)* [Licensing](c9)* [Future Directions](c10)  1 IntroductionDifferentiation is a central operation in science and engineering, used in various settings from financial markets to physics labs. For example, to optimize a function we must find the minima or maxima by using differentiation tests. There are multiple computational mechanisms to find derivatives. For example, numerical approximations like the Finite-Differences method involve utilizing the definition of a derivative and plugging in small values of h to evaluate the function. However, this method suffers from poor accuracy and stability issues due to floating point errors. Symbolic computation is another method, more accurate than numerical approximations, but often not applicable due to its unweildly and complex overhead. Here, we focus on automatic differentiation (AD), which overcomes the drawbacks of these two issues: it has less overhead than symbolic differentiation while calculating derivatives at machine precision. This makes it an integral part of machine learning algorithms, which demand computational efficiency as well as high precision.  2 BackgroundAutomatic differentiation can be approached in two ways, known as $\textit{forward mode}$ and $\textit{ reverse mode}$. Both methods depend on the decomposition of a complex function into elementary operations. Then, we can take advantage of the chain rule to compute the derivatives of complex functions from trivially computed derivatives of these elementary functions.  2.1 Computational GraphThe first step of AD is to construct a computational graph which relates the elementary inputs to the final function. Each node represents an intermediate computation $v_n$, and edges connect nodes via elementary functions $V_n$. For example, given a multivariate scalar function $f(x_1, x_2) = sin(x_1 + x_2)$, the computational graph looks like the following:  Elementary functions include but are not limited to the following: algebraic operations, exponentials, trigonometrics, and their inverses. AD takes advantage of the fact that derivatives of elementary functions are trivial to compute, allowing the differentiation of complex functions via the chain rule.  2.2 Chain RuleSuppose we want to calculate $\frac{\partial f}{\partial x_1}$. From the computational graph, we can rewrite $f$ in terms of nested intermediate functions:\begin{equation*}    f(x_1, x_2) = V_2(V_1(V_{-1}(x_1), V_0(x_2)))\end{equation*}where $V_n$ is the function that gives the value of node $v_n$. The chain rule states that\begin{equation*}    \frac{\partial f}{\partial x_1} = \frac{\partial V_2}{\partial V_1} \frac{\partial V_1}{\partial V_{-1}} \frac{\partial V_{-1}}{\partial x_1}\end{equation*}For a multivariate function $f(g_1(x), g_2(x), ..., g_n(x))$, the chain rule generalizes to:\begin{equation*}    \frac{\partial f}{\partial x} = \sum_{i = 1}^{n} \Big( \frac{\partial f}{\partial g_i} \frac{\partial g_i}{\partial x} \Big)\end{equation*}  2.3 Forward ModeForward mode, which is the conceptually simpler method, involves computing the chain rule from inside to outside. For the above example $f(x_1, x_2) = sin(x_1 + x_2)$, the $\textit{evaluation trace}$ at the point $(1, 2)$ is generated using the computational graph as a guide: The evaluation trace essentially computes a value (primal trace) and derivative (tangent trace) for each node of the computational graph. The desired final results are $\frac{\partial f}{\partial x_1}$ and $\frac{\partial f}{\partial x_2}$, shown in the last row of the trace using the gradient symbol $\nabla$. Forward mode is advantageous when there are few inputs ($m$) and many outputs ($n$). The opposite is true for reverse mode, which motivates its use in deep learning applications with many inputs.  2.4 Gradient & JacobianThe gradient of a multivariate scalar-valued function $f(g_1, g_2, ..., g_m)$ gives all of its partial derivatives in a vector of size $m \times 1$:\begin{equation*}    \nabla f = \nabla f_g = \begin{bmatrix}    \frac{\partial f}{\partial g_1} \\    \frac{\partial f}{\partial g_2} \\    \vdots \\    \frac{\partial f}{\partial g_m}    \end{bmatrix}\end{equation*}If the inputs $g_i$ are themselves multivariate functions with inputs $x \in \mathbb{R}^m$, the chain rule gives\begin{equation*}    \nabla_x f = \sum_{i = 1}^{n} \Bigg( \frac{\partial f}{\partial g_i} \nabla g_i (x) \Bigg)\end{equation*}Finally, we reach the most general form of the gradient: for multivariate vector-valued functions with $m$ inputs and $n$ outputs, the derivative of outputs with respect to inputs can be represented in a $\textit{Jacobian matrix}$ of dimensions $n \times m$:\begin{gather*}    f(x_1, x_2, ..., x_m) = \begin{bmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{bmatrix} \\    J_f    =    \begin{bmatrix} \frac{\partial y_1}{\partial x_1} & \ldots & \frac{\partial y_1}{\partial x_m} \\    \vdots & \ddots & \vdots \\    \frac{\partial y_n}{\partial x_1} & \ldots & \frac{\partial y_n}{\partial x_m} \end{bmatrix}\end{gather*}  2.5 Dual NumbersDual numbers offer a convenient way to store and compute the primal and tangent traces in forward mode AD. A dual number has the form\begin{equation*}    a + b \epsilon,\end{equation*}where $a, b \in \mathbb{R}$ and $\epsilon^2 = 0$. Combined with Taylor series expansion of a function $f$ centered at $a$, dual numbers give rise to the property\begin{equation*}    f(a + b \epsilon) = f(a) + f' (a) b \epsilon \end{equation*}because all terms with powers of $\epsilon$ higher than $2$ are equal to zero. Thus, the Taylor series approximation of a function about a dual number gives both its value and derivative, in a single calculation. In forward mode AD it is convenient to represent both the primal and tangent trace of each node as a dual number, where the real part corresponds to the primal and the dual part corresponds to the tangent. As shown below, functions (operations) on these nodes preserve the relationship between the real and dual parts:\begin{equation*}    g\big( f(a + \epsilon) \big) = g \big( f(a) + f' (a) \epsilon \big) = g \big( f(a) \big) + g' \big( f(a) \big) f'(a) \epsilon\end{equation*}Note that the coefficient of $\epsilon$ is simply the derivative of the real part $g(f(a))$, as given by the chain rule.  3 Software Organization  3.1 Directory StructureOur directory structure looks like the following:    cs107-FinalProject/    ├── LICENSE    ├── README.md    ├── pyproject.toml    ├── requirements.txt    ├── setup.py    ├── docs/    │   ├── images/    │   │   ├── comp_graph.png    │   │   └── table1.png     │   ├── milestone1.pdf    │   ├── milestone2 progress.pdf    │   ├── milestone2.pdf    │   ├── Documentation.ipynb    │   └── Documentation.html    ├── tests/    │   ├── run_test.sh    │   └── test_forward.py    │   └── test_optimize.py    └── src/        └── cs107_salad/            ├── __init__.py            ├── Forward/            │   ├── __init__.py            │   ├── salad.py            │   └── utils.py            └── Optimization/                ├── __init__.py                └── optimize.py  3.2 ModulesThree external modules that we will rely on in our implementation will be NumPy, pytest, and coverage.+ NumPy offers support for large, multi-dimensional arrays as well a wide range of optimized mathematical functions and operators for these arrays. Numerical calculation of function values in Salad will rely extensively on this module.+ pytest is a feature rich testing framework in python that makes it easy to write unit tests that scale for full libraries. + coverage is a commonly-used library used to compute code coverage during test execution.  3.3 Test SuiteThe test suite will live inside the $\texttt{tests}$ directory.To run the test suite, find $\texttt{run_tests.sh}$ in the $\texttt{tests}$ folder. Run the following command to run test and generate coverage report.```sh run_tests.sh```Local test result:                                    Module                statements                missing                excluded                coverage                                                        /home/ddv/projects/CS107/final-project/cs107-FinalProject/src/cs107_salad/Forward/__init__.py                0                0                0                100%                                        /home/ddv/projects/CS107/final-project/cs107-FinalProject/src/cs107_salad/Forward/salad.py                380                31                0                92%                                        /home/ddv/projects/CS107/final-project/cs107-FinalProject/src/cs107_salad/Forward/utils.py                28                2                0                93%                                        /home/ddv/projects/CS107/final-project/cs107-FinalProject/src/cs107_salad/Optimization/__init__.py                0                0                0                100%                                        /home/ddv/projects/CS107/final-project/cs107-FinalProject/src/cs107_salad/Optimization/optimize.py                128                5                0                96%                                        /home/ddv/projects/CS107/final-project/cs107-FinalProject/src/cs107_salad/__init__.py                0                0                0                100%                                        test_forward.py                984                29                0                97%                                        test_optimize.py                95                3                0                97%                                                        Total                1615                70                0                96%                          Test coverage is currently at 96%.  3.4 Packaging and DistributionWe have distributed the package through the Test Python Package Index (PyPi). Users will be able to install our package just by doing:```pip install cs107-salad --extra-index-url=https://test.pypi.org/simple/```Because the package is on test PyPi, the extra index url flag allows numpy, one of package dependecies, to be installed from main PyPi instead of test PyPi. Test PyPi does not have the most up to date version of numpy.Eventually, when the package has even more features, we would aim to distribute the package on the main PyPi.  4 How to Use $\texttt{Salad}$  4.1 InstallationThe users will pull $\texttt{Salad}$ from PyPI using command line<jupyter_code>!pip install cs107-salad --extra-index-url=https://test.pypi.org/simple/<jupyter_output>Looking in indexes: https://pypi.org/simple, https://test.pypi.org/simple/
Collecting cs107-salad
  Downloading https://test-files.pythonhosted.org/packages/52/7b/66409627c663197c0d185dfcf992744bffb2d86982f9aa496a9654854080/cs107_salad-1.0.1-py3-none-any.whl (12 kB)
Requirement already satisfied: numpy in /usr/local/lib/python3.7/dist-packages (from cs107-salad) (1.19.5)
Installing collected packages: cs107-salad
Successfully installed cs107-salad-1.0.1<jupyter_text>Then import $\texttt{salad}$ and other dependencies (e.g. $\texttt{numpy})$ in python script<jupyter_code>import cs107_salad.Forward.salad as ad
import numpy as np<jupyter_output><empty_output><jupyter_text>4.2 MethodsUsers can use our package to deal with 3 different scenarios: 1 input to 1 output; multiple inputs to 1 output; multiple inputs to multiple outputs.  4.2.1 1 input to 1 outputUsers may specify a single function with one input variable. See the example below with $f(x) = \exp(2x)$ with $x=2$.<jupyter_code>x = ad.Variable(2, label="x") # specifies the input variable
f = ad.exp(2 * x) # specifies functions involving the variable<jupyter_output><empty_output><jupyter_text>We make two observations at this stage: First, arithmetic operations ($\exp$ in our example) need to be implemented within $\texttt{salad}$. Second, $\texttt{f}$ is of type $\texttt{Variable}$.To retrieve the value and derivative, call<jupyter_code>val = f.val # array(exp(4)), numpy.array
val
der = f.der["x"] # 2exp(4), float
der<jupyter_output><empty_output><jupyter_text>The input value of the variable can also be a list consisting of multiple points. See the example below with $f(x) = \exp(2x)$ with $x=1,2,3,4$.<jupyter_code>x = ad.Variable([1,2,3,4], label="x") # specifies the input variable consisting of 4 points
f = ad.exp(2 * x) # specifies functions involving the variable<jupyter_output><empty_output><jupyter_text>The observations are stored in an array of the same length.<jupyter_code>val = f.val # value: exp(2x)
val
der = f.der["x"] # derivatives: 2exp(2x)
der<jupyter_output><empty_output><jupyter_text>4.2.2 $m$ inputs to 1 outputWhen there are multiple input variables, users can specify the name of the variable when initializing the instance. The name usually matches the variable name on the left hand side.  See the example below with $f(x, y) = \exp(x+2y)$, with $x=1, y=1$.<jupyter_code>x = ad.Variable(1, label='x') # specifies input variable x
y = ad.Variable(1, label='y') # specifies input variable y
f = ad.exp(x + 2 * y) # specifies functions involving the variables
val = f.val # array(exp(3)), numpy.array
val<jupyter_output><empty_output><jupyter_text>To get the derivative with respect to $x$, one simply uses "$x$" as the key to get the value from the dict $\texttt{f.der}$. Calling $\texttt{f.der}$ gives the gradient with respect to all input variables inside a dictionary.<jupyter_code>dx = f.der["x"] # exp(3), float
dx
dy = f.der["y"] # 2exp(3), float
dy
grad = f.der # {"x": exp(3), "y": 2exp(3)}, dict
grad<jupyter_output><empty_output><jupyter_text>Similarly, the users can calculate the values and derivatives of the function at multiple data points by using a list as input. See the example below with $f(x, y) = \exp(x+2y)$ at $(1, 1)$ and $(3, 0)$<jupyter_code>x = ad.Variable([1,3], label='x') # specifies input variable x
y = ad.Variable([1,0], label='y') # specifies input variable y
f = ad.exp(x + 2 * y) # specifies functions involving the variables
val = f.val # array(exp(3), exp(3)), numpy.array
val
grad = f.der # {"x": [exp(3), exp(3)], "y": [2exp(3), 2exp(3)]}, dict
grad<jupyter_output><empty_output><jupyter_text>4.2.3 $m$ inputs to $n$ outputsWe recognize that users may sometimes want to compute Jacobians. Consider the following case, $x = [x_1~x_2]^T, y = [y_1~y_2]^T, f(x, y) = [2x + \exp(y) ~ 3x+2\sin(y)]^T$.We provide another class $\texttt{Forward}$ to do the calculation.<jupyter_code>variables = {'x': 3, 'y': 5} # specifies input values
functions = ['2*x + exp(y)', '3*x + 2*sin(y)'] # specifies functions
f = ad.Forward(variables, functions)
print(f)<jupyter_output>Function: 2*x + exp(y), Value: 154.4131591025766, Derivative: {'x': array(2.), 'y': 148.4131591025766}
Function: 3*x + 2*sin(y), Value: 7.082151450673723, Derivative: {'x': array(3.), 'y': array(0.56732437)}<jupyter_text>In this case, all computed results will be saved in $\texttt{f.results}$, which is a list of values and derivatives of the two functions. The user can obtain the value and the derivatives of the $i^{th}$ function by calling<jupyter_code>val1 = f.results[0].val # array(11)
val1
der1 = f.results[0].der # {'x': array(2.), 'y': array(1.)}
der1
val2 = f.results[1].val # array(19)
val2
der2 = f.results[1].der # {'x': array(3.), 'y': array(2.)}
der2<jupyter_output><empty_output><jupyter_text>The input can also be a list of values.<jupyter_code>variables = {'x': [3,4], 'y': [2,3]} # specifies input values
functions = ['2*x + exp(y)', '3*x + 2*sin(y)'] # specifies functions
f = ad.Forward(variables, functions)
print(f)<jupyter_output>Function: 2*x + exp(y), Value: [13.3890561  28.08553692], Derivative: {'x': array([2., 2.]), 'y': array([ 7.3890561 , 20.08553692])}
Function: 3*x + 2*sin(y), Value: [10.81859485 12.28224002], Derivative: {'x': array([3., 3.]), 'y': array([-0.83229367, -1.97998499])}<jupyter_text>5 Implementation  5.1 Data StructureThe core data structures used in our implementation are $\texttt{list}$, $\texttt{numpy.array}$, and $\texttt{dict}$. See the [m inputs to n outputs](s1) example from above as an illustration. The input variables are in a $\texttt{dict}$. The input functions are in a $\texttt{dict}$. The Jacobian with respect to each vector-form input is in $\texttt{numpy.array}$. To easily match the Jacobian to the corresponding input vector, we use a $\texttt{dict}$, where the key is the name of the corresponding input variable, and the value is a $\texttt{numpy.array}$ of the Jacobian.  5.2 Classes and MethodsThe baseline implementation currently includes two classes: 1) $\texttt{Variable}$, 2) $\texttt{Forward}$  5.2.1 $\texttt{Variable}$In $\texttt{Variable}$, we re-define dunder methods and other arithmetic operations. See below for the structure of the class and method definition.<jupyter_code>class Variable(object):
    counter = 0

    def __init__(self, val, der=None, label=None, ad_mode="forward", increment_counter=True):
        pass

    def __add__(self, other):
        pass
    def __str__(self):
        pass
    def __radd__(self, other):
        pass
    def __sub__(self, other):
        pass
    def __rsub__(self, other):
        pass
    def __mul__(self, other):
        pass
    def __rmul__(self, other):
        pass
    def __truediv__(self, other):
        pass
    def __rtruediv__(self, other):
        pass
    def __neg__(self, other):
        pass
    def __pow__(self, other):
        pass
    def __rpow__(self, other):
        pass
    def __eq__(self, other):
        pass
    def __ne__(self, other):
        pass
    def __lt__(self, other):
        pass
    def __le__(self, other):
        pass
    def __gt__(self, other):
        pass
    def __ge__(self, other):
        pass<jupyter_output><empty_output><jupyter_text>On a high level, each instance of this class corresponds to a node in the computation graph in the background section. We keep track of the primal trace ($\texttt{self.val}$) and the tangent trace ($\texttt{self.der}$) as attributes for all instances.We also implement elementary functions like $\texttt{sin}$, $\texttt{sqrt}$, $\texttt{log}$, and $\texttt{exp}$.<jupyter_code>def sin(x):
    pass
def cos(x):
    pass
def tan(x):
    pass
def exp(x):
    pass
def log(x):
    pass
def ln(x):
    pass
def sqrt(x):
    pass
def logistic(x):
    pass
def arcsin(x):
    pass
def arccos(x):
    pass
def arctan(x):
    pass
def sinh(x):
    pass
def cosh(x):
    pass
def tanh(x):
    pass<jupyter_output><empty_output><jupyter_text>Users can use these functions for simple numeric calculations:<jupyter_code>x = 3
ans = ad.sin(x)
print(f'x = {ans}, \ntype of x: {type(ans)}')<jupyter_output>x = 0.1411200080598672, 
type of x: <class 'numpy.float64'><jupyter_text>Or, users can use these functions to compute function derivatives. In this case, the inputs should be of the type $\texttt{Variable}$.<jupyter_code>x = ad.Variable(3, label="x")
ans = ad.sin(x)
print(f'x = {ans.val}, x.der = {ans.der}')
print(type(ans))<jupyter_output><class 'cs107_salad.Forward.salad.Variable'><jupyter_text>More arithemetic operations and functions to be included if time allows and as we see fit. 5.2.2 $\texttt{Forward}$We also include a $\texttt{Forward}$ class for easy handling of multiple function input. The variables are input in a $\texttt{dict}$ and the functions are input in a $\texttt{list}$. The $\texttt{__init__}$ function within $\texttt{Forward}$ simply loops through and evaluate each functions one by one. See the [m inputs to n outputs](s1) example from above as an illustration.  6 ExtensionsWe implement an optimization toolkit as the extension feature. These optimization methods are useful in numerous contexts such as statistics, machine learning, operations research, and economics. The `Optimize` module includes the following optimization schemes:1. Gradient Descent (GD)2. Broyden-Fletcher-Goldfarb-Shanno (BFGS)3. Stochastic Gradient Descent (SGD)The intended use of the first two methods, GD and BFGS is as follows:1. users input an arbitrary function supported by `Salad`, $f(\textbf{x})$, where $\textbf{x} \in \mathbb R^n$2. users supply a starting value, $\textbf{x}_0 \in \mathbb R^n$, at which the method starts to search for a local minimum3. the optimization method returns $\mathbf{x}^* \in \mathbb R^n$ that achieves a local minimum of the function $f(x)$.The third method, SGD, is implemented in the context of a linear regression, $\mathbf{y}$ regressed on $\mathbf{X}$.The intended uses of SGD is as follows:1. users input covariate matrix $\mathbf{X} \in \mathbb R^{nxp}$ with $n$ rows ($n$ observations) and $p$ columns ($p$ covariates); and input response vector $\mathbf{y} \in \mathbb R^n$ with $n$ entries ($n$ observations).2. users supply the batch size used in each SGD iteration 3. user specify the starting position of the search for coefficient estimates4. SGD returns the regression coefficient $\mathbf{b}^* \in \mathbb R^b$. Users first need to import the optimization toolkit at `cs107_salad.Optimization.optimize`.<jupyter_code>import cs107_salad.Optimization.optimize as optimize<jupyter_output><empty_output><jupyter_text>6.1 [Gradient Descent (GD)](https://en.wikipedia.org/wiki/Gradient_descent)Gradient descent is an iterative optimization algorithm to find a local minimum of some function. The method iteratively calculates the gradient at each given point, and goes toward the direction of steepest descent. Specifically, consider a differentiable multi-variable function $f(\textbf{x})$. Starting with an initial guess $\textbf{x}_0$, at each step $k$, the method updates $x_{k+1}$ in the following way:$$\textbf{x}_{k+1} = \textbf{x}_k - \lambda \nabla F(\textbf{x}_k)$$The algorithm terminates when the magnitude of the update $\lambda \nabla F(\textbf{x}_k)$ is smaller than a tolerance level or the maximum number of iterations is reached. $\lambda$ is the learning rate that adjusts for how much the method goes in the negative gradient direction. Below is an demonstration of calling GD to find the minimin of a convex function $x^2 + 2y^2$. `min_params` is the $\mathbf{x}$ value that gives a local minimum, `val` is the function value at this local minimum, and `der` is the gradient at this local minimum.<jupyter_code>f = "x**2 + 2*y**2"
starting_pos = {"x": 5, "y": 2}

GD = optimize.GradientDescent()
min_params, val, der, hist = GD.optimize(f, starting_pos, full_history=True)
print('optimized parameter:', min_params)
print('function value at minimum', val)
print('gradient at minimum:', der)<jupyter_output>optimized parameter: {'x': 4.8769661418382606e-09, 'y': 1.240470176258791e-18}
function value at minimum 2.3784798748636768e-17
gradient at minimum: {'x': 9.952992126200531e-09, 'y': array(5.16862573e-18)}<jupyter_text>Users may ask `GD.optimize` to output the search history by including a `full_history` flag. The full history allows users to easily visualize the search path.<jupyter_code>import matplotlib.pyplot as plt

def f(x):
  x, y = x
  return x**2 + y**2

def plot_f(x,y,canvas, title):

    X, Y = np.meshgrid(x,y)
    Z = f([X,Y])
    canvas.contour(X,Y,Z,levels=np.linspace(0,100,80), alpha=0.3)
    canvas.plot([0],[0],'ro',lw=2, label='global min')
    canvas.set_xlabel('x')
    canvas.set_ylabel('y')
    canvas.set_title(title)
    canvas.legend(loc=3)
_, axs = plt.subplots(1,1,figsize=(8,6))
x = np.linspace(-6,6,100)
y = np.linspace(-6,6,100)
plot_f(x,y,axs,title='Gradient Descent Path')
axs.plot(hist['x'][0], hist['y'][0], 'bo', lw=1, label='starting position');
axs.plot(hist['x'], hist['y'], 'b--', lw=1, label='search path');
axs.legend();<jupyter_output><empty_output><jupyter_text>6.2 [Broyden–Fletcher–Goldfarb–Shanno (BFGS)](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm)Broyden–Fletcher–Goldfarb–Shanno (BFGS) is a popular quasi-Newton Method that iteratively finds a local minimum of some functions using approximated Hessian. Specifically, consider a differentiable multi-variable function $f(\textbf{x})$. BFGS starts with an initial guess $\textbf{x}_0$, and $\mathbf{B}_0 = \mathbf{I}$ before the first iteration. The algorithm computes the following in iteration $k$:1. Solve for $\mathbf{B}_k$ that satisfies $\mathbf{B}_k \mathbf{s}_k = -\nabla f(\mathbf{x}_k)$2. Update $\mathbf{x}_{k+1} = \mathbf{x}_k + \mathbf{s}_k$3. Compute $\mathbf{y}_k = \nabla f(\mathbf{x}_{k+1}) - \nabla f(\mathbf{x}_k)$4. Update $\mathbf{B}_{k+1} = \mathbf{B}_k + \Delta \mathbf{B}_k$where $\Delta B_k = \frac{\textbf{y}_{k}\textbf{y}_{k}^T}{\textbf{y}_{k}^T\textbf{s}_{k}} - \frac{B_k\textbf{s}_{k}\textbf{s}_{k}^TB_k}{\textbf{s}_{k}^TB_k\textbf{s}_{k}}$The algorithm terminates when the update $\mathbf{s}_k$ is smaller than a tolerance or if the maximum number of iterations has been reached.  The example below finds the local (global) minimum for a convex function $x^2 + 2y^2$. `min_params` stores the $\mathbf{x}$ value that attains the local (global) minimum, `val` stores the function value at the minimum, and `hist` stores the entire search path of $\mathbf{x}$.<jupyter_code>f = "x**2 + 2*y**2"
starting_pos = {"x": 5.0, "y": 2.0}

BFGS = optimize.BFGS()
min_params, val, hist = BFGS.optimize(f, starting_pos, full_history=True)
print('optimized parameter:', min_params)
print('function value at minimum', val)
import matplotlib.pyplot as plt

def f(x):
  x, y = x
  return x**2 + y**2

def plot_f(x,y,canvas, title):

    X, Y = np.meshgrid(x,y)
    Z = f([X,Y])
    canvas.contour(X,Y,Z,levels=np.linspace(0,100,80), alpha=0.3)
    canvas.plot([0],[0],'ro',lw=2, label='global min')
    canvas.set_xlabel('x')
    canvas.set_ylabel('y')
    canvas.set_title(title)
    canvas.legend(loc=3)
_, axs = plt.subplots(1,1,figsize=(8,6))
x = np.linspace(-6,6,100)
y = np.linspace(-6,6,100)
plot_f(x,y,axs,title='Gradient Descent Path')
axs.plot(hist['x'][0], hist['y'][0], 'bo', lw=1, label='starting position');
axs.plot(hist['x'], hist['y'], 'bx--', lw=1, label='search path');
axs.legend();<jupyter_output><empty_output><jupyter_text>6.3 [Stochastic Gradient Descent (SGD)](https://en.wikipedia.org/wiki/Stochastic_gradient_descent)Stochastic gradient descent (SGD) algorithm iteratively finds a local minimum of some functions with suitable smoothness properties. It is a more generalized version of gradient descent that involves randomness in gradient calculation. In each iteration of SGD,  it randomly samples a subset (batch) of the entire input dataset, and conducts one step of gradient descent using the gradient calculated from the subset to reduce computational burden. SGD terminates when the maximum number of iterations is reached. Our implementation of SGD has a limited use case: linear regression. In particular, users specify a covariate matrix $\mathbf{X}$ and a response vector $\mathbf{y}$, our SGD function aims to minimize the following function:$$\mathbf{b}^* = \arg\min_{\mathbf{b}} = \frac{1}{n} \sum_{i=1}^n \left(y_i - \mathbf{X}_i \mathbf{b}\right)^2$$where $\mathbf{X} \in \mathbb R^{n\times p}, \mathbf{y} \in \mathbb R^n, \mathbf{b}\in \mathbb R^p$. Note that the objective as a function of the estimated coefficients $\mathbf{b}$ is convex. Therefore, SGD in this case actually leads to the global minimum. $\mathbf{b}^*$ is the estimated coefficient by minimizing mean squared error.  The folloing example is to calculate the coefficient of regressing $\mathbf{y}$ on $\mathbf{X}$. `min_params` is the estimated cofficients, `val` is the loss attained at the estiamted coefficients, and `der` is the gradient at the minimum loss.<jupyter_code>X = np.random.rand(100, 3)
y = X @ np.array([2, 0, 3]) + 0.001 * np.random.rand(X.shape[0]) # ground truth data generating process 

SGD = optimize.StochasticGradientDescent(X, y, batch_size=10)
min_params, val, der = SGD.optimize([0, 0, 0], max_iter=5000, learning_rate=0.01)
print('optimized parameter:', min_params)<jupyter_output>optimized parameter: {'b0': 2.0003212074294807, 'b1': 0.0006186558704257924, 'b2': 2.999974911074871}<jupyter_text>7 Broader ImpactOur extension involves several optimization methods: gradient descent, stochastic gradient descent, BFGS, and Newton’s method. By itself, these are facially neutral algorithms – they are not discriminatory on their face – yet they have clear potential discriminatory applications and effects. Our subpackages could easily be extended and combined with other software modules to create bad, biased, or unethical AI models. For instance, Amazon was forced to scrap its automated recruiting tool when it was revealed that the algorithm was biased against women. The algorithm was biased because it was trained on resumes previously submitted to Amazon. However, reflecting unequal trends pervasive throughout the tech industry, most candidates that made it through the recruitment process were men. Amazon’s model, when optimized to this data, learned that resumes from men were preferable to resumes from women. One could imagine how our package, when used to optimize on biased data from any setting would have negative disparate impacts.  8 Software InclusivityThe salad package warmly welcomes all people regardless of background to contribute to the development of our package. Our package is open-source, easily accessible through Github and PyPi, and licensed under the MIT License. Nonetheless, we recognize there are many structural barriers for the average person to interact with and contribute to our software package. For instance, while we aimed to have high code quality with clear comments, our comments are all in English. As a result, non-native English speakers are at a disadvantage in contributing to our package. Moreover, computer science/programming education and education in general are not shared equally across society. BIPOC, the poor, women, and other minorities are all not afforded the same educational opportunities as the cis-male White and Asian men that currently make up the majority of programmers in the US. Such underrepresented minorities face more constraints in contributing to open source software projects like ours.    9 LicensingWe will use the MIT License. Future developers are welcome to extend our library and distribute closed source versions. Because we are using NumPy, a copyright licensed library, we are unable to use a license like GNU GPLv3 (as it requires all libraries used by the library to be copyleft).  10 Future DirectionsThere is much room for improvement of our package from supporting higher and mixed derivatives to building out support for differential programming.The four areas that we are most interested in, however, are the following:1. Support for reverse mode. Our package currently supports computing derivatives to machine precision of functions $f:\mathbb{R}^m\to\mathbb{R}^n$. However, forward mode is most efficient when $n \gg m$, ie the number of functions we are evaluating is much greater than the number of inputs. When $n \ll m$, i.e. the number of functions we are evaluating is much smaller than the number of inputs, reverse mode is actually more efficient.    - Support for backpropagation. Backpropagation is a special case of reverse mode that is used for scalar objective functions where our objective function represents the error between our output and the true value. Implementing backpropagation would allow us to extend our package to be used for developing fully connected neural networks.2. Building out support for more elementary functions such as the inverse hyperbolic functions.3. Building out more user friendly tools such as a GUI that allows the user to visualize the computational graph of the functions they have entered.4. The optimization toolkit is also preliminary at the moment. A couple ideas for future extensions include:  a. Implementation of the Newton's Method: Newton's Method for finding the local minimum of a function uses the Hessian, which requires a further extension to our `salad` implementation. We could include the computation of the Hessian in our implementation, potentially as a third variable within the `Variable` class.  b. Improvement of the current implementation of SGD: Currently, the only use case of our SGD is for regression (i.e. minimize the mean squared error). The method could benefit from  - Allowing users to specify any loss function of their choice  - Automic hyper-parameter tuning   - Repeated search from multiple independent starting point to achieve a lower local minimum  c. Implementation of other optimization fucntions such as `Adam`, `RMSprop`, and `FTRL` which are more stable and more efficient optimizers commonly used for deep learning and neural networks.<jupyter_code><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10882536289963488, 'fraction_numerical': 0.028349214830646837, 'mean_word_length': 3.9519329707481994, 'pattern_counts': {'":': 8, '<': 88, '<?xml version=': 0, '>': 88, 'https://': 8, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 41, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '22', 'dir': 'jupyter-structured-clean-dedup', 'id': '505702', 'license': 'MIT', 'n_tokens_mistral': 10254, 'n_tokens_neox': 9589, 'n_words': 3924, 'path': 'docs/Documentation.ipynb', 'repo_name': 'cs107-DaLand/cs107-FinalProject', 'size': '33714'}