Resolución numérica para IA: usa solve sin inversas explícitas

Un paso de calibración de IA puede necesitar recuperar dos parámetros desconocidos a partir de dos respuestas medidas. Las ecuaciones son lo bastante pequeñas para escribirlas en una pizarra, pero el código de producción aún tiene que decidir cómo resolverlas y cómo comprobar la respuesta. Un vector que parece plausible no demuestra que los datos lo hayan determinado de forma fiable.

Esta lección desarrolla esa decisión con NumPy. Resolverás un sistema lineal cuadrado, compararás numpy.linalg.solve con la ruta de la inversa explícita, calcularás un residuo, medirás el condicionamiento y reconocerás el límite singular en el que no existe una estimación única de los parámetros.

Resultado de aprendizaje: Usa numpy.linalg.solve y explica por qué normalmente se evita la inversión explícita.

La lección 22: Identidad y transformaciones inversas mostró qué significa una inversa desde el punto de vista matemático. La lección 24: La forma escalonada por filas revela la estructura mostró cómo los pivotes y las variables libres exponen lo que un sistema puede determinar. Aquí conectamos ese álgebra exacta con la computación de precisión finita. El sitio web oficial de apoyo de Mathematics for Machine Learning sitúa este trabajo en la sección 2.3 del PDF del libro; el escenario de calibración, los números, el código, la visualización y las evaluaciones que siguen son originales.

Un sistema lineal es un contrato de recuperación de parámetros

Escribe el vector de parámetros desconocidos como

x = \begin{bmatrix}x_1\\x_2\end{bmatrix}\in\mathbb{R}^2.

El vector tiene dos coordenadas, por lo que cada fila de una matriz de coeficientes $A$ debe contener dos coeficientes. Si dos observaciones producen el vector de términos independientes $b$ , un sistema de recuperación de parámetros tiene la forma

A x=b, \qquad A\in\mathbb{R}^{2\times 2}, \quad b\in\mathbb{R}^{2}, \quad x\in\mathbb{R}^{2}.

En una pequeña cabeza de calibración, $x_1$ podría ser una escala aplicada a una señal y $x_2$ una corrección parecida a un desplazamiento aplicada a otra. De forma más general, la misma notación describe una capa lineal pequeña, una corrección de sensor o un subproblema de recuperación de parámetros dentro de una canalización de IA más grande. Las filas de $A$ son los patrones de coeficientes observados; $b$ contiene las respuestas que esos patrones deberían producir.

Si $A$ es cuadrada e invertible, el álgebra de la lección anterior da

x=A^{-1}b.

Esa identidad es útil para razonar. No significa que un programa deba construir primero $A^{-1}$ cada vez que necesita $x$ .

Ejemplo trabajado: recuperar dos parámetros de calibración

Considera el sistema deliberadamente pequeño

A=\begin{bmatrix}3&1\\1&2\end{bmatrix}, \qquad b=\begin{bmatrix}8\\1\end{bmatrix}.

Las dos ecuaciones son

3x_1+x_2=8, \qquad x_1+2x_2=1.

Resuélvelo por sustitución antes de pedirle a NumPy que haga el mismo trabajo. La primera ecuación da $x_2=8-3x_1$ . Sustituye esa expresión en la segunda:

x_1+2(8-3x_1)=1 \Longrightarrow -5x_1=-15 \Longrightarrow x_1=3.

Después $x_2=8-3(3)=-1$ , así que

x=\begin{bmatrix}3\\-1\end{bmatrix}.

Comprueba las dos ecuaciones originales, no solo la que has reordenado:

A x = \begin{bmatrix}3&1\\1&2\end{bmatrix} \begin{bmatrix}3\\-1\end{bmatrix} = \begin{bmatrix}8\\1\end{bmatrix}=b.

El mismo cálculo en NumPy

import numpy as np

A = np.array([[3.0, 1.0],
              [1.0, 2.0]])
b = np.array([8.0, 1.0])

x_solve = np.linalg.solve(A, b)
x_inverse = np.linalg.inv(A) @ b

print(x_solve)
print(x_inverse)
print(np.linalg.norm(A @ x_solve - b))
print(np.linalg.norm(A @ x_inverse - b))

En el entorno de ejecución incluido, con NumPy 2.3.5, la salida verificada es:

[ 3. -1.]
[ 3. -1.]
2.220446049250313e-16
9.155133597044475e-16

Los valores diminutos distintos de cero son normas de residuos de punto flotante. En aritmética exacta ambos residuos son cero y ambos métodos representan la misma respuesta matemática. La ligera diferencia en los últimos dígitos recuerda que los resultados de punto flotante deben compararse con una tolerancia, no mediante igualdad exacta.

Por qué `solve` suele ser el paso intermedio adecuado

La expresión matemática directa $A^{-1}b$ oculta dos tareas computacionales diferentes:

Ruta	Qué calcula	Por qué importa para un término independiente
`np.linalg.solve(A, b)`	Factoriza $A$ y aplica los factores a $b$	Calcula solo el vector solicitado.
`np.linalg.inv(A) @ b`	Construye cada columna de $A^{-1}$ y después la multiplica por $b$	Calcula mucha más información de la que necesita esta única resolución.

Para un solo término independiente, solve evita materializar una inversa completa. Para varios términos independientes almacenados como columnas de $B$ , usa np.linalg.solve(A, B) para que la misma matriz de coeficientes se factorice una vez mientras se calculan todas las soluciones solicitadas. Conserva una inversa solo cuando la aplicación realmente necesita la inversa; no es un sustituto predeterminado de una resolución.

La documentación actual de NumPy para solve dice que la matriz de coeficientes debe ser cuadrada y de rango completo, y documenta que la rutina usa una implementación LAPACK _gesv. La documentación de inv también advierte que una matriz mal condicionada puede producir resultados inexactos aunque la inversión no genere un error. Por eso la elección de implementación y los diagnósticos deben formar parte del mismo modelo mental.

La diferencia de tiempo de ejecución depende del hardware y de las bibliotecas, así que no es un benchmark portable. Como ilustración, la siguiente ejecución en el host de autoría usó una matriz aleatoria fija de $64\times64$ , cinco repeticiones de medición y veinte llamadas por repetición:

from timeit import repeat

rng = np.random.default_rng(7)
A_timing = rng.normal(size=(64, 64)) + 64 * np.eye(64)
b_timing = rng.normal(size=64)

solve_time = min(
    repeat(lambda: np.linalg.solve(A_timing, b_timing), repeat=5, number=20)
) / 20
inverse_time = min(
    repeat(lambda: np.linalg.inv(A_timing) @ b_timing, repeat=5, number=20)
) / 20

print(solve_time * 1e6, inverse_time * 1e6, inverse_time / solve_time)

La ejecución verificada imprimió aproximadamente 15.2 µs, 52.2 µs y 3.43× para el tiempo de resolución, el tiempo de la ruta de la inversa y la proporción. Considera esas cifras como una forma reproducible de la comparación, no como una promesa de rendimiento: la razón algorítmica para preferir solve es que evita calcular una matriz que no se usa, y una biblioteca BLAS o un tamaño de matriz diferente pueden cambiar la proporción.

Un residuo es una comprobación de ajuste, no un certificado de confianza

Para el ejemplo trabajado, calcula el residuo de la ruta directa:

r=b-A\hat{x} = \begin{bmatrix}8\\1\end{bmatrix} - \begin{bmatrix}8\\1\end{bmatrix} = \begin{bmatrix}0\\0\end{bmatrix}.

El programa de punto flotante informa de una norma cercana a $10^{-15}$ en vez de exactamente cero. Es una buena evidencia de que el vector devuelto ajusta este $b$ . Sin embargo, imagina que los coeficientes o las respuestas proceden de mediciones de sensores redondeadas. Un solucionador puede encontrar un vector que ajuste muy bien esos valores redondeados y, aun así, cambiar sustancialmente cuando las mediciones varían dentro de su incertidumbre.

El condicionamiento mide la sensibilidad

El número de condición compara cuánto puede amplificar una matriz los cambios relativos de entrada. Para una matriz invertible y la norma 2,

\kappa_2(A)=\lVert A\rVert_2,\lVert A^{-1}\rVert_2 =\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)},

donde $\sigma_{\max}$ y $\sigma_{\min}$ son los valores singulares mayor y menor. Un número de condición cercano a $1$ indica direcciones equilibradas. Un valor grande significa que algún cambio pequeño en las entradas puede provocar un cambio relativo mucho mayor en la solución. Si el valor singular menor es cero, la matriz es singular y el número de condición es infinito.

Para un término independiente perturbado, una cota útil tiene la forma

\frac{\lVert\delta x\rVert}{\lVert x\rVert} \lesssim \kappa(A)\frac{\lVert\delta b\rVert}{\lVert b\rVert},

bajo las hipótesis habituales para un sistema no singular. La desigualdad advierte de una posible amplificación; no promete que toda perturbación alcance el peor caso.

Ahora usa filas casi dependientes:

B=\begin{bmatrix}1&1\\1&1.0001\end{bmatrix}, \qquad x_{\text{true}}=\begin{bmatrix}2\\-1\end{bmatrix}, \qquad b=B x_{\text{true}}=\begin{bmatrix}1\\0.9999\end{bmatrix}.

La segunda fila es casi una copia de la primera. Con el $b$ exacto que se muestra, la resolución recupera $[2,-1]$ . Cambia solo la segunda respuesta en $10^{-6}$ y vuelve a resolver. NumPy informa de lo siguiente:

condition number: 40002.000074915224
original solution: [ 2. -1.]
perturbed solution: [ 1.99 -0.99]
relative input change: 7.071421374297804e-07
relative solution change: 0.006324555320519327
perturbed residual norm: 0.0

La entrada cambió aproximadamente $7.1\times10^{-7}$ de forma relativa, mientras que el vector de parámetros cambió aproximadamente $6.3\times10^{-3}$ de forma relativa. Es aproximadamente una amplificación de cuatro órdenes de magnitud expresada por el número de condición. La solución perturbada sigue teniendo residuo cero para los datos perturbados: resolvió correctamente las nuevas ecuaciones, pero los datos no identificaron un vector de parámetros estable.

Mapa estático de estabilidad

Un residuo pequeño puede ocultar una solución sensible

El residuo comprueba la ecuación resuelta. El número de condición avisa cuánto puede moverse la solución cuando cambian los datos de entrada.

Los tres casos pueden producir un residuo pequeño para su segundo miembro. El condicionamiento determina si ese ajuste también es una estimación de parámetros fiable.

Cantidades exactas representadas en el mapa de estabilidad
Caso	Número de condición	Cambio relativo de entrada → solución	Residuo e interpretación
Bien condicionado	2.62	1.2 × 10⁻⁷ → 2.0 × 10⁻⁷	Residuos de alrededor de 10⁻¹⁵; ambas rutas coinciden en este ejemplo bien condicionado.
Mal condicionado	4.0002 × 10⁴	7.1 × 10⁻⁷ → 6.3 × 10⁻³	Residuo 0 para los datos perturbados, pero un cambio de entrada de 7,1 × 10⁻⁷ mueve x en 6,3 × 10⁻³.
Singular	∞	No aplicable	κ₂(A) es infinito; solve e inverse generan un error de matriz singular en lugar de devolver una respuesta única.

Lee ambos diagnósticos juntos: el residuo pregunta “¿este x ajusta b?”; el condicionamiento pregunta “¿cuánto cambiaría x si b o A fueran ligeramente distintos?”

La visualización separa las dos preguntas. El residuo pregunta si el candidato ajusta las ecuaciones que están actualmente en memoria. El número de condición pregunta si un pequeño cambio en las ecuaciones podría mover ese candidato. Un sistema de calibración de IA necesita ambas respuestas.

Caso límite: singular significa “no hay una solución única”

Considera

C=\begin{bmatrix}1&2\\2&4\end{bmatrix}.

La segunda fila es el doble de la primera, así que $C$ tiene rango $1$ , no $2$ . Con $b=[3,6]^\mathsf{T}$ , las ecuaciones son compatibles pero tienen infinitas soluciones: sirve cualquier $(x_1,x_2)$ que satisfaga $x_1+2x_2=3$ . Con $b=[3,7]^\mathsf{T}$ , son incompatibles porque la segunda ecuación exige el doble del primer miembro izquierdo, pero proporciona un término independiente diferente.

Ninguna de las dos situaciones tiene un vector único que np.linalg.solve pueda devolver. En el entorno incluido, tanto np.linalg.solve(C, [3, 6]) como np.linalg.inv(C) generan numpy.linalg.LinAlgError: Singular matrix. La documentación de NumPy para solve especifica el requisito de rango completo y remite a numpy.linalg.lstsq cuando una solución de mínimos cuadrados es adecuada para un problema no cuadrado o deficiente en rango. Esa es una cuestión distinta de resolver un sistema cuadrado invertible, así que no sustituyas silenciosamente un método por el otro.

Un flujo de trabajo disciplinado de NumPy para sistemas de IA

Para un pequeño subproblema de recuperación de parámetros, un flujo de trabajo defendible es:

Comprueba el contrato. Confirma que coinciden las dimensiones de $A$ , $b$ y el vector de parámetros previsto.
Resuelve directamente. Usa np.linalg.solve(A, b) para un sistema cuadrado de rango completo; usa un método diseñado para mínimos cuadrados cuando el problema no esté bien determinado.
Mide el residuo. Calcula np.linalg.norm(A @ x - b) y compáralo con la escala de $b$ , no solo con un umbral absoluto arbitrario.
Inspecciona el condicionamiento. Calcula np.linalg.cond(A) y compara la posible sensibilidad con la precisión de las mediciones y con el número de dígitos que piensas usar después.
Trata el fallo como información. Un error de singularidad, un número de condición enorme o un residuo grande con respecto a $b$ debería activar una investigación de los datos, las características, el escalado o el diseño del modelo.

En una tarea de calibración o recuperación de parámetros, este flujo convierte un vector de pesos en un resultado auditable. Un residuo pequeño puede respaldar la afirmación de que la implementación resolvió las ecuaciones suministradas. Un número de condición moderado y una comprobación de perturbaciones aportan evidencia de que esas ecuaciones restringen los parámetros con suficiente fuerza para el uso previsto. Ningún diagnóstico demuestra que el sistema de IA más amplio sea exacto; responden a preguntas de implementación más concretas que se pueden comprobar.

La próxima lección se aparta de los algoritmos numéricos y pregunta qué operaciones algebraicas son válidas sobre el sistema numérico que utilizan los tensores de IA habituales. Por ahora, conserva esta regla práctica: resuelve el sistema que tienes, comprueba el residuo que obtienes y pregunta al condicionamiento si un sistema cercano contaría una historia de parámetros diferente.