Jacobiana y hessiana, lo justo

Este capítulo es corto a propósito 🌸

Jacobiana y hessiana tienen fama de tema avanzado, y para lo que tú necesitas son dos ideas simples. Te las doy y seguimos.

Volvemos a abrir el archivo

Cada capítulo arranca de cero, así que estas líneas se repiten. Son las mismas del capítulo 1 y las vas a ver en todos.

import pandas as pd
import numpy as np

df = pd.read_csv('ventas-miss-yera.csv')
for c in ('unidades', 'monto', 'satisfaccion'):
    df[c] = pd.to_numeric(df[c], errors='coerce')
X = df[['unidades', 'monto', 'satisfaccion']].dropna().to_numpy()
print('X.shape:', X.shape)

X.shape: (2238, 3)

La jacobiana: cuando la función devuelve varios números

El gradiente del capítulo 14 es de una función que devuelve un número, el error. Si la función devuelve varios, hace falta una derivada por cada entrada y cada salida. Eso es una matriz, y se llama jacobiana.

s = X.std(axis=0)
J = np.diag(1 / s)
print('desviaciones:', np.round(s, 4))
print('jacobiana del escalado:')
print(np.round(J, 6))

desviaciones: [  5.7856 753.6678   1.424 ]
jacobiana del escalado:
[[0.172842 0.       0.      ]
 [0.       0.001327 0.      ]
 [0.       0.       0.702235]]

Esa es la jacobiana de escalar tus columnas: entran tres números y salen tres. Sale diagonal porque cada columna se escala sola, sin mirar a las otras.

Y hay un caso que ya conoces: la jacobiana de multiplicar por una matriz es esa misma matriz. Si tu función es Wx, su jacobiana es W. Por eso en una red los pesos aparecen tal cual en la retropropagación del capítulo 15: la cadena va multiplicando jacobianas, y las jacobianas de las capas lineales son los propios pesos 🐣

La hessiana: la curvatura del valle

La derivada dice la pendiente. La derivada segunda dice cuánto cambia esa pendiente, o sea si el valle es cerrado y empinado o abierto y plano.

y = X[:, 1]
A = np.column_stack([np.ones(len(y)), X[:, 0]])
H = 2 * (A.T @ A) / len(y)
print('hessiana:')
print(np.round(H, 4))

hessiana:
[[  2.      23.3521]
 [ 23.3521 339.6077]]

Es simétrica, como X^⊤X en el capítulo 6, y por el mismo motivo. Y eso quiere decir que le podemos sacar autovalores con eigh, como en el capítulo 10.

Sus autovalores dicen dos cosas

vals = np.linalg.eigvalsh(H)
print('autovalores        :', np.round(vals, 4))
print('todos positivos    :', bool(np.all(vals > 0)))
print('numero de condicion:', round(float(vals.max() / vals.min()), 2))

autovalores        : [  0.3924 341.2153]
todos positivos    : True
numero de condicion: 869.55

Todos positivos quiere decir que en todas las direcciones la superficie sube. O sea que ahí hay un fondo de verdad y no una silla de montar. Cuando eso pasa en todo el espacio, el problema se llama convexo, y convexo quiere decir que hay un solo mínimo y no te puedes quedar atrapada.

El cociente entre el mayor y el menor es el número de condición, el mismo del capítulo 8. Aquí vale 869, y traducido quiere decir que el valle es 869 veces más empinado en una dirección que en la otra.

Imagínate un cañón largo y estrecho. Si caminas siempre en la dirección de máxima pendiente, rebotas de una pared a la otra y avanzas poquito hacia el fondo. Eso es el zigzag del descenso de gradiente, y ese 869 es el número que lo predice antes de que pase 🙃

Y por eso existe el método de Newton

w = np.array([500.0, 30.0])
g = -2 * (A.T @ (y - A @ w)) / len(y)
paso = np.linalg.solve(H, g)
nuevo = w - paso
optimo = np.linalg.solve(A.T @ A, A.T @ y)

print('gradiente    :', np.round(g, 4))
print('un paso da   :', round(float(nuevo[0]), 4), round(float(nuevo[1]), 4))
print('el optimo es :', round(float(optimo[0]), 4), round(float(optimo[1]), 4))

gradiente    : [  93.7905 3069.8176]
un paso da   : 797.5076 0.5035
el optimo es : 797.5076 0.5035

Un paso. No cien, ni mil: uno, y clavado.

Newton corrige el gradiente con la curvatura, así que sabe estirar el paso donde el valle es plano y acortarlo donde es empinado. Sin zigzag.

Y entonces la pregunta obvia: si esto es tan bueno, ¿por qué las redes se entrenan con descenso de gradiente y no con Newton?

Por el tamaño. La hessiana de un modelo con un millón de pesos tiene un billón de casillas, y además hay que resolver un sistema con ella en cada paso. No cabe ni en la memoria ni en el presupuesto.

Por eso existen los métodos que están en medio, como L-BFGS, que se inventan una hessiana aproximada, o como Adam, que guarda un resumen barato de la curvatura por cada peso. Todos intentan lo mismo: algo de Newton sin pagar Newton 🌟

Lo que te llevas

• La jacobiana es la matriz de derivadas primeras cuando salen varios números.
• La jacobiana de multiplicar por W es W.
• La hessiana son las derivadas segundas y mide la curvatura.
• Autovalores todos positivos: hay un mínimo de verdad y el problema es convexo.
• El número de condición dice cuánto vas a zigzaguear.
• Newton llega en un paso y no se usa porque la hessiana no cabe.

Comprueba que se entendió

Ejercicios

y = X[:, 1]
A = np.column_stack([np.ones(len(y)), X[:, 0]])

def error(w):
    return float(((y - A @ w) ** 2).mean())

w = np.array([500.0, 30.0])
h = 1e-4
Hn = np.zeros((2, 2))
for i in range(2):
    for j in range(2):
        pp, pm, mp, mm = w.copy(), w.copy(), w.copy(), w.copy()
        pp[i] += h; pp[j] += h
        pm[i] += h; pm[j] -= h
        mp[i] -= h; mp[j] += h
        mm[i] -= h; mm[j] -= h
        Hn[i, j] = (error(pp) - error(pm) - error(mp) + error(mm)) / (4 * h * h)

print(np.round(Hn, 4))
print('coincide con la exacta:',
      np.allclose(Hn, 2 * (A.T @ A) / len(y), rtol=1e-3))

[[  1.9994  23.35  ]
 [ 23.35   339.6097]]
coincide con la exacta: True

y = X[:, 1]
u = X[:, 0]
A = np.column_stack([np.ones(len(y)), u])
uz = (u - u.mean()) / u.std()
Az = np.column_stack([np.ones(len(y)), uz])

for nombre, M in (('sin escalar', A), ('escalado  ', Az)):
    H = 2 * (M.T @ M) / len(y)
    v = np.linalg.eigvalsh(H)
    print(nombre, '-> autovalores', np.round(v, 4),
          '| condicion', round(float(v.max() / v.min()), 2))

sin escalar -> autovalores [  0.3924 341.2153] | condicion 869.55
escalado   -> autovalores [2. 2.] | condicion 1.0

Silla = np.array([[2.0, 0.0], [0.0, -2.0]])
print('autovalores:', np.linalg.eigvalsh(Silla))
print('hay minimo :', bool(np.all(np.linalg.eigvalsh(Silla) > 0)))

autovalores: [-2.  2.]
hay minimo : False