Resolución de sistemas lineales

Un sistema de ecuaciones lineales se puede escribir en forma matricial como Ax = b, donde A es la matriz de coeficientes, x es el vector de incógnitas y b es el vector del segundo miembro. Por ejemplo, dos ecuaciones con dos incógnitas se pueden expresar como una ecuación matricial 2×2. Resolver el sistema significa encontrar x tal que A multiplicada por x sea igual a b. np.linalg.solve(A, b) de NumPy lo gestiona de forma eficiente y precisa.

Matemáticamente, la solución de Ax = b es x = A⁻¹b. Sin embargo, calcular np.linalg.inv(A) @ b es más lento y menos preciso que llamar directamente a np.linalg.solve(A, b). solve() utiliza internamente la descomposición LU, que requiere menos operaciones y acumula menos errores de coma flotante. Regla práctica: si está resolviendo x en Ax=b, use siempre solve, nunca inv.

np.linalg.solve(A, b) requiere que A sea cuadrada y no singular. Si A es singular (su determinante es cero), la función genera np.linalg.LinAlgError: Singular matrix. A es singular cuando sus filas son linealmente dependientes; por ejemplo, cuando una ecuación es simplemente un múltiplo de otra. Compruebe siempre que el sistema esté bien planteado antes de llamar a solve. Si A no es cuadrada, use np.linalg.lstsq().

np.linalg.solve aplica internamente la descomposición LU: factoriza A en una matriz triangular inferior L y una matriz triangular superior U, de modo que A = LU. Resolver el sistema se convierte entonces en resolver dos sistemas triangulares más sencillos (mediante sustitución hacia delante y hacia atrás), lo que es mucho más eficiente que repetir la eliminación gaussiana desde cero. Este es el mismo algoritmo que utilizan MATLAB, Fortran LAPACK y prácticamente todas las bibliotecas de cálculo científico.

np.linalg.solve(A, B) puede gestionar una matriz B con varias columnas y resolver cada columna simultáneamente. Si B tiene dimensión (n, k), el resultado X tiene dimensión (n, k), donde cada columna de X es la solución correspondiente a la columna de B. Esto resulta útil en simulaciones físicas y aprendizaje automático, cuando necesita resolver la misma matriz de coeficientes con muchos segundos miembros diferentes.

Cuando el sistema está sobredeterminado (hay más ecuaciones que incógnitas), no existe una solución exacta. np.linalg.lstsq(A, b) encuentra el valor de x que minimiza la suma de los cuadrados de los residuos ||Ax - b||²; esto es exactamente lo que calcula la regresión lineal. Devuelve la solución x, los residuos, el rango y los valores singulares. El argumento rcond=None suprime una advertencia de obsolescencia y establece el umbral para tratar los valores singulares como cero.

Las ecuaciones normales de los mínimos cuadrados ordinarios son A.T @ A @ x = A.T @ b, que constituyen un sistema cuadrado resoluble con np.linalg.solve. Esto produce el mismo resultado que lstsq cuando A.T @ A está bien condicionada. En la práctica, se prefiere lstsq porque utiliza SVD y es estable incluso cuando las columnas son casi colineales, pero las ecuaciones normales ayudan a comprender las matemáticas subyacentes de la regresión.

Antes de resolver el sistema, es recomendable comprobar el número de condición con np.linalg.cond(A). Un número de condición cercano a 1 significa que el sistema está bien condicionado y que la solución es fiable. Un número de condición muy grande (por ejemplo, 1e12) significa que la matriz es casi singular: pequeños cambios en b producen grandes cambios en x y la solución puede no ser fiable numéricamente. En ese caso, considere la regularización (regresión ridge) o el uso de un solucionador más robusto.

En sistemas muy grandes donde A tiene principalmente entradas cero (una matriz dispersa), almacenar la matriz completa desperdicia memoria y usar np.linalg.solve denso es lento. SciPy proporciona scipy.sparse y scipy.sparse.linalg.spsolve, que aprovechan la dispersión para trabajar con sistemas enormes, como los de modelos de elementos finitos o problemas de redes. Comprender cuándo su matriz es dispersa es una consideración importante de rendimiento en el cálculo científico.

Una vez que tenga el vector solución x de np.linalg.solve(A, b), verifíquelo siempre calculando A @ x y comparándolo con b mediante np.allclose(A @ x, b). En la aritmética de punto flotante, rara vez se obtiene una igualdad exacta; allclose comprueba que la diferencia esté dentro de una pequeña tolerancia. El residuo ||Ax - b|| debería estar cerca del épsilon de máquina multiplicado por la escala del problema en un sistema bien condicionado.

Los solucionadores lineales aparecen en todo el aprendizaje automático: la regresión lineal resuelve las ecuaciones normales, los procesos gaussianos resuelven sistemas de matrices de kernel y los filtros de Kalman resuelven ecuaciones de actualización de estados. Al comprender np.linalg.solve y lstsq, podrá implementar estos algoritmos desde cero y depurarlos cuando los resultados no sean los esperados. El módulo linalg de NumPy le proporciona las mismas herramientas numéricas que se utilizan en bibliotecas científicas de producción.

En esta lección ha aprendido que np.linalg.solve(A, b) resuelve eficazmente sistemas lineales cuadrados mediante la descomposición LU, np.linalg.lstsq() encuentra la solución de mínimos cuadrados para sistemas sobredeterminados y el número de condición mide la estabilidad numérica; los valores grandes advierten de matrices casi singulares. A continuación, exploraremos los valores propios y la SVD, fundamentales para PCA y muchos algoritmos de aprendizaje automático.