A <- matrix(c(4, 1, -1, 1, 6, 2, -1, 2, 5), nrow = 3, byrow = TRUE)
b <- c(6, 15, 12)
# Resolver Ax = b
x <- solve(A, b)
print(x)[1] 1.651685 1.516854 2.123596# Verificación (debería dar b)
A %*% x [,1]
[1,] 6
[2,] 15
[3,] 122 Resolución de grandes sistemas de ecuaciones
2.1 El problema de la escala
En el análisis de datos moderno, la modelización predictiva y la ingeniería, raramente nos enfrentamos a sistemas de 3x3 ecuaciones. Lo habitual es tratar con sistemas donde el número de variables, \(n\), puede ser de miles, millones o incluso miles de millones.
Resolver un sistema \(A\symbfit{x} = \symbfit{b}\) con \(n=10.000\) es trivial para un ordenador moderno si la matriz \(A\) es densa (casi todas sus entradas son no nulas). Pero, ¿qué ocurre si \(n=10.000.000\)?
Si \(A\) es una matriz densa de \(10^7 \times 10^7\), necesitaríamos almacenar \((10^7)^2 = 10^{14}\) valores. Si cada valor ocupa 8 bytes (un numeric en R), el espacio de memoria requerido sería:
\(10^{14} \text{ valores} \times 8 \text{ bytes/valor} \approx 8 \times 10^{14} \text{ bytes} \approx 800 \text{ Terabytes}\)
Ningún ordenador personal puede almacenar esa matriz en RAM. Además, los métodos directos clásicos, como la eliminación de Gauss, tienen un coste computacional de \(O(n^3)\) operaciones. Para \(n=10^7\), esto es \(10^{21}\) operaciones, una cantidad intratable.
La solución a este problema proviene de dos frentes complementarios. Por un lado, aprovechamos que la mayoría de los sistemas del mundo real son matrices dispersas, lo que nos permite usar estructuras de datos que almacenan solo los valores no nulos. Por otro lado, recurrimos a métodos numéricos avanzados, como los iterativos (Jacobi, Gauss-Seidel), que refinan una solución aproximada en lugar de intentar resolver el sistema de forma directa y costosa.
En este capítulo, exploraremos las herramientas teóricas y computacionales en R para abordar ambos frentes.
2.2 Métodos directos y factorización
Los métodos directos resuelven el sistema en un número finito y predecible de pasos. El más conocido es la Eliminación de Gauss, que transforma el sistema en uno triangular superior equivalente, fácil de resolver mediante sustitución regresiva.
NotaCarl Friedrich Gauss
Carl Friedrich Gauss (1777-1855) no inventó el método (ya se conocía en China en el siglo I a.C.), pero lo popularizó y formalizó en Europa. Su contribución clave fue el análisis sistemático de los errores de medición, lo que le llevó al método de mínimos cuadrados, intrínsecamente ligado a la resolución de sistemas.
En R, la función solve() es una “caja negra” altamente optimizada que implementa un método directo.
2.2.1 El coste de solve(): factorización LU
La función solve(A, b) no calcula la inversa \(A^{-1}\) y luego multiplica \(A^{-1}\symbfit{b}\), ya que calcular la inversa es computacionalmente más caro e inestable.
En su lugar, la mayoría de los solvers modernos (incluido el de R) se basan en una factorización de la matriz. La más común es la factorización LU.
Definición 2.1 (Factorización LU) Para una matriz \(A\) (bajo ciertas condiciones), la factorización LU encuentra una matriz triangular inferior \(L\) (Lower), una matriz triangular superior \(U\) (Upper) y (a menudo) una matriz de permutación \(P\) (que reordena las filas para estabilidad numérica) tales que:
\[PA = LU\]
Resolver \(A\symbfit{x} = \symbfit{b}\) se convierte en \(PA\symbfit{x} = P\symbfit{b}\), es decir, \(LU\symbfit{x} = P\symbfit{b}\). Este problema se resuelve en dos pasos (mucho más rápidos):
Resolver el sistema \(LU\symbfit{x} = P\symbfit{b}\) se convierte en un proceso de dos pasos muy eficiente. Primero, realizamos una sustitución progresiva para resolver \(L\symbfit{y} = P\symbfit{b}\) y encontrar \(\symbfit{y}\). Acto seguido, usamos ese resultado en una sustitución regresiva para resolver \(U\symbfit{x} = \symbfit{y}\) y hallar finalmente \(\symbfit{x}\).
En R, usamos el paquete Matrix para obtener esta descomposición explícitamente.
Ejemplo 2.1 (Ejemplo: Factorización LU en R)
library(Matrix)
A <- matrix(c(2, 2, 1, 4, 3, 3, 8, 7, 7), nrow = 3, byrow = TRUE)
b <- c(1, 1, 1)
# Obtenemos la descomposición LU (con pivoteo)
decomp_lu <- lu(A)
print(decomp_lu)LU factorization of Formal class 'denseLU' [package "Matrix"] with 4 slots
..@ x : num [1:9] 8 0.5 0.25 7 -0.5 -0.5 7 -0.5 -1
..@ perm : int [1:3] 3 2 3
..@ Dim : int [1:2] 3 3
..@ Dimnames:List of 2
.. ..$ : NULL
.. ..$ : NULL# Extraemos las matrices L, U y P directamente del objeto
res <- Matrix::expand(decomp_lu)
L <- res$L
U <- res$U
P <- res$P # P es una matriz de permutación
# Verificación: P %*% A debe ser igual a L %*% U
print(P %*% A)3 x 3 Matrix of class "dgeMatrix"
[,1] [,2] [,3]
[1,] 8 7 7
[2,] 4 3 3
[3,] 2 2 1print(L %*% U)3 x 3 Matrix of class "dgeMatrix"
[,1] [,2] [,3]
[1,] 8 7 7
[2,] 4 3 3
[3,] 2 2 1# --- Resolución manual usando la factorización ---
# 1. Aplicamos la permutación a b
b_perm <- P %*% b
# 2. Resolvemos Ly = b_perm (sustitución progresiva)
# Usamos 'solve' que es inteligente y detecta que L es triangular
# o bien, como L es triangular inferior, podemos usar la función
# `forwardsolve()`
y <- solve(L, b_perm)
print(y)3 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 1.0
[2,] 0.5
[3,] 1.0# 3. Resolvemos Ux = y (sustitución regresiva)
# `solve()` o bien `backsolve()` son las funciones a recordar
x <- solve(U, y)
print(x)3 x 1 Matrix of class "dgeMatrix"
[,1]
[1,] 1
[2,] 0
[3,] -1# Comparar con la solución directa de R
solve(A, b)[1] 1 0 -12.2.2 Factorización de Cholesky
Un caso especial muy importante en estadística y machine learning (ej. en regresión de Mínimos Cuadrados Ordinarios) es cuando la matriz \(A\) es simétrica y definida positiva.
Recordemos que hemos visto en la asignatura de Estructuras algebraicas lo que se entiende por matriz (simétrica) y definida positiva. La definición es: decimos que una matriz \(A\in\mathcal M_n(\mathbb R)\) simétrica es definida positiva si \(\symbfit{x}^T A \symbfit{x} > 0\) para todo \(\symbfit{x}\neq\symbfit{0}\).
Tenemos varios criterios para comprobar la definición positiva de una matriz: si todos sus autovalores son estrictamente positivos, si todos los menores principales son estrictamente positivos, etc.
Definición 2.2 (Factorización de Cholesky) Si \(A\) es simétrica y definida positiva, existe una única matriz triangular inferior \(L\) tal que: \[A = LL^T\] Esto es computacionalmente más rápido (aproximadamente la mitad de operaciones que LU) y más estable numéricamente.
NotaAndré-Louis Cholesky
André-Louis Cholesky (1875-1918) fue un oficial militar y matemático francés. Desarrolló su famosa factorización para resolver problemas de topografía y geodesia mientras servía en el ejército francés. Su método no se publicó hasta después de su muerte en la Primera Guerra Mundial, gracias al esfuerzo de su colega el Comandante Benoit.
En R, la función chol() calcula la factorización de Cholesky.
Ejemplo 2.2 (Ejemplo: Cholesky en R)
# A es simétrica y definida positiva
A <- matrix(c(4, 2, -2, 2, 10, 2, -2, 2, 5), nrow = 3, byrow = TRUE)
b <- c(6, 15, 12)
# 1. Calcular la descomposición
R <- chol(A) # R es triangular SUPERIOR
print(R) [,1] [,2] [,3]
[1,] 2 1 -1.000000
[2,] 0 3 1.000000
[3,] 0 0 1.732051# Verificación: R^T * R debe ser A
t(R) %*% R [,1] [,2] [,3]
[1,] 4 2 -2
[2,] 2 10 2
[3,] -2 2 5# --- Resolución manual usando Cholesky ---
# A*x = b => (R^T * R) * x = b => R^T * (R*x) = b
# Paso 1: Resolver R^T * y = b (sustitución progresiva)
# (Usamos `forwardsolve()` en R base, o solve)
y <- forwardsolve(t(R), b)
# Paso 2: Resolver R * x = y (sustitución regresiva)
x <- backsolve(R, y)
print(x)[1] 3.2777778 0.1111111 3.6666667# Comparar con la solución directa
solve(A, b)[1] 3.2777778 0.1111111 3.66666672.2.2.1 Aplicación Estelar: Regresión Lineal por Mínimos Cuadrados
La factorización de Cholesky no es solo una curiosidad teórica; es el motor de cálculo de millones de modelos estadísticos cada día. Cuando ajustamos una regresión lineal \(\symbfit{y} = X\symbfit{\beta} + \symbfit{\epsilon}\), buscamos el vector de coeficientes \(\hat{\symbfit{\beta}}\) que minimiza el error cuadrático.
La solución analítica viene dada por las Ecuaciones Normales: \[(X^T X) \hat{\symbfit{\beta}} = X^T \symbfit{y}\]
La matriz \(A = X^T X\) es simétrica y (generalmente) definida positiva. Por tanto, para hallar \(\hat{\symbfit{\beta}}\), no invertimos \(A\). Hacemos Cholesky sobre \(A\).
# 1. Generamos datos sintéticos: y = 2 + 1.5*x + ruido
set.seed(123)
n <- 100
x <- rnorm(n)
y <- 2 + 1.5 * x + rnorm(n, sd = 0.5)
# 2. Construimos la Matriz de Diseño X (incluyendo la columna de 1s para el intercepto)
X <- cbind(Intercept = 1, Slope = x)
# 3. Calculamos la matriz A = X^T * X y el vector b = X^T * y
A <- t(X) %*% X
b <- t(X) %*% y
# 4. Resolvemos el sistema A * beta = b usando Cholesky
# A = R^T * R
R <- chol(A) # R es triangular superior
# R^T * y_temp = b -> solve(t(R), b)
y_temp <- solve(t(R), b)
# R * beta = y_temp -> solve(R, y_temp)
beta_hat <- solve(R, y_temp)
print("Coeficientes calculados manualmente con Cholesky:")[1] "Coeficientes calculados manualmente con Cholesky:"print(t(beta_hat)) Intercept Slope
[1,] 1.948598 1.473764# 5. Comparación con la función lm() de R (la forma estándar)
modelo <- lm(y ~ x)
print("Coeficientes con lm():")[1] "Coeficientes con lm():"print(coef(modelo))(Intercept) x
1.948598 1.473764 ¡Exactamente el mismo resultado! Así es como funcionan las estadísticas “bajo el capó”.
2.2.3 Factorización QR
Definición 2.3 (Factorización QR) Para una matriz \(A\) (no necesariamente cuadrada), la factorización QR encuentra una matriz ortogonal \(Q\) y una matriz triangular superior \(R\) tal que: \[A = QR\]
Recordemos que una matriz \(Q\) es ortogonal si sus columnas forman un sistema de vectores ortonormales, es decir, si \(Q^T Q = I\).
Existen muchos métodos para la construcción de la matriz \(Q\), pero desde el punto de vista algebraico, podemos recurrir al método de Gram-Schmidt que se vio en la asignatura de Estructuras Algebraicas.
¿Cómo se resuelve el sistema \(A\symbfit{x} = \symbfit{b}\) con QR?
Si \(A = QR\), entonces el sistema \(A\symbfit{x} = \symbfit{b}\) se convierte en \(QR\symbfit{x} = \symbfit{b}\). Multiplicando por \(Q^T\) en ambos lados y aprovechando que \(Q^T Q = I\):
\[R\symbfit{x} = Q^T \symbfit{b}\]
Como \(R\) es una matriz triangular superior, este sistema se resuelve de forma muy eficiente mediante sustitución regresiva.
A <- matrix(c(4, 1, -1, 1, 6, 2, -1, 2, 5), nrow = 3, byrow = TRUE)
b <- c(6, 15, 12)
# 1. Calcular la descomposición QR
decomp_qr <- qr(A)
# 2. Resolver directamente
# (método preferido por estabilidad numérica)
x_qr <- qr.solve(decomp_qr, b)
print(x_qr)[1] 1.651685 1.516854 2.123596# --- Resolución manual paso a paso ---
## Éstas son las funciones para obtener Q y R
Q <- qr.Q(decomp_qr)
R <- qr.R(decomp_qr)
# Paso 1: y = Q^T * b
y <- t(Q) %*% b
# Paso 2: R * x = y (sustitución regresiva)
x_manual <- backsolve(R, y)
print(x_manual) [,1]
[1,] 1.651685
[2,] 1.516854
[3,] 2.1235962.3 Análisis detallado de matrices con estructura
Como mencionamos al inicio, el principal cuello de botella en sistemas grandes es la memoria. La solución es no almacenar los ceros. En el mundo real, la mayoría de los sistemas masivos que provienen de la discretización de fenómenos físicos (calor, fluidos, elasticidad) o del análisis de redes (sociales, de internet) son dispersos.
Definición 2.4 (Matriz Dispersa (Sparse)) Una matriz dispersa es una matriz en la que la mayoría de sus elementos son cero. Una matriz es “suficientemente dispersa” si vale la pena utilizar algoritmos y estructuras de datos especiales para evitar almacenar y operar con esos ceros.
El nivel de dispersión (o sparsity) se mide como: \[\text{Sparsity} = 1 - \frac{\text{Nº elementos no nulos (NNZ)}}{\text{Nº total de elementos (} n \times m \text{)}}\] Una matriz con un sparsity del 99.9% es común.
2.3.1 El paquete Matrix: el ecosistema disperso
R base no maneja matrices dispersas de forma nativa. Un matrix en R base siempre almacena todos sus elementos. El paquete Matrix es la herramienta estándar y fundamental para el álgebra lineal numérica en R.
Ejemplo 2.3 (Ejemplo: Creando Matrices Dispersas) Comparemos una matriz densa vs. una dispersa de \(5000 \times 5000\).
library(Matrix)
n <- 5000
# 1. Matriz densa (R base)
A_densa <- matrix(0, nrow = n, ncol = n)
A_densa[1, 1] <- 5
A_densa[n / 2, n / 2] <- 5
A_densa[n, n] <- 5
# 2. Matriz dispersa (Paquete Matrix)
# Creamos una matriz especificando solo las posiciones (i, j) y los valores (x)
A_dispersa <- sparseMatrix(
i = c(1, n / 2, n),
j = c(1, n / 2, n),
x = 5,
dims = c(n, n)
)
# Comparamos el tamaño en memoria
print(object.size(A_densa), units = "MB")190.7 Mbprint(object.size(A_dispersa), units = "KB")21 Kb¡La matriz dispersa ocupa solo 21 Kb frente a los 190.7 Mb de la densa! Esto es una diferencia de más de 9000 veces en el uso de memoria para almacenar la misma información.
2.3.2 Optimizaciones para matrices triangulares
El primer tipo de matriz estructurada que encontramos son las triangulares, que son la base de las factorizaciones LU y Cholesky.
Resolver un sistema triangular \(T\symbfit{x} = \symbfit{b}\) no requiere un método \(O(n^3)\) como Gauss. Se puede hacer directamente por sustitución (regresiva si es triangular superior, progresiva si es inferior) en \(O(n^2)\).
Las funciones backsolve() y forwardsolve() de R base están optimizadas para esto.
Ejemplo 2.4 (Ejemplo: Eficiencia de Solvers Triangulares) Comparemos el tiempo de un solver genérico (solve) contra uno especializado (backsolve) para un sistema triangular superior.
library(microbenchmark)
library(Matrix)
n <- 1000 # Un sistema de 1000x1000
# 1. Creamos un sistema triangular superior denso
R_densa <- matrix(0.1 * runif(n * n), n, n)
diag(R_densa) <- 1
R_densa[lower.tri(R_densa)] <- 0 # Hacemos que sea triangular superior
b <- rnorm(n)
# Convertimos R_densa a una matriz dispersa triangular
R_sparse <- as(R_densa, "dgCMatrix")
# 2. Comparamos los tiempos de resolución
mbm <- microbenchmark(
generico_denso = solve(R_densa, b),
especializado_denso = backsolve(R_densa, b),
generico_disperso = solve(R_sparse, b),
times = 10
)
print(mbm, unit = "ms") # mostramos en milisegundosUnit: milliseconds
expr min lq mean median uq
generico_denso 97.857037 99.328650 110.0334630 104.5604345 119.339848
especializado_denso 0.240301 0.256824 0.3180944 0.3146135 0.365966
generico_disperso 0.258915 0.331731 1.7539103 0.3619275 0.386876
max neval
132.761854 10
0.402538 10
14.220358 10De los resultados, observamos que backsolve es significativamente más rápido que solve en matrices densas de R base. Pero aún más importante, el método solve del paquete Matrix es también extremadamente rápido, ya que detecta la estructura triangular y aplica el solver más eficiente posible (sustitución progresiva/regresiva) sobre la estructura de datos dispersa.
2.3.3 Matrices banda y tridiagonales
El caso más importante de matrices estructuradas en la práctica son las matrices banda.
Definición 2.5 (Matriz Banda) Una matriz \(A\) es una matriz banda si todos sus elementos no nulos se encuentran en una banda alrededor de la diagonal principal. Formalmente, \(a_{ij} = 0\) si \(|i-j| > k\), donde \(k\) es el semiancho de banda (bandwidth).
Una matriz \(A\) es una matriz banda si todos sus elementos no nulos se encuentran confinados en una franja alrededor de la diagonal principal. Dependiendo del ancho de dicha franja (\(k\)), recibimos nombres específicos: si \(k=0\) estamos ante una matriz diagonal, mientras que \(k=1\) define una matriz tridiagonal y \(k=2\) una pentadiagonal.
Las matrices tridiagonales son ubicuas en la modelización matemática, apareciendo frecuentemente en la discretización de ecuaciones diferenciales (como la Ecuación del Calor), en el cálculo de splines para interpolación estadística y en diversos modelos de dinámica de poblaciones.
El Algoritmo de Thomas Para los sistemas tridiagonales \(A\symbfit{x} = \symbfit{b}\), no es necesario usar Eliminación de Gauss (\(O(n^3)\)) ni siquiera un solver de banda general. Existe un algoritmo específico, llamado Algoritmo de Thomas, que es una versión simplificada de la eliminación de Gauss que resuelve el sistema en coste de tiempo y memoria \(O(n)\).
NotaLlewellyn Thomas
Llewellyn Thomas (1903-1992) fue un físico y matemático británico-estadounidense. Es famoso por la “precesión de Thomas” en mecánica cuántica. Desarrolló este algoritmo en 1949 para resolver problemas de física de fluidos. En la Unión Soviética, se conocía como el “método de progonka”.
Cuando usamos solve(A, b) sobre una matriz tridiagonal del paquete Matrix, éste detecta la estructura de banda y aplica un solver altamente optimizado (como una versión del Algoritmo de Thomas o una factorización LU específica para bandas) que es, en efecto, \(O(n)\).
Ejemplo 2.5 (Creación y Resolución de un Sistema Tridiagonal) Vamos a crear una matriz tridiagonal (con \(n=5\)) y a resolver un sistema de ecuaciones con ella.
\(A = \begin{pmatrix} 2 & -1 & 0 & 0 & 0 \\\ -1 & 2 & -1 & 0 & 0 \\\ 0 & -1 & 2 & -1 & 0 \\\ 0 & 0 & -1 & 2 & -1 \\\ 0 & 0 & 0 & -1 & 2 \end{pmatrix}, \quad \symbfit{b} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 1 \end{pmatrix}\)
n <- 5
# 1. Definir las diagonales
diag_principal <- rep(2, n)
sub_diag <- rep(-1, n - 1)
super_diag <- rep(-1, n - 1)
# 2. Definir las posiciones (i, j)
i_idx <- c(1:n, 2:n, 1:(n - 1))
j_idx <- c(1:n, 1:(n - 1), 2:n)
x_val <- c(diag_principal, sub_diag, super_diag)
# 3. Crear la matriz dispersa
A_tridiag <- sparseMatrix(
i = i_idx,
j = j_idx,
x = x_val,
dims = c(n, n)
)
# 4. Crear el vector b
b <- rep(0, n)
b[n] <- 1
print(A_tridiag)5 x 5 sparse Matrix of class "dgCMatrix"
[1,] 2 -1 . . .
[2,] -1 2 -1 . .
[3,] . -1 2 -1 .
[4,] . . -1 2 -1
[5,] . . . -1 2print(b)[1] 0 0 0 0 1# 5. Resolver
x <- solve(A_tridiag, b)
print(x)[1] 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333Ejemplo 2.6 (Comparativa de Rendimiento: Densa vs. Dispersa) Vamos a comparar la resolución de un sistema tridiagonal de \(n=5000\) usando una matriz densa (R base) frente a una matriz dispersa (Matrix).
n <- 5000
# --- 1. Matriz Densa ---
A_densa <- matrix(0, n, n)
diag(A_densa) <- 2
# Rellena las sub/super-diagonales
# (¡muy lento en R!)
# La notación `[` no permite las asignaciones
# en bloques como nos gustaría.
for (i in 1:(n - 1)) {
A_densa[i, i + 1] <- -1
A_densa[i + 1, i] <- -1
}
b <- rep(0, n)
b[n] <- 100
# --- 2. Matriz Dispersa ---
i_idx <- c(1:n, 2:n, 1:(n - 1))
j_idx <- c(1:n, 1:(n - 1), 2:n)
x_val <- c(rep(2, n), rep(-1, (n - 1) * 2))
A_dispersa <- sparseMatrix(
i = i_idx,
j = j_idx,
x = x_val,
dims = c(n, n)
)
# --- 3. Comparativa de Tiempos ---
# NOTA: Ejecutamos solo 1 vez porque es MUY lento.
print(system.time({
x_densa <- solve(A_densa, b)
})) user system elapsed
12.030 0.107 12.398 print(system.time({
x_dispersa <- solve(A_dispersa, b)
})) user system elapsed
0.001 0.000 0.001 El resultado es dramático. En una máquina estándar:
- La solución densa puede tardar más de 10 segundos (y consumir ~200 MB de RAM).
- La solución dispersa tarda menos de 0.01 segundos (y consumir ~150 KB de RAM).
Para la modelización predictiva y los sistemas complejos, la diferencia no es solo de velocidad, sino de viabilidad. Un problema de \(n=50.000\) sería imposible con matrices densas, pero se resolvería en menos de un segundo con matrices dispersas.
Ejercicio 2.1
- Adapta el código del Ejemplo 2.5 para crear una matriz \(A\) de \(100 \times 100\) que sea pentadiagonal.
- La diagonal principal debe tener el valor ‘10’.
- Las primeras sub/super-diagonales deben tener el valor ‘-2’.
- Las segundas sub/super-diagonales (es decir, \(a_{i, i+2}\) y \(a_{i, i-2}\)) deben tener el valor ‘1’.
- Visualiza la estructura de la matriz resultante con
image(A).
Ejercicio 2.2
- Crea una matriz triangular superior \(R\) de \(1000\times1000\) con números aleatorios en y sobre la diagonal. (Pista: para generar números aleatorios, puedes usar la función
rnormorunif, y a continuación puedes usar la funciónmatrixpara convertirlos en una matriz. A continuación, puedes usar las funcionesupper.triylower.trisabiamente para hacer 0 por debajo de la diagonal principal). - Crea un vector \(\symbfit{b}\) de 1000 números aleatorios.
- Resuelve el sistema \(R\symbfit{x} = \symbfit{b}\) de dos maneras:
- Usando el solver genérico (
solve()). - Usando el solver específico para sistemas triangulares (
backsolve()).
- Usando el solver genérico (
- Usa el paquete
microbenchmarkpara comparar el tiempo de ambas operaciones:microbenchmark(solve(R, b), backsolve(R, b), times = 10). - ¿Qué conclusiones extraes sobre usar la función adecuada a la estructura de la matriz?
2.4 Estabilidad numérica y condicionamiento
Hasta ahora, hemos asumido que la solución que nos da el ordenador es fiable. Sin embargo, en el mundo de la aritmética de punto flotante, los errores de redondeo pueden acumularse y destruir la precisión si el sistema es inestable.
Definición 2.6 (Número de Condición) El número de condición de una matriz, denotado como \(\kappa(A) = ||A|| \cdot ||A^{-1}||\), mide cuánto se amplifica un pequeño cambio en los datos de entrada (\(\symbfit{b}\)) en la solución (\(\symbfit{x}\)).
- Si \(\kappa(A) \approx 1\), el sistema es bien condicionado.
- Si \(\kappa(A)\) es muy grande (ej. \(10^{16}\)), el sistema es mal condicionado o singular.
En R, podemos estimar el recíproco del número de condición con rcond() o calcularlo (para matrices pequeñas) con kappa().
A <- matrix(c(1, 2, 3, 4), nrow = 2)
# Uno es el inverso del otro
cat("Recíproco del número de condición (rcond):\n",
rcond(A), "\n")Recíproco del número de condición (rcond):
0.04761905 cat("Número de condición (kappa):\n",
kappa(A), "\n")Número de condición (kappa):
18.77778 Un número de condición elevado advierte que la solución es numéricamente inestable: pequeñas variaciones en los coeficientes o en el vector \(\symbfit{b}\) (debidas a ruido en los datos o a la precisión limitada del punto flotante) provocan cambios desproporcionados en \(\symbfit{x}\). Como regla general, si \(\kappa(A) = 10^k\), se pueden perder hasta \(k\) dígitos de precisión.
Ejemplo 2.7 (Ejemplo de Mal Condicionamiento) Veamos un ejemplo de sistema mal condicionado. El sistema que vamos a estudiar es \(Ax = b\), donde: \[ A = \begin{pmatrix} 1 & 1 \\ 1 & 1.0001 \end{pmatrix}, \quad b = \begin{pmatrix} 2 \\ 2.0001 \end{pmatrix} \] es decir, el sistema: \[\left\{ \begin{array}{rcrcl} x_1 & + & x_2 & = & 2 \\ x_1 & + & 1.001x_2 & = & 2.0001 \end{array} \right.\] La solución original es \(x = (1, 1)\). Veamos lo que ocurre si alteramos \(b\) mínimamente.
# Definimos A y b
A <- matrix(c(1, 1, 1, 1.0001), nrow = 2)
b <- c(2, 2.0001)
# Solución original: x = (1, 1)
x <- solve(A, b)
cat("Solución original:\n ", x, "\n")Solución original:
1 1 # Alteramos b mínimamente (0.0001 de diferencia)
b_alt <- c(2, 2.0002)
x_alt <- solve(A, b_alt)
cat("Solución con b ligeramente alterado:\n ", x_alt, "\n")Solución con b ligeramente alterado:
0 2 # Observamos el mal condicionamiento
cat("Número de condición (kappa):\n ",
kappa(A, exact = TRUE), "\n")Número de condición (kappa):
40002 En este ejemplo, un cambio de \(10^{-4}\) en un elemento de \(\symbfit{b}\) ha cambiado la solución de \((1, 1)\) a \((0, 2)\), lo que demuestra que el sistema no es fiable para cálculos de alta precisión.
2.5 Métodos iterativos para sistemas grandes
Cuando una matriz \(A\) es tan masiva que ni siquiera la factorización LU/Cholesky es viable (pensemos en \(n=10^8\)), o cuando la matriz es densa pero muy grande, los métodos directos fallan.
La alternativa son los métodos iterativos. Estos métodos no dan la solución exacta, sino que construyen una sucesión de vectores \(\symbfit{x}^{(k)}\) que converge a la solución real \(\symbfit{x}\).
Definición 2.7 (Método Iterativo) Un método iterativo reescribe \(A\symbfit{x} = \symbfit{b}\) en una forma de punto fijo equivalente: \(\symbfit{x} = T\symbfit{x} + \symbfit{c}\).
Comenzando con una estimación inicial \(\symbfit{x}^{(0)}\) (a menudo \(\symbfit{0}\)), se genera la sucesión:
\[\symbfit{x}^{(k+1)} = T\symbfit{x}^{(k)} + \symbfit{c}\]
Si el método converge, \(\lim_{k \to \infty} \symbfit{x}^{(k)} = \symbfit{x}\).
Para descomponer \(A\), usamos la notación \(A = D + L + R\), donde:
- \(D\) es la matriz diagonal de \(A\).
- \(L\) es la matriz triangular estrictamente inferior de \(A\).
- \(R\) es la matriz triangular estrictamente superior de \(A\).
2.5.1 Método de Jacobi
El método de Jacobi se deriva de reescribir \(A\symbfit{x} = \symbfit{b}\) como \((D+L+R)\symbfit{x} = \symbfit{b}\), y despejando la \(x\) de la parte diagonal:
\(D\symbfit{x} = -(L+R)\symbfit{x} + \symbfit{b}\)
\(\symbfit{x} = -D^{-1}(L+R)\symbfit{x} + D^{-1}\symbfit{b}\)
Esto nos da la fórmula iterativa de Jacobi:
\[\symbfit{x}^{(k+1)} = -D^{-1}(L+R)\symbfit{x}^{(k)} + D^{-1}\symbfit{b}\]
En la práctica, esto significa que para calcular cada componente \(x_i^{(k+1)}\), usamos solo los valores de la iteración anterior \(\symbfit{x}^{(k)}\).
Esto es equivalente a escribir este método iterativo: \[\left. \begin{array}{ll} x_{1}^{(k+1)} = &\frac{1}{a_{1,1}}\left( b_{1} -a_{1,2}x_{2}^{(k)} -a_{1,3}x_{3}^{(k)} - \ldots -a_{1,n}x_{n}^{(k)}\right) \\ x_{2}^{(k+1)} = &\frac{1}{a_{2,2}}\left( b_{2} -a_{2,1}x_{1}^{(k)} -a_{2,3}x_{3}^{(k)} - \ldots -a_{2,n}x_{n}^{(k)}\right) \\ \vdots & \vdots \\ x_{i}^{(k+1)} = &\frac{1}{a_{i,i}}\left( b_{i} -a_{i,1}x_{1}^{(k)} -\ldots -a_{i,i-1}x_{i-1}^{(k)}-a_{i,i+1}x_{i+1}^{(k)}-\ldots-a_{i,n}x_{n}^{(k)}\right)\\ \vdots&\vdots\\ x_{n}^{(k+1)} = &\frac{1}{a_{n,n}}\left( b_{n} -a_{n,1}x_{1}^{(k)} -a_{n,2}x_{2}^{(k)} - \ldots -a_{n,n-1}x_{n-1}^{(k)}\right) \\ \end{array} \right\}\]
NotaCarl Jacobi
Carl Gustav Jacob Jacobi (1804-1851) fue un matemático alemán, considerado uno de los más grandes de su tiempo. Aunque es famoso por sus trabajos en funciones elípticas y dinámica, el método iterativo que lleva su nombre es hoy esencial para resolver sistemas masivos de forma paralela, ya que cada componente puede calcularse de forma independiente.
El algoritmo de Jacobi es simple y fijo, pero puede converger lentamente para matrices con valores propios muy desiguales.
Ejemplo 2.8 (Uso del método de Jacobi en R) El método de Jacobi es simple de implementar, pero puede ser ineficiente para matrices grandes. En R, existen paquetes como pracma o Rlinsolve que implementan métodos iterativos eficientes.
library(Rlinsolve)** ------------------------------------------------------- **** Rlinsolve** - Solving (Sparse) System of Linear Equations**** Version : 0.3.3 (2026)** Maintainer : Kisung You (kisung.you@outlook.com)**** Please share any bugs or suggestions to the maintainer.** ------------------------------------------------------- **# Sistema de un ejemplo anterior
A <- matrix(c(4, 1, -1, 1, 6, 2, -1, 2, 5), nrow = 3, byrow = TRUE)
b <- c(6, 15, 12)
x0 <- c(0, 0, 0) # Empezamos desde cero
lsolve.jacobi(A, b, xinit = x0)* lsolve.jacobi : Initialiszed.* lsolve.jacobi : computations finished.$x
[,1]
[1,] 1.651648
[2,] 1.516892
[3,] 2.123553
$iter
[1] 25
$errors
[,1]
[1,] 1.329212e-01
[2,] 5.762669e-02
[3,] 3.927065e-02
[4,] 2.728179e-02
[5,] 1.889019e-02
[6,] 1.305733e-02
[7,] 9.020086e-03
[8,] 6.229908e-03
[9,] 4.302543e-03
[10,] 2.971393e-03
[11,] 2.052070e-03
[12,] 1.417174e-03
[13,] 9.787103e-04
[14,] 6.759039e-04
[15,] 4.667838e-04
[16,] 3.223640e-04
[17,] 2.226267e-04
[18,] 1.537475e-04
[19,] 1.061790e-04
[20,] 7.332795e-05
[21,] 5.064077e-05
[22,] 3.497285e-05
[23,] 2.415248e-05
[24,] 1.667987e-05
[25,] 1.151923e-05
[26,] 7.955254e-06Ejercicio 2.3 (Método de Jacobi) Resolver el siguiente sistema utilizando el método de Jacobi, dando 300 iteraciones como máximo y partiendo del vector nulo: \[ \left\{ \begin{array}{cccccc} 10x & -y & +2z & & = & 6 \\ -x & +11y & -z & +3t & = & 25 \\ 2x & -y &+10z & -t & = & -11 \\ & 3y & -z & +8t & = & 15 \end{array} \right. \]
2.5.2 Método de Gauss-Seidel
El método de Gauss-Seidel es una optimización simple de Jacobi. Cuando estamos calculando \(x_i^{(k+1)}\), ya hemos calculado \(x_1^{(k+1)}, \dots, x_{i-1}^{(k+1)}\) en la misma iteración. ¿Por qué no usar esos valores nuevos y supuestamente mejores en lugar de los antiguos de \(x^{(k)}\)?
La fórmula se deriva de \((D+L)\symbfit{x} = -R\symbfit{x} + \symbfit{b}\):
\[\symbfit{x}^{(k+1)} = -(D+L)^{-1}R\symbfit{x}^{(k)} + (D+L)^{-1}\symbfit{b}\]
NotaPhilipp Ludwig von Seidel
Philipp Ludwig von Seidel (1821-1896) fue un matemático y astrónomo alemán. Su refinamiento del método de Jacobi, conocido hoy como Gauss-Seidel, acelera la convergencia al usar los valores más recientes tan pronto como están disponibles. Seidel también hizo contribuciones fundamentales a la óptica y al análisis de errores en astronomía.
O, más explícitamente: \[\left. \begin{array}{cc} x_{1}^{(k+1)} = &\frac{1}{a_{1,1}}\left( b_{1} -a_{1,2}x_{2}^{(k)} -a_{1,3}x_{3}^{(k)} - \ldots -a_{1,n}x_{n}^{(k)}\right) \\ x_{2}^{(k+1)} = &\frac{1}{a_{2,2}}\left( b_{2} -a_{2,1}\textcolor{blue}{x_{1}^{(k+1)}} -a_{2,3}x_{3}^{(k)} - \ldots -a_{2,n}x_{n}^{(k)}\right) \\ \vdots & \vdots \\ x_{i}^{(k+1)} = &\frac{1}{a_{i,i}}\left( b_{i} -a_{i,1}\textcolor{blue}{x_{1}^{(k+1)}} -\ldots -a_{i,i-1}\textcolor{blue}{x_{i-1}^{(k+1)}}-a_{i,i+1}x_{i+1}^{(k)}-\ldots-a_{i,n}x_{n}^{(k)}\right)\\ \vdots&\vdots\\ x_{n}^{(k+1)} = &\frac{1}{a_{n,n}}\left( b_{n} -a_{n,1}\textcolor{blue}{x_{1}^{(k+1)}} -a_{n,2}\textcolor{blue}{x_{2}^{(k+1)}} - \ldots -a_{n,n-1}\textcolor{blue}{x_{n-1}^{(k+1)}}\right) \\ \end{array} \right\}\]
Ejemplo 2.9 (Uso del método de Gauss-Seidel en R) Vamos a utilizar el método lsolve.gs de la librería Rlinsolve para resolver el sistema anterior.
# Nota: por defecto, si no le especificamos
# el vector inicial, lo inicializa
# cercano al vector cero,
# de forma aleatoria
lsolve.gs(A, b)* lsolve.gs : Initialiszed.* lsolve.gs : computations finished.$x
[,1]
[1,] 1.651642
[2,] 1.516885
[3,] 2.123574
$iter
[1] 8
$errors
[,1]
[1,] 1.501422e-01
[2,] 3.297790e-02
[3,] 8.791833e-03
[4,] 2.656862e-03
[5,] 7.884058e-04
[6,] 2.348104e-04
[7,] 6.988591e-05
[8,] 2.080260e-05
[9,] 6.192061e-06Generalmente, Gauss-Seidel converge más rápido que Jacobi, podemos ver el número de iteraciones que ha realizado cada método hasta llegar a la aproximación final.
2.5.3 Método del gradiente conjugado (CG)
Para matrices simétricas y definidas positivas (muy comunes en física y optimización), el método del Gradiente Conjugado es generalmente superior a Jacobi y Gauss-Seidel.
A diferencia de los anteriores, que “suavizan” el error localmente, CG busca la solución a lo largo de direcciones “conjugadas” (ortogonales respecto a la métrica inducida por \(A\)). Teóricamente, en aritmética exacta, converge en \(n\) pasos. En aritmética flotante, se utiliza como método iterativo y suele converger mucho antes para una tolerancia dada.
Ejemplo 2.10 (Uso del método de Gradiente Conjugado en R) Vamos a utilizar el método lsolve.cg de la librería Rlinsolve para resolver el sistema anterior.
lsolve.cg(A, b)* lsolve.cg : Initialiszed.* lsolve.cg : preprocessing finished ...* lsolve.cg : convergence was well achieved.* lsolve.cg : computations finished.$x
[,1]
[1,] 1.651685
[2,] 1.516854
[3,] 2.123596
$iter
[1] 2
$errors
[,1]
[1,] 0.81756550
[2,] 0.22477011
[3,] 0.072404372.5.4 Convergencia de métodos iterativos
No todos los sistemas convergen con estos métodos. Una condición suficiente (pero no necesaria) muy útil es la de diagonal dominante.
Definición 2.8 (Matriz de Diagonal (Estrictamente) Dominante) Una matriz \(A\) es de diagonal dominante (por filas) si para cada fila \(i\), el valor absoluto del elemento de la diagonal es mayor o igual que la suma de los valores absolutos de los demás elementos de esa fila.
\[|a_{ii}| \geq \sum_{j \neq i} |a_{ij}| \quad \text{para todo } i\]
Si la desigualdad es estricta, es decir, \(|a_{ii}| > \sum_{j \neq i} |a_{ij}|\), se dice que la matriz es de diagonal estrictamente dominante.
Teorema 2.1 (Convergencia de métodos iterativos)
- Si \(A\) es de diagonal dominante, el método de Jacobi converge a la única solución del sistema.
- Si \(A\) es de diagonal estrictamente dominante, o si \(A\) es simétrica y definida positiva, el método de Gauss-Seidel y el método del gradiente conjugado convergen a la única solución del sistema.
En todos los casos, el método converge sin importar el vector inicial \(\symbfit{x}^{(0)}\).
2.6 Práctica: Simulación de Red de Calor
Objetivo: Modelar un problema físico simple que genera un sistema de ecuaciones grande y disperso, y comparar la eficiencia de los métodos de resolución.
Problema: Imaginemos una barra de metal de 1 metro. Mantenemos el extremo izquierdo a \(0^\circ C\) y el extremo derecho a \(100^\circ C\). Queremos encontrar la temperatura en equilibrio en \(n\) puntos interiores de la barra.
La física (usando la ecuación del calor en estado estacionario) nos dice que la temperatura en un punto \(i\) es el promedio de la temperatura de sus vecinos \(i-1\) e \(i+1\).
\[T_i = \frac{T_{i-1} + T_{i+1}}{2} \implies -T_{i-1} + 2T_i - T_{i+1} = 0\]
Si discretizáramos la barra en \(n=10\) puntos interiores (\(T_1, \dots, T_{10}\)), con los bordes \(T_0 = 0\) y \(T_{11} = 100\), generamos un sistema de ecuaciones:
- \(i=1\): Despejamos y aplicamos la condición de contorno: \[-T_0 + 2T_1 - T_2 = 0 \implies 2T_1 - T_2 = T_0 \implies 2T_1 - T_2 = 0\]
- \(i=2\): \(-T_1 + 2T_2 - T_3 = 0\)
- …
- \(i=10\): Igual que en el primer caso, despejamos y usamos la condición de contorno \(T_{11} = 100\). \[-T_9 + 2T_{10} - T_{11} = 0 \implies -T_9 + 2T_{10} = T_{11} \implies -T_9 + 2T_{10} = 100\]
Esto genera una matriz \(A\) de \(10 \times 10\) tridiagonal y un vector \(b\):
\(A = \begin{pmatrix} 2 & -1 & 0 & \dots \\ -1 & 2 & -1 & \dots \\ 0 & -1 & 2 & \dots \\ \vdots & \vdots & \ddots & -1 \\ \dots & 0 & -1 & 2 \end{pmatrix}, \quad \symbfit{b} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 100 \end{pmatrix}\)
Tu Tarea:
Crear la Matriz y el Vector (n=1000): No lo hagas con \(n=10\), ¡hazlo con \(n=1000\)! Usa
sparseMatrix()para crear la matriz \(A\) (\(1000\times1000\)) y el vector \(b\) (\(1000\times1\)).Resolver y Comparar Tiempos: Usa
microbenchmark::microbenchmark()para comparar el tiempo de ejecución de los siguientes métodos para resolver \(A\symbfit{x} = \symbfit{b}\) (con \(n=1000\)):- El “solver” genérico de R.
- Realización de la factorización de Cholesky.
- Factorización QR.
- Método de Jacobi.
- Método de Gauss-Seidel.
- Método de Gradiente Conjugado.
Analizar los Resultados:
- Escribe una breve conclusión. ¿Es \(A\) de diagonal dominante? ¿Por qué convergen los métodos iterativos?
- ¿Cuál fue el método más rápido?
- ¿Qué crees que pasaría si \(n = 100.000\) (no lo ejecutes, solo razónalo)? ¿Qué método(s) fallaría(n) por problemas de memoria al usar
as.matrix()?
2.7 Ejercicios Adicionales
Ejercicio 2.4
Usando sparseMatrix(), construye una matriz \(A\) de \(100 \times 100\) que tenga la diagonal principal llena de ‘5’, las segundas subdiagonales (k=-2) y superdiagonales (k=2) llenas de ‘1’, y el resto ceros. Visualízala con image(A).
Ejercicio 2.5
La matriz \(A\) del Ejemplo 2.2 es simétrica y definida positiva. Resuelve el sistema \(A^2 \symbfit{x} = \symbfit{b}\) (es decir, \(A(A\symbfit{x}) = \symbfit{b}\)).
Pista: No calcules \(A^2\). Usa la factorización \(A=R^T R\) y resuelve \(R^T(R(R^T(R \symbfit{x}))) = \symbfit{b}\) en cuatro pasos de sustitución progresiva/regresiva.