¿Por qué el gradiente conjugado funciona con este preacondicionador no simétrico?

En este hilo anterior, se sugirió la siguiente forma multiplicativa de combinar preacondicionadores simétricos y para el sistema simétrico : \ begin {align} P_ \ text {combo} ^ {- 1}: = & P_1 ^ {- 1} + P_2 ^ {- 1} (I - A P_1 ^ {- 1}) \\ = & P_1 ^ {- 1} + P_2 ^ {- 1} - P_2 ^ {- 1} A P_1 ^ {- 1}. \ end {align}P 2 A x = b P - 1 combo : = P - 1 1 + P - 1 2 ( I - A P - 1 1 ) = P - 1 1 + P - 1 2 - P - 1 2 A P - 1 1$P_1$$P_2$$Ax=b$

$$\begin{align} P_\text{combo}^{-1} :=& P_1^{-1} + P_2^{-1}(I - A P_1^{-1})\\ =&P_1^{-1} + P_2^{-1} - P_2^{-1} A P_1^{-1}. \end{align}$$

Este preacondicionador combinado no es simétrico. Sin embargo, he intentado usarlo en gradiente conjugado de todos modos en varios contextos diferentes, y el método siempre parece converger bien. ¿Por qué es esto?

---

Ejemplo 1: matrices aleatorias.

% Testing multiplicative combination of preconditioners
n = 500;
[U,S,~] = svd(randn(n,n)); A = U*S*U';
[W1,S1,~] = svd(randn(n,n)); noise1 = W1*S1*W1';
[W2,S2,~] = svd(randn(n,n)); noise2 = W2*S2*W2';
P1 = A + 0.5 * noise1;
P2 = A + 0.5 * noise2;

solve_P = @(x) P1\x + P2\x - P2\(A*(P1\x));

b = randn(n,1);
x_true = A\b;
pcg(A,b,1e-6,n);
pcg(A,b,1e-6,n,P1);
x = pcg(A,b,1e-6,n,solve_P);
norm(x_true - x)/norm(x_true)

Aquí está la salida que obtengo:

pcg converged at iteration 127 to a solution with relative residual 9.9e-07.
pcg converged at iteration 62 to a solution with relative residual 6.8e-07.
pcg converged at iteration 51 to a solution with relative residual 8.1e-07.
relative error= 4.23e-07

---

Ejemplo 2: Combinando multigrid con cholesky incompleto para resolver la ecuación de Poisson.

n = 50; N = n^2;
A = gallery('poisson',n); %Laplacian on n-by-n grid, zero dirichlet BC

L = ichol(A);
solve_P1 = @(x) L'\(L\x);
% Combinatorial multigrid: http://www.cs.cmu.edu/~jkoutis/cmg.html
solve_P2 = cmg_sdd(A); 

solve_P = @(x) solve_P1(x) + solve_P2(x) - solve_P1(A*solve_P2(x));

b = randn(N,1);
x_true = A\b;
pcg(A,b,1e-6,N);
pcg(A,b,1e-6,N,solve_P1);
pcg(A,b,1e-6,N,solve_P2);
x = pcg(A,b,1e-6,N,solve_P);
disp(['relative error= ', num2str(norm(x_true - x)/norm(x_true),3)])

Para esto obtengo los resultados:

pcg converged at iteration 131 to a solution with relative residual 8.4e-07.
pcg converged at iteration 44 to a solution with relative residual 6e-07.
pcg converged at iteration 19 to a solution with relative residual 7e-07.
pcg converged at iteration 12 to a solution with relative residual 4.7e-07.
relative error= 5.2e-07

---

Notas:

• También he visto el mismo comportamiento cualitativo en matrices que surge en situaciones más complicadas / realistas.
• Dada una solución incorrecta a , el error el residual están relacionados por la ecuación de error . Uno puede ver este preacondicionador combinado como aproximadamente resolver la ecuación original usando lugar de , luego aproximadamente resolver la ecuación de error con lugar de , y luego agregar el error aproximado para corregir la solución aproximada original.$Ax=b$r A e = r P 1 A P 2$e$$r$$Ae=r$$P_1$$A$$P_2$$A$

En resumen, la ortogonalización de los vectores de Krylov ocurre con respecto al operador, pero no con respecto al preacondicionador.

Muy bien, así que decir que queremos resolver con acondicionador previo . La iteración CG precondicionada es básicamente:$Ax=b$$B$

$$\begin{align*} \hat{v}_1=\tilde{v}_1 =& Bb\\ v_1 =& \tilde{v}_1 / c_1\\ \\ \hat{v}_i =& BAv_{i-1}\\ \tilde{v}_i=&\hat{v}_i-\sum\limits_{j=1}^{i-1} \langle Av_k,\hat{v}_i\rangle v_k\\ v_i =& \tilde{v}_i/c_i\\ c_i =& \sqrt{\langle \tilde{v}_i,\tilde{v}_i\rangle} \end{align*}$$ En otras palabras

• $\hat{v}_i$ - Nuevo vector de Krylov
• $\tilde{v}_i$ - Vector de Krylov ortogonalizado
• $v_i$ - Vector de Krylov normalizado

Ahora, lo interesante de esto es que obtenemos una serie de propiedades. Los dos que realmente te importan son

1. $\textrm{span}\{v_i\}_{i=1}^m = \textrm{span}\{(BA)^{i-1}(Bb)\}_{i=1}^m$
2. $V^TAV = I$

La primera propiedad se deduce de cómo se genera . Estamos restando cosas y normalizando, pero básicamente estamos obteniendo un montón de donde es una escala de . La segunda propiedad se sigue de . Aquí, no usamos en la ortogonalización. La ortogonalización depende únicamente de que sea ​​simétrico positivo definido.$\hat{v}_i$$BAv_{i-1}$$v_1$$Bb$$\tilde{v}_i$$B$$A$

Ahora, digamos que queremos resolver el sistema . Si elegimos en nuestro espacio Krylov, resolvemos en Sin embargo, esto probablemente esté sobredeterminado. Por lo tanto, resolvemos las ecuaciones normales, de modo que Desde arriba, sabemos que , por lo que realmente tenemos que o que Tenga en cuenta, esto es sin tener en cuenta el preacondicionador . De hecho, podríamos tomar vectores aleatorios y -ortogonalizarlos uno contra el otro y lo mismo sería válido. Realmente, lo único que afecta es el espacio de Krylov. Mientras$Ax=b$$x$$\beta$

$$AV\beta = b$$ $$V^TAV\beta = V^Tb$$ $V^TAV=I$ $$\beta = V^Tb$$ $$x = VV^Tb.$$ $B$$A$$B$$B$ ayuda a agrupar los valores propios de , el algoritmo funciona bien. Sin embargo, si no es simétrico, es posible que tengamos que mirar algo así como los pseudopectra para descubrirlo. No estoy del todo seguro. Sobre todo, agrupe los espectros.$A$$B$

¿Perdemos algo con una no asimétrica ? Probablemente. Normalmente, en CG, tenemos propiedades como para todo . Estoy bastante seguro de que ya no se sostiene. Además, normalmente tenemos para . Eso puede o no puede sostenerse. Realmente, cualquier resultado que tenga sería sospechoso.⟨ B r i , r j ⟩ = 0 i ≠ j ⟨ B A v i , A v j ⟩ = 0 | i - j | ≥ 2 $B$$\langle Br_i,r_j\rangle = 0$$i\neq j$$\langle BAv_i,Av_j\rangle=0$$|i-j|\geq 2$$B$

De todos modos, estoy bastante seguro de que es posible codificar CG de una manera que rompería todo con una no asimétrica . Sin embargo, si lo codificamos usando las ecuaciones anteriores, no hay nada que realmente le prohíba trabajar con una no simétrica . Lo hago todo el tiempo y parece funcionar bien.$B$$B$

Lea http://www.sciencedirect.com/science/article/pii/S1877050915010492 "Precondicionamiento no simétrico para gradientes conjugados y métodos de descenso más pronunciados"