# Diferenciação automática

## Formas mais convencionais de calcular o gradiente

- Manual
- Simbólica
- Numérica

## Forma manual

Cálculo segundo regras de derivação de cada uma das funções e operações envolvidas. Exemplo:

\[y = \frac{x+3}{x+4}\]

Usando a regra do quociente

\[\left(\frac{f}{g}\right)' = \frac{f' \cdot g - f \cdot g'}{g^2},\]

obtemos

\[y' = \frac{1 \cdot (x+4) - (x+3) \cdot 1}{(x+4)^2} = \frac{1}{(x+4)^2}.\]

## Forma simbólica

- Da mesma forma que a manual, mas usando o computador

- No Python, pode ser usado o SymPy: https://www.sympy.org/

    from sympy import symbols, diff

    x, y  = symbols('x, y')
    y = (x + 3) / (x + 4)

    y_p = diff(y, x)
    y_p

\(\displaystyle - \frac{x + 3}{\left(x + 4\right)^{2}} + \frac{1}{x + 4}\)

Vale notar que a expressão pode ser simplificada para \(\frac{1}{(x+4)^2}\)

## Forma simbólica

- Avaliando em \(x = 6\)

    x_v = 6
    values = {x:x_v}
    print("Valor de y")
    y.subs(values)

    Valor de y

\(\displaystyle \frac{9}{10}\)

    print("Valor de y_p")
    y_p.subs(values)

    Valor de y_p

\(\displaystyle \frac{1}{100}\)

## Forma numérica

\[y' = f'(x) = {\rm lim}_{\Delta \rightarrow 0}\frac{f(x + \Delta) - f(x)}{\Delta}\]

    x = 6

    # f(x)
    y = (x + 3) / (x + 4)

    Delta = 1e-8

    # f(x +  Delta)
    y_Delta = ((x + Delta) + 3) / ((x + Delta) + 4)

    y_p = (y_Delta - y)/Delta
    y_p

    0.01000000082740371

## Forma numérica

\[y' = f'(x) = {\rm lim}_{\Delta \rightarrow 0}\frac{f(x + \Delta) - f(x)}{\Delta}\]

    x = 6

    # f(x)
    y = (x + 3) / (x + 4)

    # Delta = 1e-8
    Delta = 1e-10

    # f(x +  Delta)
    y_Delta = ((x + Delta) + 3) / ((x + Delta) + 4)

    y_p = (y_Delta - y)/Delta
    y_p

    0.009999778782798785

## Forma numérica

\[y' = f'(x) = {\rm lim}_{\Delta \rightarrow 0}\frac{f(x + \Delta) - f(x)}{\Delta}\]

    x = 6

    # f(x)
    y = (x + 3) / (x + 4)

    # Delta = 1e-8
    # Delta = 1e-10
    Delta = 1e-15

    # f(x +  Delta)
    y_Delta = ((x + Delta) + 3) / ((x + Delta) + 4)

    y_p = (y_Delta - y)/Delta
    y_p

    0.0

## Diferenciação automática

- Nome mais adequado: diferenciação algorítmica

- Cálculo exato

- A ideia é definir formas das derivadas de operações “atômicas”

- É a técnica usada em frameworks como PyTorch e Tensorflow

- Duas formas mais conhecidas:
  - Progressiva
  - Regressiva

## Problema da forma progressiva

\(L\) é calculado para um batch de dados \(\left[ \begin{matrix} x_1 \\ x_2 \\ \vdots \\ x_N \end{matrix}\right]\)

\(L = L_1 + L_2 + \cdots + L_N\), \(\;\;\;\;\; L_i = f_i(\{w_k\}, x_i)\)

Podemos dizer que a loss é função do vetor de losses:

\(L = f\left(\left[\begin{matrix} L_1 \\ L_2 \\ \vdots \\ L_N \end{matrix}\right]\right) = f(\mathbf{L})\), \(\;\;\;\; \mathbf{L} = \left[\begin{matrix} L_1 \\ L_2 \\ \vdots \\ L_N \end{matrix}\right] = \left[\begin{matrix} f_1(\{w_k\}, x_1) \\ f_2(\{w_k\}, x_2) \\ \vdots \\ f_N(\{w_k\}, x_N) \end{matrix}\right]\)

## Problema da forma progressiva

Estamos interessados:

Na função vetorial \(\mathbf{L}\):

\(\mathbf{L} = \left[\begin{matrix} L_1 \\ L_2 \\ \vdots \\ L_N \end{matrix}\right] = \left[\begin{matrix} f_1(\{w_k\}, x_1) \\ f_2(\{w_k\}, x_2) \\ \vdots \\ f_N(\{w_k\}, x_N) \end{matrix}\right]\)

E no Jacobiano de \(\mathbf{L}\):

\(J_{\mathbf{L}} = \left[\begin{matrix} \frac{\partial L_1}{\partial w_1} & \cdots & \frac{\partial L_1}{\partial w_K} \\ \vdots & \ddots & \vdots \\ \frac{\partial L_N}{\partial w_1} & \cdots & \frac{\partial L_N}{\partial w_K} \end{matrix}\right]\)

O método progressivo calcula uma coluna de \(J_{\mathbf{L}}\) por vez

O problema é que \(K >> N \Rightarrow\) Há muito mais pesos do que elementos no batch

O método regressivo se propõe a calcular uma linha de cada vez

- Necessita armazenar informação do grafo na construção do tipo

## Implementação da diferenciação automática - forma progressiva

\(\;\;\)

\(v_1 = x = 6\)

\(\dot{v}_1 \triangleq \partial v_1 / \partial x\)

\(\dot{v}_1 = 1\)

## Implementação da diferenciação automática - forma progressiva

\(\;\;\)

\(v_1 = x = 6\)

\(v_2 = 3\)

\(\dot{v}_1 \triangleq \partial v_1 / \partial x\)

\(\dot{v}_1 = 1\)

\(\dot{v}_2 = 0\)

Vamos para a constante do numerador agora. Então tem o 3 aqui, que aparece no numerador. Ele vai ser o nosso v_2. Então v_2 é 3 e a derivada de v_2 em relação à v_1, o v_2 ponto é zero, porque v_2 é constante.

## Implementação da diferenciação automática - forma progressiva

\(\;\;\)

\(v_1 = x = 6\)

\(v_2 = 3\)

\(v_3 = v_1 + v_2 = 9\)

\(\dot{v}_1 \triangleq \partial v_1 / \partial x\)

\(\dot{v}_1 = 1\)

\(\dot{v}_2 = 0\)

\(\dot{v}_3 = \dot{v}_1 + \dot{v}_2 = 1\)

## Implementação da diferenciação automática - forma progressiva

\(\;\;\)

\(v_1 = x = 6\)

\(v_2 = 3\)

\(v_3 = v_1 + v_2 = 9\)

\(v_4 = 4\)

\(\dot{v}_1 \triangleq \partial v_1 / \partial x\)

\(\dot{v}_1 = 1\)

\(\dot{v}_2 = 0\)

\(\dot{v}_3 = \dot{v}_1 + \dot{v}_2 = 1\)

\(\dot{v}_4 = 0\)

## Implementação da diferenciação automática - forma progressiva

\(\;\;\)

\(v_1 = x = 6\)

\(v_2 = 3\)

\(v_3 = v_1 + v_2 = 9\)

\(v_4 = 4\)

\(v_5 = v_1 + v_4 = 10\)

\(\dot{v}_1 \triangleq \partial v_1 / \partial x\)

\(\dot{v}_1 = 1\)

\(\dot{v}_2 = 0\)

\(\dot{v}_3 = \dot{v}_1 + \dot{v}_2 = 1\)

\(\dot{v}_4 = 0\)

\(\dot{v}_5 = \dot{v}_1 + \dot{v}_4 = 1\)

## Implementação da diferenciação automática - forma progressiva

\(\;\;\)

\(v_1 = x = 6\)

\(v_2 = 3\)

\(v_3 = v_1 + v_2 = 9\)

\(v_4 = 4\)

\(v_5 = v_1 + v_4 = 10\)

\(\displaystyle v_6 = \frac{v_3}{v_5} = \frac{9}{10} = y\)

\(\dot{v}_1 \triangleq \partial v_1 / \partial x\)

\(\dot{v}_1 = 1\)

\(\dot{v}_2 = 0\)

\(\dot{v}_3 = \dot{v}_1 + \dot{v}_2 = 1\)

\(\dot{v}_4 = 0\)

\(\dot{v}_5 = \dot{v}_1 + \dot{v}_4 = 1\)

\(\dot{v}_6 = \frac{\dot{v}_3 v_5 - v_3 \dot{v}_5}{v_5^2} = 0,01\)

## Sobrecarga de operadores

Tipo para representar uma certa quantidade de minutos e segundos:

    class Tempo():
        def __init__(self, m, s):
            self.m = m
            self.s = s
        def __repr__(self):
            return f"{self.m} minutos e {self.s} segundos"

    t1 = Tempo(2, 20)
    print(t1)

    2 minutos e 20 segundos

    t2 = Tempo(3, 45)
    print(t2)

    3 minutos e 45 segundos

    #t1 + t2 leva a exceção

    ---------------------------------------------------------------------------
    TypeError                                 Traceback (most recent call last)
    Cell In[9], line 1
    ----> 1 t1 + t2

    TypeError: unsupported operand type(s) for +: 'Tempo' and 'Tempo'

## Sobrecarga de operadores

É possível programar o comportamento dos operadores (ex. `+`, `-`, etc) quando utilizados com tipos customizados como o tipo `Tempo` que definimos. A esse processo dá-se no nome de sobrecarga de operadores.

    class Tempo():
        def __init__(self, m, s):
            self.m = m
            self.s = s
        def __repr__(self):
            return f"{self.m} minutos e {self.s} segundos"
        def __add__(self, other):
            ms = self.m + other.m
            ss = self.s + other.s
            if ss > 60:
                ms += 1
                ss -= 60
            return Tempo(ms, ss)

    t1 = Tempo(2, 20)
    t2 = Tempo(3, 45)
    soma = t1 + t2
    print(soma)

    6 minutos e 5 segundos

O Python permite fazer sobrecarga em diversos operadores: https://docs.python.org/3/reference/datamodel.html#emulating-numeric-types

## Sobrecarga de operadores

Cuidados:

- Operando com um tipo de dado não suportado:

    t1 = Tempo(2, 20)
    t2 = 3

    #t1 + t2 # Leva à exceção

    ---------------------------------------------------------------------------
    AttributeError                            Traceback (most recent call last)
    Cell In[16], line 1
    ----> 1 t1 + t2 # Leva à exceção

    Cell In[11], line 9, in Tempo.__add__(self, other)
          8 def __add__(self, other):
    ----> 9     ms = self.m + other.m
         10     ss = self.s + other.s
         11     if ss > 60:

    AttributeError: 'int' object has no attribute 'm'

## Sobrecarga de operadores

- É necessário converter os tipos adequadamente:

    class Tempo():
        def __init__(self, m, s):
            self.m = m
            self.s = s
        def __repr__(self):
            return f"{self.m} minutos e {self.s} segundos"
        def __add__(self, other):
            if not isinstance(other, Tempo):
                other = Tempo(other, 0)
            ms = self.m + other.m
            ss = self.s + other.s
            if ss > 60:
                ms += 1
                ss -= 60
            return Tempo(ms, ss)

    t1 = Tempo(2, 20)
    t2 = 3

    t1 + t2

    5 minutos e 20 segundos

## Sobrecarga de operadores

- Operações com tipos simples à esquerda:

    t1 = 3
    t2 = Tempo(2, 20)

    #t1 + t2 # Leva à exceção

    ---------------------------------------------------------------------------
    TypeError                                 Traceback (most recent call last)
    Cell In[23], line 1
    ----> 1 t1 + t2 # Leva à exceção

    TypeError: unsupported operand type(s) for +: 'int' and 'Tempo'

## Sobrecarga de operadores

- Existem métodos para especificar o comportamento com operações à esquerda:

    class Tempo():
        def __init__(self, m, s):
            self.m = m
            self.s = s
        def __repr__(self):
            return f"{self.m} minutos e {self.s} segundos"
        def __add__(self, other):
            if not isinstance(other, Tempo):
                other = Tempo(other, 0)
            ms = self.m + other.m
            ss = self.s + other.s
            if ss > 60:
                ms += 1
                ss -= 60
            return Tempo(ms, ss)
        def __radd__(self, lhs):
            return self + lhs

    t1 = Tempo(2, 20)
    t2 = 3

    t1 + t2

    5 minutos e 20 segundos

## Implementação do Tipo Dual

`y = (x + 3) / (x + 4)`

    class D:
        def __init__(self, x, xd):
            self.x = x
            self.xd = xd

        def __add__(self, other):
            # Self é sempre do tipo D
            # other pode ser uma constante
            # Derivada de constante é zero
            if not isinstance(other, D):
                other = D(other, 0)
            return D(self.x + other.x, self.xd + other.xd)

        def __truediv__(self, other):
            # Self é sempre do tipo D
            # other pode ser uma constante
            # Derivada de constante é zero
            if not isinstance(other, D):
                other = D(other, 0)
            return D(self.x / other.x, (self.xd * other.x - self.x * other.xd) / (other.x**2))

## Implementação do Tipo Dual 2

`y = (x + 3) / (x + 4)`

    class D:
        def __init__(self, x, xd):
            self.x = x
            self.xd = xd

        def __add__(self, other):
            # Self é sempre do tipo D
            # other pode ser uma constante
            # Derivada de constante é zero
            if not isinstance(other, D):
                other = D(other, 0)
            return D(self.x + other.x, self.xd + other.xd)

        def __radd__(self, lhs):
            # Se o método de soma do tipo à esquerda falha,
            # o Python chama o método __radd__ do tipo à direita
            return self + lhs

        def __truediv__(self, other):
            # Self é sempre do tipo D
            # other pode ser uma constante
            # Derivada de constante é zero
            if not isinstance(other, D):
                other = D(other, 0)
            return D(self.x / other.x, (self.xd * other.x - self.x * other.xd) / (other.x**2))

        def __rtruediv__(self, lhs):
            # Se o método de divisão do tipo à esquerda falha,
            # o Python chama o método __rtruediv__ do tipo à direita
            if not isinstance(lhs, D):
                lhs = D(lhs, 0)
            return lhs / self

        def __repr__(self):
            return f"{self.x} + {self.xd}ϵ"

## Tipo Dual - Teste

`y = (x + 3) / (x + 4)`

    x = 6
    x_D = D(x, 1)

    y_D = (x_D + 3) / (x_D + 4)

    print(y_D)

    0.9 + 0.01ϵ

Outra forma:

    def f(x):
        return (x + 3) / (x + 4)

    def diff(f, x):
        return f(D(x, 1)).xd

    x = 6
    f(x)

    0.9

    diff(f, x)

    0.01

## Tipo Dual - funções complicadas

- Raiz quadrada: \(\displaystyle y_2=\sqrt{x} \Rightarrow \frac{dy_2}{dx} = \frac{1}{2\sqrt{x}}\)

    import math
    x = 6
    y_2 = math.sqrt(x)
    print(y_2)

    2.449489742783178

    y_2_p = 1/(2*math.sqrt(x))
    print(y_2_p)

    0.20412414523193154

    x_D = D(x, 1)
    #y_2_D = math.sqrt(x_D)

Leva à exceção:

    TypeError                                 Traceback (most recent call last)
    Cell In[54], line 4
          1 x_D = D(x, 1)
    ----> 2 y_2_D = math.sqrt(x_D)

    TypeError: must be real number, not D

## Tipo Dual - funções complicadas

- Método babilônio (ref.: https://blogs.sas.com/content/iml/2016/05/16/babylonian-square-roots.html):

    def babylonian(x, N = 10):
        t = (1 + x) / 2
        for i in range(2, N):
            t = (t + x / t) / 2
        return t

- Usa apenas somas e divisões

    y_2_b = babylonian(x)
    print(y_2_b)

    2.449489742783178

    print(y_2) # calculado com a função math.sqrt

    2.449489742783178

    x_D = D(x, 1)
    y_2_D = babylonian(x_D)
    print(y_2_D)

    2.449489742783178 + 0.2041241452319315ϵ

    print(y_2_p) # calculado com 1/(2*math.sqrt(x))

    0.20412414523193154

## Diferenciação progressiva - aplicação em redes neurais

- Função custo, por exemplo: \(\displaystyle L = \sum_i e_i^2\)

- Queremos calcular \(\displaystyle \frac{\partial L}{\partial w_k}\), para todos os \(w_k\)

- Para simplificar:
  - Rede com dois parâmetros \(w_1\) e \(w_2\)
  - Função custo (só com as operações que implementamos) \(\displaystyle L = e_1 + e_2\), em que:
    - \(x_i = 1\) (entrada constante), \(\;\;\;\) \(y_{1} = \frac{w_1}{w_2}\cdot x_1 = \frac{w_1}{w_2}\)
    - \(y_{2} = \frac{w_1}{w_2}\cdot x_2 + w_1 = \frac{w_1}{w_2} + w_1\)
    - \(e_i = \frac{d_i}{y_i}\)

## Diferenciação progressiva - aplicação em redes neurais

    def L_func(w1, w2, d1, d2):
        y1 = w1/w2
        y2 = w1/w2 + w1
        e1 = d1/y1
        e2 = d2/y2
        L = e1 + e2
        L_aux = (e1, e2)
        return L, L_aux

Supondo valores para \(w_1\), \(w_2\), \(d_1\) e \(d_2\), calculamos o valor da função:

    w1 = 1
    w2 = 2
    d1 = 3
    d2 = 2
    L, L_aux = L_func(w1, w2, d1, d2)
    print(L)

    7.333333333333333

## Diferenciação progressiva - aplicação em redes neurais

Para calcular o valor da derivada em relação à \(w_1\)

    w1 = D(1, 1) # derivada igual a 1
    w2 = D(2, 0) # constante -> derivada igual a zero
    L, L_aux = L_func(w1, w2, d1, d2)
    print(L)

    7.333333333333333 + -7.333333333333333ϵ

    w1 = D(1, 0) # constante -> derivada igual a zero
    w2 = D(2, 1) # derivada igual a 1
    L, L_aux = L_func(w1, w2, d1, d2)
    print(L)

    7.333333333333333 + 3.2222222222222223ϵ

    print(L_aux)

    (6.0 + 3.0ϵ, 1.3333333333333333 + 0.2222222222222222ϵ)

- De uma só vez obtemos:
  - Derivada em relação a um parâmetro, para vários exemplos de um batch (1 e 2 nesse caso)
    - \(\displaystyle \frac{\partial L}{\partial w_2} = \frac{\partial e_1}{\partial w_2} + \frac{\partial e_2}{\partial w_2}\)