Transpuesta e inversa, y por qué casi nunca se invierte

Dos operaciones que van juntas por un motivo 🌸

La transpuesta ya la viste de pasada en el capítulo 1. La inversa es nueva. Van en el mismo capítulo porque la fórmula más famosa de todo machine learning las usa a las dos seguidas, y al final del capítulo la vas a poder leer entera.

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 transpuesta gira la tabla

print('X.shape  :', X.shape)
print('X.T.shape:', X.T.shape)
print('(X.T).T es X:', np.array_equal(X.T.T, X))

X.shape  : (2238, 3)
X.T.shape: (3, 2238)
(X.T).T es X: True

Girarla dos veces te devuelve la de partida, que es lo que uno esperaría de algo que se llama "girar".

Y para qué sirve: sirve para que las formas encajen. En el capítulo 5 viste que W @ X reventaba. La transpuesta es lo que arregla eso, y el 90% de las veces que la vas a usar es por eso mismo.

X traspuesta por X, la matriz que aparece en todas partes

Vamos a trabajar con dos columnas para que los números se lean:

U = X[:, [0, 2]]
G = U.T @ U
print('U.shape:', U.shape)
print('G = U.T @ U')
print(G)

U.shape: (2238, 2)
G = U.T @ U
[[380021.  78592.]
 [ 78592.  24915.]]

Fíjate en la forma: entra (2, 2238) por (2238, 2) y sale (2, 2). Los 2.238 se gastaron.

Eso es lo que hace X^⊤X y por eso está en todas partes: convierte una tabla enorme en una matriz chiquita que solo tiene el tamaño de tus columnas. Da igual que tengas mil filas o un millón, esa matriz sigue siendo de 2 por 2.

Y fíjate en otra cosa: los dos 78592 de las esquinas son el mismo número. Sale simétrica siempre, y no es casualidad. Es el producto punto de la columna 0 con la columna 1, que da lo mismo en los dos sentidos.

La inversa deshace

Gi = np.linalg.inv(G)
print('inversa de G:')
print(np.round(Gi, 8))
print('G @ Gi redondeado a 10:')
print(np.round(G @ Gi, 10))

inversa de G:
[[ 7.5700e-06 -2.3880e-05]
 [-2.3880e-05  1.1545e-04]]
G @ Gi redondeado a 10:
[[ 1. -0.]
 [-0.  1.]]

La inversa es el dividir de las matrices. Con números, 5 por un quinto da 1. Con matrices, A por su inversa da la identidad, que es el 1 de este mundo.

Y ojo con ese -0. de la salida. No es un error de NumPy: es coma flotante. El resultado no es exactamente cero, es algo como diez a la menos diecisiete, y al redondear sale con signo. Por eso nunca compares matrices con ==:

print(np.allclose(G @ Gi, np.eye(2)))

True

np.allclose pregunta "¿son iguales dentro de lo razonable?" y esa es la pregunta correcta cuando hay decimales de por medio.

Cuándo NO hay inversa, y por qué da más miedo de lo que parece

Le pego a X una columna que son las unidades multiplicadas por dos, o sea información que ya estaba:

D = np.column_stack([X[:, 0], X[:, 0] * 2])
GD = D.T @ D
print('GD ='); print(GD)
print('rango:', np.linalg.matrix_rank(GD), 'de 2')

GD =
[[ 380021.  760042.]
 [ 760042. 1520084.]]
rango: 1 de 2

Rango 1 de 2 quiere decir que ahí solo hay una columna de información. Esta matriz no tiene inversa: se dice que es singular.

Y ahora la parte incómoda:

inv = np.linalg.inv(GD)
print('inv ='); print(inv)
print('GD @ inv, que tendria que ser la identidad:')
print(np.round(GD @ inv, 4))

inv =
[[-1.71798692e+10  8.58993459e+09]
 [ 8.58993459e+09 -4.29496730e+09]]
GD @ inv, que tendria que ser la identidad:
[[ 1.     -0.1376]
 [ 2.     -0.2752]]

No dio error. Te devolvió una matriz con toda la tranquilidad del mundo.

Pero mírala: números en el orden de diez mil millones, cuando la matriz de partida tenía cientos de miles. Y la comprobación no miente: eso no es la identidad ni de lejos. Debajo tendría que salir uno, cero, cero, uno, y sale un dos y un menos cero coma veintisiete.

Es basura. Basura que se puede imprimir, meter en un informe y presentar en una reunión 🙃

A veces sí avisa. Si la dependencia sale tan limpia que el determinante da cero exacto, NumPy corta:

Xd = np.column_stack([X, X[:, 0] * 2])
np.linalg.inv(Xd.T @ Xd)

LinAlgError: Singular matrix

Pero no puedes contar con eso. Depende de los números que te tocaron, y en datos de verdad el aviso casi nunca llega.

Si en tu tabla metes "precio en soles" y "precio en dólares" con un tipo de cambio fijo, acabas de hacer exactamente esto. La forma de detectarlo antes, sin depender de que algo reviente, es mirar el rango, que es el capítulo 8.

Y ahora la fórmula famosa

Léela de derecha a izquierda: agarra tus datos y tu objetivo, gíralos para que encajen, y deshaz. Lo que sale son los pesos.

Eso es una regresión lineal resuelta de una vez, sin entrenar ni iterar nada. Se llaman las ecuaciones normales.

Pero en código no se invierte

y = X[:, 1]
U = X[:, [0, 2]]
G = U.T @ U

print('con solve:', np.round(np.linalg.solve(G, U.T @ y), 4))
print('con inv  :', np.round(np.linalg.inv(G) @ (U.T @ y), 4))

con solve: [ 30.937  118.0057]
con inv  : [ 30.937  118.0057]

Dan lo mismo, y aun así usa siempre solve. Dos motivos: hace menos cuentas, y arrastra menos error de redondeo. Cuando la matriz está casi en el borde de ser singular, inv te devuelve números que parecen razonables y no lo son.

La regla práctica: si escribiste inv en código que va a producción, casi seguro querías solve 🙂

Lo que te llevas

• La transpuesta cambia filas por columnas y se usa sobre todo para que las formas encajen.
• X^⊤X convierte una tabla enorme en una matriz del tamaño de tus columnas, y sale simétrica.
• La inversa deshace, y la identidad es el 1 de las matrices.
• Compara con np.allclose, nunca con ==.
• Singular quiere decir que sobra información.
• En código: solve, no inv.

Comprueba que se entendió

Ejercicios

G = X.T @ X
print('es simetrica:', np.allclose(G, G.T))
print('forma       :', G.shape)

es simetrica: True
forma       : (3, 3)

I = np.eye(3)
print(np.array_equal(X @ I, X))
print(np.array_equal(I @ X.T, X.T))

True
True

D = np.column_stack([X[:, 0], X[:, 0] * 2])
print('las dos primeras filas:')
print(D[:2])
print('determinante:', float(np.linalg.det(D.T @ D)))

las dos primeras filas:
[[10. 20.]
 [10. 20.]]
determinante: -8.848053403198714e-05

Preguntas frecuentes