Cómo sabe la red hacia dónde corregir

En el capítulo 4 la red aprendió y yo te dije que no miraras las líneas del "hacia atrás". Aquí las miramos 🔍

Y de paso pago lo que quedé debiendo en el capítulo 1, cuando dije que la razón de escalar estaba en el gradiente. Está, y se ve clarísima.

Primero: cuánto te equivocas

Para corregir hace falta un número que diga cómo de mal vas. Ese número es la pérdida, y para un sí o no se usa esta:

import numpy as np

def sigmoide(z):
    return 1 / (1 + np.exp(-z))

def perdida(p, y):
    p = np.clip(p, 1e-12, 1 - 1e-12)      # para que log(0) no reviente
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))

y = np.array([1., 1., 0., 0.])
print('acertando de lleno  :', round(perdida(np.array([.99, .99, .01, .01]), y), 4))
print('sin idea (0,5 a todo):', round(perdida(np.array([.5, .5, .5, .5]), y), 4))
print('equivocada de lleno :', round(perdida(np.array([.01, .01, .99, .99]), y), 4))

acertando de lleno  : 0.0101
sin idea (0,5 a todo): 0.6931
equivocada de lleno : 4.6052

Se llama entropía cruzada y funciona así: castiga estar seguro y equivocado mucho más que estar dudando.

Fíjate en los números. Dudar cuesta 0,6931, que ya salió en el capítulo 4 y es ln(2). Y estar seguro y equivocado cuesta 4,6052, o sea casi siete veces más.

Esa asimetría es a propósito. Un modelo que dice "70%" y falla es un modelo que se puede arreglar; uno que dice "99%" y falla está roto 😬

Segundo: hacia dónde crece el error

El gradiente es una lista de números, uno por peso, que dice cuánto sube la pérdida si subes ese peso un poquito.

Para una neurona con sigmoide y entropía cruzada, sale una fórmula sorprendentemente corta:

rng = np.random.default_rng(0)
X = rng.normal(0, 1, (200, 3))
y = (rng.random(200) < sigmoide(X @ np.array([1., -0.5, 0.3]))).astype(float)

w = np.array([0.2, -0.1, 0.4])
b = 0.05

error = sigmoide(X @ w + b) - y
gradiente = X.T @ error / len(y)
print('gradiente:', np.round(gradiente, 6))

gradiente: [-0.139486  0.142932  0.048427]

Esa línea, X.T @ (predicho - real) / n, es todo el cálculo. Sale de derivar la pérdida y se simplifica muchísimo porque la sigmoide y la entropía cruzada están hechas la una para la otra.

Ahora, ¿cómo sabemos que esa fórmula está bien? Se comprueba moviendo el peso a mano y midiendo cuánto cambia la pérdida:

def perdida_de(w, b):
    return perdida(sigmoide(X @ w + b), y)

h = 1e-5
numerico = np.zeros(3)
for i in range(3):
    arriba = w.copy(); arriba[i] += h
    abajo = w.copy(); abajo[i] -= h
    numerico[i] = (perdida_de(arriba, b) - perdida_de(abajo, b)) / (2 * h)

print('por fórmula:', np.round(gradiente, 8))
print('moviendo    :', np.round(numerico, 8))
print('diferencia  :', float(np.abs(gradiente - numerico).max()))

por fórmula: [-0.1394857   0.14293227  0.0484271 ]
moviendo    : [-0.1394857   0.14293227  0.0484271 ]
diferencia  : 6.669387264679472e-12

6,7e-12 de diferencia, o sea que la fórmula está bien 🎯

Esto tiene nombre, gradient checking, y es lo que se hace cuando escribes una capa nueva a mano y no estás segura de la derivada. Es lento (hay que evaluar la pérdida dos veces por peso) así que solo se usa para comprobar, nunca para entrenar.

Te lo enseño porque es la herramienta que convierte "creo que la derivada está bien" en "está bien", y esa diferencia vale mucho.

Tercero: dar el paso

El gradiente apunta hacia donde el error crece. Así que para mejorar hay que ir al revés:

print('pérdida ahora     :', round(perdida_de(w, b), 6))
print('un paso a favor   :', round(perdida_de(w - 0.5 * gradiente, b), 6))
print('un paso al revés  :', round(perdida_de(w + 0.5 * gradiente, b), 6))

pérdida ahora     : 0.648283
un paso a favor   : 0.628465
un paso al revés  : 0.67071

Restando baja de 0,6483 a 0,6285 y sumando sube a 0,6707. Por eso todas las líneas de corrección del libro llevan un menos 🔽

Y ese 0,5 que multiplica es el paso, o tasa de aprendizaje, que en el capítulo 4 ya vimos que puede arruinarlo todo por los dos lados.

Entrenar es repetir esas tres cosas hasta que deje de mejorar. Nada más.

Lo que quedó pendiente en el capítulo 1

Ahora sí. Armo un problema con dos columnas de escalas muy distintas, que es lo que pasa siempre en datos de verdad: un monto en cientos de soles y una satisfacción del 1 al 5.

n = 500
monto = rng.normal(800, 300, n)
satisfaccion = rng.normal(3.5, 1.0, n)
z = 0.002 * (monto - 800) + 0.8 * (satisfaccion - 3.5)
objetivo = (rng.random(n) < sigmoide(z)).astype(float)

crudo = np.column_stack([monto, satisfaccion])
escalado = (crudo - crudo.mean(axis=0)) / crudo.std(axis=0)

def baja(datos, paso, vueltas=200):
    w = np.zeros(datos.shape[1])
    b = 0.0
    inicial = perdida(sigmoide(datos @ w + b), objetivo)
    for _ in range(vueltas):
        e = sigmoide(datos @ w + b) - objetivo
        w -= paso * (datos.T @ e) / len(objetivo)
        b -= paso * e.mean()
    return inicial, perdida(sigmoide(datos @ w + b), objetivo)

for nombre, datos in [('sin escalar', crudo), ('escalado', escalado)]:
    for paso in [0.00001, 0.0001, 0.01, 0.1]:
        ini, fin = baja(datos, paso)
        estado = 'DIVERGE' if fin > ini else f'{fin:.4f}'
        print(f'{nombre:12} paso {paso:<8} {estado}')

sin escalar  paso 1e-05    0.6806
sin escalar  paso 0.0001   DIVERGE
sin escalar  paso 0.01     DIVERGE
sin escalar  paso 0.1      DIVERGE
escalado     paso 1e-05    0.6930
escalado     paso 0.0001   0.6920
escalado     paso 0.01     0.6175
escalado     paso 0.1      0.5544

Ahí está la respuesta 😳

Sin escalar, desde 0,0001 la pérdida explota. El único paso que sobrevive es 0,00001, y con ese la pérdida baja de 0,6931 a 0,6806 en 200 vueltas.

Con las columnas escaladas, ese mismo paso de 0,00001 apenas mueve nada (0,6930), pero como ahí sí puedes usar 0,1, la pérdida llega a 0,5544.

Compara el avance: sin escalar ganó 0,0125 de pérdida y escalado ganó 0,1387, o sea once veces más, en las mismas 200 vueltas.

El mecanismo es este. El gradiente lleva la escala de la columna dentro: la del monto sale unas trescientas veces más grande que la de satisfacción. Un paso suficientemente chico para no hacer explotar el peso del monto es trescientas veces demasiado chico para mover el de satisfacción. Y no existe un número que sirva para los dos.

Escalar es poner las dos columnas en la misma escala para que un solo paso les valga a las dos. Por eso no es una buena práctica: es la condición para que el descenso de gradiente funcione ⚖️

Y por eso también el StandardScaler va dentro del pipeline desde el capítulo 1.

Y ahora sobre las ventas de verdad

Todo lo anterior fue con datos que fabriqué yo, para poder comparar contra la respuesta verdadera. Vamos a soltar el mismo bucle sobre el CSV de la distribuidora 🐔

import pandas as pd
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

URL = 'https://missyera.com/static/datasets/ventas-miss-yera.csv'

def carga_limpia(url):
    v = pd.read_csv(url).drop_duplicates()
    v['ciudad'] = (v['ciudad'].str.strip().str.lower()
                   .str.normalize('NFKD')
                   .str.encode('ascii', 'ignore').str.decode('utf-8'))
    v['monto'] = pd.to_numeric(v['monto'].str.replace(',', '.'))
    for col in ['fecha', 'fecha_ultima_compra']:
        f = pd.to_datetime(v[col], format='%Y-%m-%d', errors='coerce')
        falta = f.isna() & v[col].notna()
        f[falta] = pd.to_datetime(v.loc[falta, col], format='%d/%m/%Y', errors='coerce')
        v[col] = f
    return v

def prepara(v):
    v = v.sort_values(['cliente_id', 'fecha']).copy()
    v['sin_compra_previa'] = v['fecha_ultima_compra'].isna().astype(int)
    v['sin_descuento'] = v['descuento'].isna().astype(int)
    v['sin_satisfaccion'] = v['satisfaccion'].isna().astype(int)
    v['precio_unitario'] = v['monto'] / v['unidades']
    v['visita_numero'] = v.groupby('cliente_id').cumcount() + 1
    return v

NUMERICAS = ['unidades', 'monto', 'descuento', 'satisfaccion', 'precio_unitario',
             'sin_compra_previa', 'sin_descuento', 'sin_satisfaccion', 'visita_numero']
CATEGORICAS = ['ciudad', 'segmento', 'canal', 'categoria']

ventas = prepara(carga_limpia(URL))
Xv = ventas[NUMERICAS + CATEGORICAS]
yv = ventas['compro']
Xv_tr, Xv_te, yv_tr, yv_te = train_test_split(Xv, yv, test_size=0.25,
                                              random_state=42, stratify=yv)

sklearn_modelo = Pipeline([
    ('pre', ColumnTransformer([
        ('num', Pipeline([('r', SimpleImputer(strategy='median')),
                          ('e', StandardScaler())]), NUMERICAS),
        ('cat', Pipeline([('r', SimpleImputer(strategy='most_frequent')),
                          ('c', OneHotEncoder(handle_unknown='ignore'))]), CATEGORICAS),
    ])),
    ('mod', LogisticRegression(max_iter=1000, random_state=42)),
]).fit(Xv_tr, yv_tr)

T_tr = sklearn_modelo.named_steps['pre'].transform(Xv_tr)
T_te = sklearn_modelo.named_steps['pre'].transform(Xv_te)
objetivo_v = yv_tr.values.astype(float)

w_mano = np.zeros(T_tr.shape[1])
b_mano = 0.0
for _ in range(3000):
    e = sigmoide(T_tr @ w_mano + b_mano) - objetivo_v
    w_mano -= 0.5 * (T_tr.T @ e) / len(objetivo_v)
    b_mano -= 0.5 * e.mean()

print('nuestro bucle:', round(roc_auc_score(yv_te, sigmoide(T_te @ w_mano + b_mano)), 4))
print('scikit-learn :', round(roc_auc_score(
    yv_te, sklearn_modelo.predict_proba(Xv_te)[:, 1]), 4))

nuestro bucle: 0.7216
scikit-learn : 0.7214

0,7216 el bucle que acabamos de escribir, contra 0,7214 la librería. Con veintiocho columnas de ventas peruanas de verdad, con su suciedad, sus nulos y su OneHotEncoder 🎉

Y no es que hayamos empatado por poco: las dos están resolviendo exactamente el mismo problema con el mismo método. La diferencia de dos diezmilésimas es que scikit-learn usa un optimizador más listo que el nuestro y llega a un sitio ligerísimamente distinto.

Ese es el punto de haber escrito todo esto a mano. Ya no hay ninguna parte de una regresión logística que no sepas hacer tú, y lo que viene en el capítulo 6 es la misma idea con capas 💪

Ejercicios

w = np.zeros(2); b = 0.0
for i in range(201):
    p = sigmoide(escalado @ w + b)
    if i % 40 == 0:
        print(f'vuelta {i:3d}  pérdida {perdida(p, objetivo):.4f}')
    e = p - objetivo
    w -= 0.1 * (escalado.T @ e) / len(objetivo)
    b -= 0.1 * e.mean()

vuelta   0  pérdida 0.6931
vuelta  40  pérdida 0.5858
vuelta  80  pérdida 0.5639
vuelta 120  pérdida 0.5574
vuelta 160  pérdida 0.5552
vuelta 200  pérdida 0.5544

print('pesos encontrados:', np.round(w, 4))
print('sesgo            :', round(float(b), 4))
print()
print('yo generé con 0.002 por sol de monto y 0.8 por punto de satisfacción')
print('en unidades escaladas eso es:',
      np.round([0.002 * crudo[:, 0].std(), 0.8 * crudo[:, 1].std()], 4))

pesos encontrados: [0.6997 1.0115]
sesgo            : 0.171

yo generé con 0.002 por sol de monto y 0.8 por punto de satisfacción
en unidades escaladas eso es: [0.5728 0.8234]

for vueltas in [200, 1000, 5000, 20000]:
    ww = np.zeros(2); bb = 0.0
    for _ in range(vueltas):
        e = sigmoide(escalado @ ww + bb) - objetivo
        ww -= 0.1 * (escalado.T @ e) / len(objetivo)
        bb -= 0.1 * e.mean()
    print(f'{vueltas:6d} vueltas: pesos {np.round(ww, 4)}  '
          f'pérdida {perdida(sigmoide(escalado @ ww + bb), objetivo):.4f}')

   200 vueltas: pesos [0.699  1.0106]  pérdida 0.5544
  1000 vueltas: pesos [0.763  1.0924]  pérdida 0.5538
  5000 vueltas: pesos [0.763  1.0925]  pérdida 0.5538
 20000 vueltas: pesos [0.763  1.0925]  pérdida 0.5538

for filas in [500, 5000, 50000]:
    g = np.random.default_rng(0)
    m = g.normal(800, 300, filas)
    s = g.normal(3.5, 1.0, filas)
    obj = (g.random(filas) < sigmoide(0.002 * (m - 800) + 0.8 * (s - 3.5))).astype(float)
    c = np.column_stack([m, s])
    esc = (c - c.mean(axis=0)) / c.std(axis=0)

    ww = np.zeros(2); bb = 0.0
    for _ in range(20000):
        e = sigmoide(esc @ ww + bb) - obj
        ww -= 0.1 * (esc.T @ e) / filas
        bb -= 0.1 * e.mean()

    verdad = np.array([0.002 * c[:, 0].std(), 0.8 * c[:, 1].std()])
    print(f'{filas:6d} filas: encontrado {np.round(ww, 4)}  '
          f'verdadero {np.round(verdad, 4)}')

   500 filas: encontrado [0.5796 0.6495]  verdadero [0.6082 0.7504]
  5000 filas: encontrado [0.635  0.8078]  verdadero [0.5972 0.8005]
 50000 filas: encontrado [0.594  0.8179]  verdadero [0.6008 0.7991]

for vueltas in [100, 500, 1500, 3000]:
    ww = np.zeros(T_tr.shape[1]); bb = 0.0
    for _ in range(vueltas):
        e = sigmoide(T_tr @ ww + bb) - objetivo_v
        ww -= 0.5 * (T_tr.T @ e) / len(objetivo_v)
        bb -= 0.5 * e.mean()
    print(f'{vueltas:5d} vueltas: AUC {roc_auc_score(yv_te, sigmoide(T_te @ ww + bb)):.4f}')

  100 vueltas: AUC 0.7189
  500 vueltas: AUC 0.7213
 1500 vueltas: AUC 0.7216
 3000 vueltas: AUC 0.7216

perfecto = np.array([10., 0., 0.])
malo = np.array([-10., 0., 0.])
y_facil = (X[:, 0] > 0).astype(float)

for nombre, pesos in [('acierta', perfecto), ('se equivoca', malo)]:
    e = sigmoide(X @ pesos) - y_facil
    g = X.T @ e / len(y_facil)
    print(f'{nombre:12} pérdida {perdida(sigmoide(X @ pesos), y_facil):.4f}  '
          f'gradiente {np.round(g, 4)}')

acierta      pérdida 0.0628  gradiente [-0.0067 -0.0028  0.0041]
se equivoca  pérdida 8.0274  gradiente [-0.7963  0.0155  0.0704]

lote_rng = np.random.default_rng(7)
w_lote = np.zeros(2); b_lote = 0.0
for i in range(201):
    idx = lote_rng.choice(len(objetivo), 32, replace=False)
    e = sigmoide(escalado[idx] @ w_lote + b_lote) - objetivo[idx]
    w_lote -= 0.1 * (escalado[idx].T @ e) / len(idx)
    b_lote -= 0.1 * e.mean()

print('con lotes de 32 :', np.round(w_lote, 4),
      f'pérdida {perdida(sigmoide(escalado @ w_lote + b_lote), objetivo):.4f}')
print('con las 500     :', np.round(w, 4),
      f'pérdida {perdida(sigmoide(escalado @ w + b), objetivo):.4f}')

con lotes de 32 : [0.7054 1.0285] pérdida 0.5547
con las 500     : [0.6997 1.0115] pérdida 0.5543

lote_rng.choice(len(objetivo), 600, replace=False)

ValueError: Cannot take a larger sample than population when replace is False

formula = float((sigmoide(X @ w[:3] if len(w) == 3 else X @ np.array([0.2, -0.1, 0.4]) + b) - y).mean())
numerico_b = (perdida_de(np.array([0.2, -0.1, 0.4]), b + h)
              - perdida_de(np.array([0.2, -0.1, 0.4]), b - h)) / (2 * h)
print('por fórmula:', round(formula, 8))
print('moviendo   :', round(numerico_b, 8))
print('diferencia :', round(abs(formula - numerico_b), 12))

por fórmula: 0.03743041
moviendo   : 0.03743041
diferencia : 0.0

w_mal = np.zeros(2); b_mal = 0.0
for i in range(8):
    p = sigmoide(escalado @ w_mal + b_mal)
    print(f'vuelta {i}: pérdida {perdida(p, objetivo):10.4f}  pesos {np.round(w_mal, 3)}')
    e = p - objetivo
    w_mal -= 60.0 * (escalado.T @ e) / len(objetivo)
    b_mal -= 60.0 * e.mean()

vuelta 0: pérdida     0.6931  pesos [0. 0.]
vuelta 1: pérdida     2.1362  pesos [ 7.897 11.991]
vuelta 2: pérdida     1.0199  pesos [4.257 4.276]
vuelta 3: pérdida     1.6548  pesos [-2.816  0.363]
vuelta 4: pérdida     4.7848  pesos [24.973  8.635]
vuelta 5: pérdida     2.6199  pesos [10.694 14.776]
vuelta 6: pérdida     1.5089  pesos [6.047 7.431]
vuelta 7: pérdida     0.6832  pesos [0.707 1.659]

def perdida_sin_clip(p, y):
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))

print(perdida_sin_clip(np.array([1.0, 0.0]), np.array([1.0, 0.0])))

nan

Lo que te llevas

• 📏 La pérdida de entropía cruzada castiga estar seguro y equivocado (4,6052) mucho más que dudar (0,6931).
• 🧮 El gradiente de una neurona cabe en una línea: X.T @ (predicho - real) / n.
• 🔬 Se comprueba moviendo cada peso a mano, y coincide con 6,7e-12 de diferencia.
• 🔽 El gradiente apunta hacia donde el error crece, así que se resta.
• ⚖️ Sin escalar no hay ningún paso que sirva: desde 0,0001 explota, y con 0,00001 la pérdida se mueve una diezmilésima en 200 vueltas.
• ⏱️ El 99% del trabajo se hace al principio: entre 5.000 y 20.000 vueltas la pérdida solo baja 0,0005.
• 🚀 Con lotes de 32 filas se llega casi al mismo sitio mirando una quinceava parte de los datos.
• 🪃 Pesos que saltan de signo son la firma de un paso demasiado grande.

En el capítulo 6 llevamos esto a una red con capas, que es donde el cálculo se complica y aparece la retropropagación.

Que tengas lindo día! 🌸