## Milestone 2

## Introduction

Differentiation, or the process of finding a derivative, is an extremely important mathematical operation with a wide range of applications. The discovery of extrema or zeros of functions is essential in any optimization problem, and the solving of differential equations is fundamental to modern science and engineering. Differentiation is essential in nearly all quantitative disciplines: physicists may take the derivative of the displacement of a moving object with respect to time in order to find the velocity of that object, and data scientists may use derivatives when optimizing weights in a neural network. 

Naturally, we would like to compute the derivative as accurately and efficiently as possible. Two classical methods of calculating the derivative have clear shortcomings. Symbolic differentiation (finding the derivative of a given formula with respect to a specified variable, producing a new formula as its output) will be accurate, but can be quite expensive computationally. The finite difference method ($\frac{\partial f}{\partial x} = \frac{f(x+\epsilon)-f(x)}{\epsilon}$ for some small $\epsilon$) does not have this issue, but will be less precise as different values of epsilon will give different results. This brings us to automatic differentiation, a less costly and more precise approach.

## Background

Automatic differentiation breaks down the main function into elementary functions, evaluated upon one another. It then uses the chain rule to update the derivative at each step and ends in the derivative of the entire function.

To better understand this process, let's look at an example. Consider the example function

\begin{equation}
f(x) = x + 4(sin(\frac{x}{4}))
\end{equation}

We would like to compute the derivative of this function at a particular value of x. Let's say that in this case, we are interested in evaluating the derivative at $x=\pi$. In other words, we want to find $f'(\pi)$ where $f'(x) = \frac{\partial f}{\partial x}$

We know how to solve this _symbolically_ using methods that we learned in calculus, but remember, we want to compute this answer as accurately and efficiently as possible, which is why we want to solve it using automatic differentiation. 

To solve this using automatic differentiation, we need to find the decomposition of the differentials provied by the chain rule. Remember, the chain rule is a formula for computing the derivative of the composition of two or more functions. So if we have a function $h\left(u\left(t\right)\right)$ and we want the derivative of $h$ with respect to $t$, then we know by the chain rule that the derivative is $\dfrac{\partial h}{\partial t} = \dfrac{\partial h}{\partial u}\dfrac{\partial u}{\partial t}.$ The chain rule can also be expanded to deal with multiple arguments and vector inputs (in which case we would be calculating the _gradient)_.

Our function $f(x)$ is composed of elemental functions for which we know the derivatives. We will separate out each of these elemental functions, evaluating the derivative at each step using the chain rule. 

The evaluation trace for this problem looks like:


| Trace    | Elementary Operation &nbsp;&nbsp;&nbsp;| Derivative &nbsp;&nbsp;&nbsp; | $\left(f\left(a\right),  f^{\prime}\left(a\right)\right)$&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;|
| :------: | :----------------------:               | :------------------------------: | :------------------------------: |
| $x_{1}$  | $\pi$                      | $1$                | $\left(\pi, 1\right)$ |
| $x_{2}$  | $\frac{x_{1}}{4}$                               | $\frac{\dot{x}_{1}}{4}$                 | $\left(\dfrac{\pi}{4}, \frac{1}{4}\right)$ |
| $x_{3}$  | $\sin\left(x_{2}\right)$               | $\cos\left(x_{2}\right)\dot{x}_{2}$            | $\left(\dfrac{\sqrt{2}}{2}, \frac{\sqrt{2}}{8}\right)$ |
| $x_{4}$  | $4x_{3}$                               | $4\dot{x}_{3}$                  | $\left(2\sqrt{2}, \frac{\sqrt{2}}{2}\right)$ |
| $x_{5}$  | $x_{1} + x_{4}$                        | $\dot{x}_{1} + \dot{x}_{4}$ | $\left(\pi + 2\sqrt{2}, 1 + \frac{\sqrt{2}}{2}\right)$ |


This evaluation trace provides some intuition for how automatic differentiation (in the forward mode) is used to calculate the derivative of a function evaluated at a certain value ($f'(\pi) = 1 + \frac{\sqrt{2}}{2}$). 

You may notice that when we computed the derivative above, we "seeded" the derivative with a value of 1. This seed vector doesn't have to be 1, but the utility of using a unit vector becomes apparent when we consider a problem involving directional derivatives.

The definition of the directional derivative (where p is the seed vector)
$$D_{p}x_{3} = \sum_{j=1}^{2}{\dfrac{\partial x_{3}}{\partial x_{j}}p_{j}}$$

can be expanded to

\begin{align}
  D_{p}x_{3} &= \dfrac{\partial x_{3}}{\partial x_{1}}p_{1} + \dfrac{\partial x_{3}}{\partial x_{2}}p_{2} \\
             &= x_{2}p_{1} + x_{1}p_{2}
\end{align}

If we choose p to be a the unit vector, we can see how this is beneficial:

$p = \left(1,0\right)$ gives $\dfrac{\partial f}{\partial x}$

$p = \left(0,1\right)$ gives $\dfrac{\partial f}{\partial y}$

So to summarize, the forward mode of automatic differentiation is really computing the _**product of the gradient of our function with the seed vector:**_

$$D_{p}x_{3} = \nabla x_{3}\cdot p.$$ 

If our function is a vector, then the forward mode actually computes $Jp$ where $J = \dfrac{\partial f_{i}}{\partial x_{j}}, \quad i = 1,\ldots, n, \quad j = 1,\ldots,m$ is the Jacobian matrix. Often we will really only want the "action" of the Jacobian on a vector, so we will just want to compute the matrix-vector product $Jp$ for some vector $p$.

Automatic differentiation can be used to compute derivatives to machine precision of functions $f:\mathbb{R}^{m} \to \mathbb{R}^{n}$

The forward mode is more efficient when $n\gg m$.
   - This corresponds to the case where the number of functions to evaluate is much greater than the number of inputs.
   - Actually computes the Jacobian-vector product $Jp$.
   
The reverse mode is more efficient when $n\ll m$.
   - This corresponds to the case where the number of inputs is much greater than the number of functions.
   - Actually computes the Jacobian-transpose-vector product $J^{T}p$.

## How to Use autodiff32
### Installation

**1) Create a virtual environment (optional)**

From the terminal, create a virtual environment:

_(The command below will create the virtual environment in your present working directory, so consider moving to a project folder or a known location before creating the environment)_

```virtualenv env```

activate the virtual environment:

```source env/bin/activate```

if you plan to launch a jupyter notebook using this virtual environment, run the following to install and set up jupyter in your virtual environment:

```python -m pip install jupyter```

```python -m ipykernel install --user --name=env```

**3) Install the autodiff32 package**

In the terminal, type:

```pip install autodiff32```

Package dependencies will be taken care of automatically!

_(Alternatively, it is also possible to install the autodiff32 package by downloading this GitHub repository. If you choose that method, use the requirements.txt file to ensure you have installed all necessary dependencies.)_

### How to use

It is easy to use the autodiff32 package in a Jupyter notebook, as we will demonstrate here:

_(Alternatively, you can start a Python interpreter by typing ```Python``` into the terminal, or work from your favorite Python IDE.)_

_(Remember, if you are using a virtual environment, follow steps 1 through 3 above and then type ```jupyter notebook``` into your terminal to launch a notebook. Inside the notebook, switch the kernel to that of your virtual environment.)_


```python
import autodiff32 as ad
import math # only necessary for this particular example

"""
Initialize an AutoDiff object  
with the number you would like to pass into your function
"""
X = ad.AutoDiff(math.pi)

"""
Define your function of interest
"""
func = X + 4*ad.sin(X/4)
"""
Look at the derivative of your function 
evaluated at the number you gave above
"""
print("derivative:",func.der)

"""
Look at the value of your function 
evaluated at the number you gave:
"""
print("value:", func.val)
```

    derivative: 1.7071067811865475
    value: 5.970019778335983


Notice that this is the same equation used in our example above: $f(x) = x + 4(sin(\frac{x}{4}))$. Just for fun, let's see if the derivative that we calculated in the evolution trace is the same as the result using autodiff32:


```python
print("autodiff32 derivative:", func.der)
print("evolution trace derivative:", 1+math.sqrt(2)/2)
```

    autodiff32 derivative: 1.7071067811865475
    evolution trace derivative: 1.7071067811865475


We can see that the derivative calculated using autodiff32 is the same as the derivative calulated by walking manually through the evolution trace!

## Software Organization
### Directory Structure
Our structure is as follows:

    /cs207-FinalProject
        README.md
        LICENSE
        .gitignore
        .travis.yml
        setup.py
        requirements.txt
        docs/
            milestone1.ipynb
            milestone2.ipynb
        autodif32/
            __init__.py
            AutoDiffObj.py
            Elementary.py
            JacobianVectorFunc.py
            MultivariateVarCreator.py
        test/
            __init__.py
            autodiffobj_test.py
            elementary_test.py

### Basic Modules

The ```AutoDiffObj``` module creates an AutoDiff object based on the scalar value you would like to evaluate a function and its derivative at. It overloads the basic operations including multiply, add, negative, subtract, powers, division, and equality. It also includes a Jacobian method, which returns the Jacobian of the function. If the function is univariate, then the Jacobian is just the derivative. If the function is multivariate (not yet implemented), then the Jacobian will be an array.

The ```Elementary``` module implements some of the elementary functions, including exp, log, sqrt, sin, cos, tan, asin, acos, atan, sinh, cosh and tanh.

The ```JacobianVectorFunc``` module implements the Jacobian method for vector inputs.

The ```MultivariateVarCreator``` module takes in the values of each variable in a user defined multivariate function, and returns len(kwargs) number of AutoDiff class variables with derivative (seed) as an np.array. The MultivariateVarCreator class acts as a helper function for user to create multiple AutoDiff Objects conveniently (instead of manually create many AutoDiff objects) to use in the evaluation of the multivariate function. Please note that the implementation of multivariate functions is still in progress.

### Testing Suite
Our testing files live in the `test/` directory. The tests are run using pytest.


## Implementation details

The current implementation of AutoDiff32 allows for scalar univariate inputs to functions. Core classes include the AutoDiff class, . AutoDiff32 is externally dependent on numpy, and this dependency has been automatically taken care of in the released version of the package (as well as in the requirements.txt file if the user chooses to manually download the package). AutoDiff32 has implemented a number of elementary functions, as listed above in the description of the basic modules. 

We plan to continue development of the AutoDiff32 package to allow for vector inputs as well as multivariate inputs. This will require robust Jacobian matrix functionality. 

In addition to the currently implemented forward mode, we plan to build out the reverse mode of auto differentiation as an advanced feature.

DETAILS:


The core data structure (and also external dependencies) for our implementation will be numpy arrays, and the core classes we have implemented are described below:


1) An ***AutoDiff class*** which stores the current value and derivative for the current node. The class method will conatin overloading operators such as plus, minus, multiply, etc.

```python
"""
Initialize an Automatic Differentiation object which stores its current value and derivative

Note that the derivative needs to not be a scalar value

For multivariable differentiation problems of a scalar function, the derivative will be a vector
"""

def __init__(self,value,der=1) :
  self.val = val
  self.der = der
  
def check_dimension(self):
  """
  This method will be integrated in other class methods as well to check the dimensions when performing vector computations and catch appropriate errors
  """
  """
  This method has not been implemented yet
  """

#overloading operators to enable basic operations between classes and numbers (not exhaustive):

"""
These methods will differentiate cases in which other is a scalar, a vector, a class, or any child class

All of these operators will return AutoDiff classes
"""
def __mult__(self,other):
def __rmult__(self,other):
def __radd__(self,other): 
def __add__(self,other):
  
# compute Jacobian for the multivariate function
def Jacobian(self):
  return self.der


#For univairate functions
X = AutoDiff(3)
func = 3*X + 2
print(func.val)
print（func.der)

```

Multivariate functions (both scalar and vector) can be also evaluated using the AutoDiff class as below, and the resulting Jacobian will be a vector array:

```python
X = AutoDiff(1,np.array([1,0]))
Y =  AutoDiff(1,np.array([0,1]))

func = X + 2*Y
```

While this way of defining and evaluating multivariate functions is feasible, it is very inconvenient for the users to have to keep track of the dimensionality of the derivatives. For,example, the above func definition will raise an error if Y's derivative is defined as np.array([0,0,1]). This potential problem in dimensionality will also be a cause difficulties in the error handling process in the code. 

The way we tackle this problem is to create a helper class called **Multi_AutoDiff_Creator** as described below:

2) A ***Multi_AutoDiff_Creator class*** which helps the user create variable objects from a multivariable function

```python

"""
This class helps users initialize different variables from a multivariable function without explicitly specifying them using separate AutoDiff classes

It will need to import AutoDiff Class
"""
	def __init__(self,*args,**kwargs):
		deri = np.identity(len(kwargs))
		i = 0 
		self.Vars =[]
		for key, value in kwargs.items():
			self.List.append(AutoDiff(value,deri[i]))
			i+=1

#Demo and comparison

'''initiate a multivariable function using Multi_AutoDiff_Creator class'''
X,Y= Multi_AutoDiff_Creator(X = 1., Y=3.).Vars
func = X + 2*Y*X 
```

Notice that this class only serves as a helper class to ensure that every variable created has the correct format of dimensionality. The class itself has no influence on the forward mode differentation process.

For better calculation of the derivatives of our elementary functions, we also introduce our elementary function methods.

3)An ***Elementary function*** file which calculate derivatives for elementary functions as described previously
```python  
    #Elementary Function Derivative (not exhaustive):
    def exp(self,ADobj):
    def tan(self,ADobj):
    def sin(self,ADobj):
    def cos(self.ADobj):
# ...

```

4) A ***Jacobian*** class which helps the user compute the Jacobian matrix for vector functions
```python 
class Jacobian:
	def __init__(self,vector_func):
		self.function = vector_func
		self.l = len(self.function)

	def value(self):
		Jacob=[]
		for i in range(self.l):
			Jacob.append(self.function[i].der)
		return np.array(Jacob)


    x, y = ad.Multi_AutoDiff_Creator(x=2, y=3).Vars

    #Deifne a vector function
    func = np.array([x+y, 2*x*y])
    Jacob = ad.Jacobian(func)
    print(Jacob.value()) # this class method output the full jacobian matrix as an np array
```

5) An ***Advanced Feature Class*** which implements more advanced derivatives (such as the Hessian/Backward propagation) which inherits from the AutoDiff class.

6) A ***Test*** class which tests existing module functionalities, which we have already implemented.


**Aspects yet to be implemented:**

1) A dimensionality check method in both the AutoDiff and the Jacobian class to help check the dimensionality of users' operations

2) Error checking for all the trig functions such as sec, tanh, cosh, etc.


3)The advanced feature of back propagation

4) Other error handling issues that might occur



## Future Features

We plan to implement the reverse mode automatic differentiation algorithm for our future advanced feature. The directory structure and basic modules used will remain the same. Additional data structures might involve trees or linked lists if necessary for our implementation of the reverse mode of automatic differentation.

Since reverse mode automatic differentiation will require the storage of all relative derivatives for the final derivative propagation (we will need to keep intermediate data in the memory during the
forward pass for used in the backpropagation), the primary challenge for implementing this new feature is to understand how to store all of the derivatives (in particular what data structure we should use to store them) with respect to their related values and how to recursively calculate the final gradient.

The implementation can be a tree like structure in which, when we are overloading operators, we find a way to store the related values and derivatives as children for the previous nodes and when all the nodes (in this case autodiff objects) are traversed, we use recursion to back propagate the "tree".

Optimally, we would like to implement a very simple neural network finally using gradient descent (or back propagation algorithm) using our self implemented reverse mode with proper layers, nodes, loss and activation function.


d7d9ecd1db88bcd2c2fffebd8fb12c1b7bef235b
ipynb
Jupyter Notebook
docs/milestone2.ipynb
ELAA207/cs207-FinalProject
6e4f265c220a2a156bab51e7717878dea252256c
MIT
2019-10-26T12:46:22.000Z
2019-10-26T12:46:22.000Z
docs/milestone2.ipynb
ELAA207/cs207-FinalProject
6e4f265c220a2a156bab51e7717878dea252256c
MIT
2019-11-19T07:33:03.000Z
2019-12-10T14:53:05.000Z
docs/milestone2.ipynb
ELAA207/cs207-FinalProject
6e4f265c220a2a156bab51e7717878dea252256c
MIT
Qwen/Qwen-72B
1. YES
2. YES