Descenso de gradiente, a mano y después con NumPy

Cuatro líneas y con eso se entrena casi todo 🌸

Ya tienes las dos piezas. En el capítulo 14 viste que el gradiente apunta cuesta arriba. En el capítulo 13 bajaste un valle a mano sin saber que lo estabas haciendo.

Esto es juntarlas y ponerle nombre.

Volvemos a abrir el archivo, escalado

Escalo las columnas desde el principio, y en el capítulo que viene vas a ver por qué no es opcional.

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()

y = X[:, 1]
F = X[:, [0, 2]]
Zf = (F - F.mean(axis=0)) / F.std(axis=0)
A = np.column_stack([np.ones(len(y)), Zf])
n = len(y)
print('A.shape:', A.shape)

A.shape: (2238, 3)

Quiero predecir el monto con las unidades y la satisfacción. Ojo con esto: el monto sale de las columnas que predicen, no está dentro de ellas. Meter la respuesta entre las preguntas es de los errores más caros que hay, y en mi libro de machine learning tiene un capítulo entero.

El algoritmo entero

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

def gradiente(w):
    return -2 * (A.T @ (y - A @ w)) / n

w = np.zeros(3)
for i in range(1, 51):
    w = w - 0.1 * gradiente(w)
    if i in (1, 2, 5, 10, 25, 50):
        print('  iter', str(i).rjust(2), '| error', round(error(w), 2))

  iter  1 | error 980984.52
  iter  2 | error 832210.25
  iter  5 | error 637056.84
  iter 10 | error 575167.98
  iter 25 | error 567732.61
  iter 50 | error 567723.39

Eso es. Una resta dentro de un bucle.

Y mira cómo baja: al principio a saltos enormes, después cada vez menos. Tiene sentido, porque cerca del fondo el gradiente es chiquito, así que el paso también. El algoritmo frena solo al llegar 🐣

¿Llegó al sitio correcto?

optimo = np.linalg.solve(A.T @ A, A.T @ y)
print('descenso :', np.round(w, 4))
print('exacto   :', np.round(optimo, 4))
print('distancia:', round(float(np.linalg.norm(w - optimo)), 6))

descenso : [803.3748   2.679  -16.8303]
exacto   : [803.3863   2.679  -16.8305]
distancia: 0.011469

Cincuenta iteraciones y está a una centésima del óptimo que el capítulo 6 calcula de una sola vez con una fórmula.

Y aquí la pregunta obvia: si la fórmula da el resultado exacto de un tirón, ¿para qué el bucle?

Porque la fórmula solo existe para modelos lineales. Para una red neuronal, para una regresión logística o para casi cualquier otra cosa, no hay fórmula. Solo queda buscar, y buscar es esto 🌟

La tasa de aprendizaje, que es lo único que hay que elegir

for tasa in (0.001, 0.01, 0.1, 0.5, 0.9, 1.01):
    w = np.zeros(3)
    for _ in range(200):
        w = w - tasa * gradiente(w)
    print('  tasa', tasa, '-> error', round(error(w), 2))

  tasa 0.001 -> error 857631.82
  tasa 0.01 -> error 567923.13
  tasa 0.1 -> error 567723.39
  tasa 0.5 -> error 567723.39
  tasa 0.9 -> error 567723.39
  tasa 1.01 -> error 29584831947.19

Lee esa tabla de arriba abajo, porque es el capítulo entero:

• 0.001: doscientas iteraciones y todavía va por 857.000. Llega, pero de aquí a mañana.
• 0.01: casi. Se queda a doscientos del fondo.
• 0.1 a 0.9: perfecto. Hay un rango cómodo y bastante ancho.
• 1.01: veintinueve mil millones. Explotó.

Entre 0.9 y 1.01 no hay casi nada de diferencia, y una funciona y la otra revienta. Por qué existe ese borde tan filoso es el capítulo 18.

Las tres variantes que vas a oír nombrar

Lo que hicimos usa todas las filas en cada paso. Se llama descenso de gradiente por lotes, y con 2.238 filas va sobrado. Con veinte millones no.

• Estocástico, o SGD: usa una fila cada vez. Cada paso es baratísimo y va dando tumbos, porque una sola fila es una opinión muy parcial sobre hacia dónde está el fondo.
• Por minilotes: usa un puñado, típicamente entre 32 y 256. Es el que se usa de verdad en todas partes.
• Con momento, Adam y compañía: le añaden memoria de los pasos anteriores para no rebotar tanto. Adam además le pone tasa propia a cada peso, y eso es lo que intenta arreglar lo que veremos de la curvatura en el capítulo 16.

Todas son la misma línea. Lo único que cambia es con cuántas filas se calcula el gradiente y qué se recuerda entre pasos 🙂

Lo que te llevas

• El descenso de gradiente es restar el gradiente por la tasa, en un bucle.
• Frena solo al acercarse, porque el gradiente se hace chiquito.
• Llega al mismo sitio que la fórmula exacta, pero funciona donde no hay fórmula.
• La tasa es lo único que eliges, y tiene un rango cómodo ancho con un borde filoso.
• SGD, minilotes y Adam son la misma línea con distinta cantidad de filas y distinta memoria.

Comprueba que se entendió

Ejercicios

y = X[:, 1]
F = X[:, [0, 2]]
Zf = (F - F.mean(axis=0)) / F.std(axis=0)
A = np.column_stack([np.ones(len(y)), Zf])
n = len(y)

def gradiente(w):
    return -2 * (A.T @ (y - A @ w)) / n

w = np.zeros(3)
for i in range(1, 41):
    g = gradiente(w)
    w = w - 0.1 * g
    if i in (1, 5, 10, 20, 40):
        print('iter', str(i).rjust(2), '| largo del gradiente',
              round(float(np.linalg.norm(g)), 4))

iter  1 | largo del gradiente 1607.1372
iter  5 | largo del gradiente 658.282
iter 10 | largo del gradiente 215.7054
iter 20 | largo del gradiente 23.1611
iter 40 | largo del gradiente 0.267

y = X[:, 1]
F = X[:, [0, 2]]
Zf = (F - F.mean(axis=0)) / F.std(axis=0)
A = np.column_stack([np.ones(len(y)), Zf])
n = len(y)

def gradiente(w):
    return -2 * (A.T @ (y - A @ w)) / n

for inicio in (np.zeros(3), np.array([5000.0, -300.0, 900.0])):
    w = inicio.copy()
    for _ in range(300):
        w = w - 0.1 * gradiente(w)
    print('desde', np.round(inicio, 1), '-> llega a', np.round(w, 4))

desde [0. 0. 0.] -> llega a [803.3863   2.679  -16.8305]
desde [5000. -300.  900.] -> llega a [803.3863   2.679  -16.8305]

y = X[:, 1]
F = X[:, [0, 2]]
Zf = (F - F.mean(axis=0)) / F.std(axis=0)
A = np.column_stack([np.ones(len(y)), Zf])
n = len(y)

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

rng = np.random.default_rng(0)
w = np.zeros(3)
for paso in range(1, 4001):
    i = int(rng.integers(n))
    g = -2 * A[i] * (y[i] - A[i] @ w)
    w = w - 0.01 * g
    if paso in (100, 500, 1000, 4000):
        print('paso', str(paso).rjust(4), '| error', round(error(w), 2))

paso  100 | error 593507.79
paso  500 | error 584273.73
paso 1000 | error 596225.78
paso 4000 | error 585674.45