La PizarraMatemática.

Máster y doctorado

Precondicionamiento y métodos de Krylov

En un sistema enorme suele ser fácil multiplicar una matriz por un vector y difícil guardar su inversa. Los métodos de Krylov construyen la solución a partir de esas multiplicaciones. Resolveremos un sistema de dos componentes a mano para ver qué significan sus pasos y por qué precondicionar puede ahorrar trabajo.

1. El residuo comprueba las ecuaciones y la condición relaciona ese fallo con el error

Para Ax=b, una aproximación x deja residuo r=b−Ax.
Si A es invertible y x* la solución, el error e=x*−x cumple Ae=r.
Por tanto ∥e∥≤∥A⁻¹∥∥r∥.

El tamaño del residuo no dice por sí solo el tamaño del error.

Con b≠0, la desigualdad relativa es ∥e∥/∥x*∥≤κ(A)∥r∥/∥b∥ en normas compatibles.

Una κ grande puede hacer poco informativa la cota, aunque el residuo sea pequeño. El algoritmo usa productos Av y una regla de parada que fija norma, escalas y tolerancia.

Un residuo pequeño y un error grande

A=diag(1,1/1000), b=(1,1/1000), aproximación x=(1,0).

  1. Solución

    x*=(1,1).
  2. Residuo

    b−Ax=(0,1/1000), norma 0,001.
  3. Error

    x*−x=(0,1), norma uno.
Resultado

La dirección débil puede ocultar un error importante.

Lo comprobamos

Aquí la inversa multiplica esa dirección por mil.

Un producto no requiere inversa

A=[[2,−1],[−1,2]] y v=(3,1).

  1. Primera fila

    2·3−1=5.
  2. Segunda

    −3+2·1=−1.
  3. Resultado

    Av=(5,−1).
Resultado

Cada producto usa las conexiones de la matriz.

Lo comprobamos

En una matriz dispersa se aprovechan sus entradas no nulas.

2. El espacio de Krylov reúne las direcciones generadas por la matriz

Desde x₀, calculamos r₀=b−Ax₀.

El espacio \mathcal K_m(A,r_0)=\operatorname{span}\{r_0,Ar_0,\ldots,A^{m-1}r_0\} contiene todas sus combinaciones lineales. La búsqueda se hace en x₀+K_{m}, no solo en uno de los vectores de la lista.

Hay m generadores, pero su dimensión puede ser menor.

Si r₀=0 ya se ha resuelto.

Si los productos generan un subespacio invariante, un método apropiado puede terminar antes de la dimensión total. Dependencia exacta y pérdida numérica de independencia se tratan de forma distinta.

Construir dos direcciones

A=diag(1,4), x₀=0, b=(1,1).

  1. Residuo

    r₀=(1,1).
  2. Producto

    Ar₀=(1,4).
  3. Independencia

    El determinante de las dos columnas es tres; K₂=ℝ².
Resultado

Dos direcciones bastan para representar cualquier candidato en este ejemplo.

Lo comprobamos

Falta elegir cuál es el mejor en el espacio.

Cuatro generadores no siempre dan dimensión cuatro

A=2I y r₀≠0.

  1. Productos

    Ar₀=2r₀, A²r₀=4r₀ y A³r₀=8r₀.
  2. Span

    Todos están en la misma recta.

  3. Dimensión

    K₄ tiene dimensión uno, con cuatro vectores en su lista generadora.

Resultado

Contar la lista y calcular el rango son tareas diferentes.

Lo comprobamos

En este caso el problema se resuelve dentro de esa recta.

3. CG minimiza una energía cuando A es simétrica positiva definida

A SPD significa A=Aᵀ y vᵀAv>0 para v≠0.
La función Φ(x)=½xᵀAx−bᵀx tiene gradiente Ax−b y mínimo único x*.
Además \Phi(x)-\Phi(x^*)=\tfrac12(x-x^*)^TA(x-x^*).

CG elige x_{m} que minimiza esa norma de error A sobre x₀+K_{m}.

Sus direcciones son conjugadas: p_{i}ᵀAp_j=0 si i≠j.

No tienen por qué ser ortogonales euclídeamente. La propiedad se garantiza aquí para SPD y aritmética exacta.

Comprobar positividad

A=diag(1,4).

  1. Forma

    vᵀAv=v₁²+4v₂².
  2. Signo

    Es positiva para cualquier vector no nulo.

  3. Energía

    ∥v∥A=√(v₁²+4v₂²).
Resultado

La segunda componente pesa más en esta norma.

Lo comprobamos

No se usa la norma euclídea sin tener en cuenta la matriz.

Conjugadas y no perpendiculares

p=(1,1), q=(4,−1) y A=diag(1,4).

  1. Euclídeo

    pᵀq=4−1=3.
  2. Con A

    Aq=(4,−4).
  3. Producto

    pᵀAq=4−4=0.
Resultado

La geometría relevante cambia con A.

Lo comprobamos

Conjugación no significa que los segmentos formen noventa grados en un dibujo euclídeo.

4. Las recurrencias calculan cada paso sin reiniciar toda la búsqueda

CG inicia p₀=r₀.

Mientras r_{k}≠0, usa \alpha_k=\frac{r_k^Tr_k}{p_k^TAp_k},\quad x_{k+1}=x_k+\alpha_kp_k,\quad r_{k+1}=r_k-\alpha_kAp_k. El denominador es positivo para una dirección no nula por SPD.

Dos pasos conjugados llegan a la solución del sistema
A=diag(1,4); b=(1,1)x₁=(2/5,2/5)x₂=(1,1/4)p₀=(1,1); p₁ ∥ (4,−1)p₀ᵀAp₁=0Conjugadas en la geometría de A.
Gradiente conjugado exacto para A diagonal uno cuatro, b uno uno e inicio cero cero. Ambos ejes de coordenadas usan la misma escala. El primer paso llega a dos quintos dos quintos; el segundo a uno un cuarto, la solución exacta. El primer segmento sigue p cero igual uno uno y el segundo una dirección proporcional a cuatro menos uno. Son conjugadas respecto a A porque (1,1) por A por (4,−1) da cero, aunque su producto euclídeo es tres. El dibujo muestra los iterados de este ejemplo de dos dimensiones en aritmética exacta, no una garantía de dos pasos para una matriz general.

Si el nuevo residuo es cero se termina.

Si no, \beta_k=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k},\qquad p_{k+1}=r_{k+1}+\beta_kp_k.

En dimensión n, CG termina en a lo sumo n pasos exactos. Con redondeo puede perder ortogonalidad; esa cifra no es un límite práctico garantizado.

Primer paso completo

A=diag(1,4), b=(1,1) y x₀=0.

  1. Dirección

    r₀=p₀=(1,1); p₀ᵀAp₀=5 y r₀ᵀr₀=2.
  2. Paso

    α₀=2/5 y x₁=(2/5,2/5).
  3. Nuevo residuo

    r₁=(1,1)−(2/5)(1,4)=(3/5,−3/5).
Resultado

El nuevo residuo es ortogonal al anterior: su producto vale cero.

Lo comprobamos

Aún no es la solución, porque r₁≠0.

Segundo paso y comprobación

Continuamos el mismo sistema.

  1. Dirección

    β₀=9/25 y p₁=(24/25,−6/25).
  2. Paso

    p₁ᵀAp₁=144/125 y r₁ᵀr₁=18/25; α₁=5/8.
  3. Llegada

    x₂=x₁+(5/8)p₁=(1,1/4); Ax₂=(1,1).
Resultado

El sistema está resuelto tras dos pasos exactos.

Lo comprobamos

Se detiene antes de intentar dividir por la norma de un residuo ya nulo.

5. La condición da una garantía y los autovalores dan más información

Para SPD, κ₂(A)=λ_max/λ_min.

Una cota clásica de CG es \|e_m\|_A\le2\left(\frac{\sqrt\kappa-1}{\sqrt\kappa+1}\right)^m\|e_0\|_A. Es una cota superior, no la trayectoria exacta del error.

Si A tiene pocos autovalores distintos, un polinomio de bajo grado puede anular todas las componentes del error. El agrupamiento espectral puede favorecer convergencia más rápida que la cota extrema. Para matrices no normales no se traslada este análisis como si todas fueran SPD.

Leer la cota

κ=9 y m=3.

  1. Raíz

    √κ=3.
  2. Cociente

    (3−1)/(3+1)=1/2.
  3. Cota

    ∥e₃∥A≤2·(1/2)³∥e₀∥A=(1/4)∥e₀∥A.
Resultado

La garantía reduce el error a como máximo una cuarta parte.

Lo comprobamos

El error real podría ser mucho menor o cero.

Muchos tamaños, solo dos autovalores

A es diagonal con entradas uno o cuatro.

  1. Polinomio

    p(λ)=(1−λ)(1−λ/4) cumple p(0)=1.
  2. Anular

    p(1)=p(4)=0.
  3. Consecuencia

    En aritmética exacta CG puede eliminar el error en a lo sumo dos pasos, aunque haya muchas componentes.

Resultado

La dimensión n no cuenta toda la historia.

Lo comprobamos

Se conservan las hipótesis SPD y los valores espectrales exactos.

6. GMRES minimiza el residuo en lugar del error de energía

Para sistemas generales invertibles, GMRES usa Arnoldi para construir una base ortonormal del espacio de Krylov. Resuelve un pequeño problema de mínimos cuadrados y elige el candidato con menor ∥b−Ax∥₂ en ese espacio afín.

Almacenar y ortogonalizar todos los vectores cuesta cada vez más. GMRES reiniciado limita la memoria, pero puede perder convergencia o estancarse. CG y GMRES pueden dar puntos distintos incluso sobre una matriz SPD porque minimizan cantidades distintas.

Comparar el primer paso

Mismo A=diag(1,4), b=(1,1) e inicio cero; x=α(1,1).

  1. Residual

    Su cuadrado es (1−α)²+(1−4α)².

  2. Minimizar

    Derivar da −10+34α=0, luego α=5/17.
  3. Comparar

    CG escogió 2/5: no son iguales.

Resultado

GMRES minimiza aquí el residual euclídeo.

Lo comprobamos

CG mantiene su garantía diferente de mínimo de energía.

Cuantificar la diferencia

Primeros residuos de ese ejemplo.

  1. GMRES

    r=(12/17,−3/17); norma cuadrada 9/17.
  2. CG

    r=(3/5,−3/5); norma cuadrada 18/25.
  3. Comparar

    9/17<18/25, aunque ambos usan la misma recta de candidatos.

Resultado

El mínimo residual no tiene por qué coincidir con el mínimo error en norma A.

Lo comprobamos

Una norma de error pequeña y un residual pequeño se relacionan mediante A.

7. PCG respeta la simetría mediante un cambio de variables

Escogemos M SPD y barato de resolver.

La matriz transformada \widehat A=M^{-1/2}AM^{-1/2} es SPD.

PCG implementa esa geometría usando z=M⁻¹r mediante resoluciones con M, sin necesitar formar la raíz o una inversa densa.

Un cambio de variables equilibra escalas distintas
A=M=diag(1,16)κ=16κ=1z=(x₁, 4x₂)M⁻¹ᐟ²AM⁻¹ᐟ²=ISe transforma la geometría del sistema.
Izquierda: una curva de nivel x uno cuadrado más dieciséis x dos cuadrado igual una constante, con ejes euclídeos igualmente escalados. Derecha: tras definir z uno igual x uno y z dos igual cuatro x dos, la misma expresión es z uno cuadrado más z dos cuadrado, con nivel circular. En este ejemplo A y M son ambas la diagonal uno dieciséis; M menos un medio A M menos un medio es la identidad y la condición pasa de dieciséis a uno. No se cambia arbitrariamente la escala del dibujo para afirmar que la condición mejora: se define el cambio de variables que transforma el sistema. En una matriz general aplicar un precondicionador perfecto puede ser tan costoso como resolver el sistema.

M⁻¹A es similar a la matriz transformada, pero puede no ser simétrica euclídeamente. Por eso aplicar CG ordinario a una matriz transformada arbitrariamente y llamarla simétrica no justifica el algoritmo. Las recurrencias de PCG usan productos rᵀz y la estructura de M.

Una transformación que equilibra todo

A=M=diag(1,16).

  1. Raíz inversa

    M⁻¹ᐟ²=diag(1,1/4).
  2. Matriz

    M⁻¹ᐟ²AM⁻¹ᐟ²=I.
  3. Geometría

    Con z=(x₁,4x₂), x₁²+16x₂² se convierte en z₁²+z₂².
Resultado

La condición del sistema transformado vale uno.

Lo comprobamos

En esta diagonal aplicar M⁻¹ es barato.

No suponer simetría del producto izquierdo

A=[[2,1],[1,2]] y M=diag(2,1), ambas SPD.

  1. Multiplicar

    M⁻¹A=[[1,1/2],[1,2]].
  2. Comparar

    Las entradas fuera de diagonal son distintas.

  3. Equivalente simétrico

    M⁻¹ᐟ²AM⁻¹ᐟ²=[[1,1/√2],[1/√2,2]], que sí es SPD.
Resultado

La justificación viene de la transformación simétrica.

Lo comprobamos

Un M positivo no garantiza por sí solo mejorar la condición: también importa su relación con A.

8. El objetivo práctico es reducir coste total, no solo iteraciones

M=A sería perfecto para la condición, pero resolver Mz=r puede ser el problema original. Un buen precondicionador combina construcción asequible, aplicación barata y mejora de convergencia. Jacobi, factorizaciones incompletas y multigrid explotan estructuras diferentes.

Con precondicionamiento izquierdo, GMRES minimiza un residual transformado. Con precondicionamiento derecho se conservan las ecuaciones originales en la minimización habitual. Hay que observar el residual que importa y contabilizar memoria, preparación y productos, además de iteraciones.

Comparar dos costes

Método A prepara en 100 unidades y usa diez iteraciones de cinco; B prepara en diez y usa treinta de cuatro.

  1. A

    100+10·5=150.
  2. B

    10+30·4=130.
  3. Elegir

    Para una sola resolución, B cuesta menos pese a más iteraciones.

Resultado

Menos iteraciones no equivale siempre a menos tiempo total.

Lo comprobamos

Con muchos términos derechos el coste de preparación puede amortizarse.

El residual transformado puede engañar si no se especifica

M=1000I y residual original r de norma uno.

  1. Izquierda

    M⁻¹r tiene norma 0,001.

  2. Original

    El fallo de Ax=b sigue teniendo norma uno.
  3. Comparación

    La cifra solo se interpreta sabiendo qué norma y transformación se ha usado.

Resultado

Se comprueba también el residual real del sistema original.

Lo comprobamos

No se deduce precisión mil veces mayor por una simple escala.

9. La parada comprueba el sistema real y las hipótesis del método

Durante la iteración se actualiza r por recurrencias para ahorrar trabajo. El redondeo puede separarlo de b−Ax; recalcular el residual real permite comprobar la parada. Una tolerancia combina a menudo una parte absoluta y otra relativa al dato o al residuo inicial.

Si el problema requiere una cota de error, se incorpora la condición o una estimación específica. También se revisan simetría, positividad y compatibilidad del precondicionador. Terminar en n pasos exactos, conservar ortogonalidad o una cota de convergencia no son afirmaciones independientes de esas condiciones.

Aplicar una regla de parada

Se pide ∥r∥≤max(10⁻⁸,10⁻⁶∥b∥), con ∥b∥=100.

  1. Parte relativa

    10⁻⁶·100=10⁻⁴.
  2. Umbral

    El máximo es 10⁻⁴.

  3. Residuo

    Con ∥r∥=5·10⁻⁵ se cumple esta regla.
Resultado

La regla garantiza el residual indicado, no automáticamente un error relativo igual.

Lo comprobamos

La condición conecta las dos afirmaciones.

Detectar una matriz fuera del marco de CG

A=diag(1,−1), dirección p=(1,1).

  1. Producto

    Ap=(1,−1).
  2. Denominador

    pᵀAp=0.
  3. Fallo

    La fórmula del paso puede dividir por cero; A no es positiva definida.

Resultado

Antes de ejecutar se verifican las hipótesis del método.

Lo comprobamos

La invertibilidad de esta matriz no la convierte en SPD.

Ahora lo resolvemos de principio a fin.

EJEMPLO RESUELTO 1

Condición diagonal

Calcula κ₂ de diag(1,100).

  1. Leer la matriz

    A=diag(1,100).
  2. Valores singulares

    En esta diagonal positiva son uno y cien.

  3. Máximo

    La norma espectral de A vale cien.

  4. Inversa

    A⁻¹=diag(1,1/100), de norma espectral uno.
  5. Condición

    κ₂(A)=100·1=100.
  6. Comprobar

    También es el cociente entre autovalor máximo y mínimo porque A es SPD.

Resultado

La condición espectral es 100.

Lo comprobamos

La disparidad de escalas puede frenar métodos iterativos.

EJEMPLO RESUELTO 2

Precondición ideal

Para esa diagonal, usa M=diag(1,100). Calcula M⁻¹A.

  1. Escoger el precondicionador

    M=diag(1,100)=A.
  2. Aplicar inversa

    Resolver Mz=r divide la segunda componente por cien.

  3. Multiplicar

    M⁻¹A=diag(1,1)=I.
  4. Condición

    La identidad tiene condición uno.

  5. Geometría

    La transformación simétrica equivalente también es la identidad.

  6. Coste

    Aquí resolver con M es barato por ser diagonal; para una matriz general elegir M=A puede ser costoso.

Resultado

M⁻¹A=I, de condición 1.

Lo comprobamos

En este caso aplicar M⁻¹ es barato porque es diagonal.

Los símbolos que hemos usado

r₀=b−Ax₀
Residuo de la aproximación inicial.
K_m(A,r₀)
Span de r₀, Ar₀, …, A^(m−1) r₀.
SPD
Matriz real simétrica y definida positiva.
‖e‖_A
Norma de energía √(eᵀAe), cuando A es SPD.
κ₂(A)
Condición en norma euclídea; para SPD es λ_max/λ_min.
M
Precondicionador cuya acción inversa se puede aplicar a coste razonable.
GMRES
Método que minimiza la norma euclídea del residuo sobre un espacio de búsqueda de Krylov.

Ahora te toca a ti.

Tienes cinco ejercicios. Inténtalos; las pistas te ayudan si te atascas.

Practicar este tema →

Para ampliar