[[2 1] [1 2]] Energia en la soluciĆ³n: [[-1.]] xnew [[0.35714286 0.71428571]]
Energia en xnew: [[-0.89285714]] xnew [[-5.55111512e-17 1.00000000e+00]]
Energia en xnew: [[-1.]]
MĆ©todos potenciales de prospecciĆ³n, FCAGLP, 2020.
MPP-2020.
$\newcommand{\vec}[1]{\mathbf{#1}}$ TambiĆ©n conocido como mĆ©todo de Stenes y Stiefel. "Es un algoritmo fantĆ”sico" (Strang) y "extraordinariamente simple" (Trefehen y Bau). Tiene una gran aplicaciĆ³n en la inversiĆ³n de problemas geofĆsicos, en particular los mĆ©todos potenciales y la sĆsmica de prospecciĆ³n.
Asumimos que $A$ $(m\times m)$ es una matriz a coeficientes reales, simĆ©trica y ademĆ”s definida positiva. Esto implica que los autovalores de $A$ son todos positivos, o equivalentemente $\vec{x}^TA\vec{x}>0$ para cada $\vec{x}\in \mathbb{R}^m$ distinto del vector $\vec{0}$. Con estas suposiciones, la funciĆ³n $||\cdot||_A$ definida por $||\vec{x}||_A=\vec{x}^TA\vec{x}$ define una norma en $\mathbb{R}^m$. Podemos llamar a esta norma la "A-norma".
Lo que sigue es minimizar la "A-norma" del vector error, entre la soluciĆ³n exacta y la aproximaciĆ³n.
El vector con A-norma de interĆ©s es $\vec{e}_n = \vec{x}^*-\vec{x}_n$, el error en el paso $n$ del algoritmo iterativo de CG respecto de la soluciĆ³n $\vec{x}^*$ (desconocida). El proceso iterativo de CG puede resumirse entonces como: "es un sistema que genera una Ćŗnica secuencia de vectores $\vec{x}_n$ con la propiedad que en cada paso $n$, $||\vec{e}_n||_A$ es minimizado".
Nota. Los iterados de CG cumplen con $\vec{x}_n \in \overline{\{\vec{b},A\vec{b},A^2\vec{b},...,A^{n-1}\vec{b}\}}$. El espacio generado por los vectores $\{\vec{b},A\vec{b},A^2\vec{b},...,A^{n-1}\vec{b}\}$ recibe el nombre de espacio de Krylov, $\mathcal{K}_n$. Es decir, la aproximaciĆ³n a la soluciĆ³n generada en el paso $n$ de CG pertenece al espacio generado por las $n-1$ potencias de $A$ aplicadas a $\vec{b}$.
$\newcommand{\vec}[1]{\mathbf{#1}}$
$\vec{x}_0 = \vec{0}$ (soluciĆ³n inicial), $\vec{r}_0 = \vec{b}$ (residuo inicial), $\vec{d}_0 = \vec{r}_0$ (direcciĆ³n de bĆŗsqueda inicial).
Para $n=1,2,3,...$
\begin{align} \alpha_n &= \big(\vec{r}_{n-1}^T \vec{r}_{n-1}\big) / \big(\vec{d}_{n-1}^TA\vec{d}_{n-1}\big); \\ \vec{x}_n &= \vec{x}_{n-1} + \alpha_n \vec{d}_{n-1}; \\ \vec{r}_n &= \vec{r}_{n-1} - \alpha_n A\vec{d}_{n-1};\\ \beta_n &= \big(\vec{r}_{n}^T \vec{r}_{n}\big)/\big(\vec{r}_{n-1}^T \vec{r}_{n-1}\big);\\ \vec{d}_n &= \vec{r}_n + \beta_n \vec{d}_{n-1}. \end{align}Como puede observarse, el algoritmo involucra sĆ³lo un producto matriz por vector: $A\vec{d}_{n-1}$.
$\alpha_n$ es el tamaƱo de paso, $\vec{x}_n$ la soluciĆ³n aproximada, $\vec{r}_n$ el residuo, $\beta_n$ es la mejora en el paso dado y $\vec{d}_n$ la nueva direcciĆ³n de bĆŗsqueda.
En el mĆ©todo iterativo encontramos la fĆ³rmula
$$\vec{x}_n = \vec{x}_{n-1} + \alpha_n \vec{d}_{n-1},$$que es una expresiĆ³n familiar en los mĆ©todos de optimizaciĆ³n. Esta fĆ³rmula nos dice que la soluciĆ³n actual $\vec{x}_{n-1}$ es actualizada a una nueva aproximaciĆ³n $\vec{x}_n$ moviĆ©ndonos una cantidad $\alpha_n$ en la direcciĆ³n $\vec{d}_{n-1}$.
Puede probarse que mientras el mƩtodo se aplique a una matriz simƩtrica definida positiva y las iteracciones no hayan convergido (i.e., $\vec{r}_n\neq \vec{0}$), el algoritmo procede sin divisiones por cero; los residuos son ortogonales,
$$\vec{r}_n^T\vec{r}_j=0, j<n;$$y las direcciones de bĆŗsqueda,
$$\vec{d}_n^TA\vec{d}_j=0,$$son "A-ortogonales" o "A-conjugadas", de donde el mƩtodo recibe su nombre.
Por Ćŗltimo, puede probarse ademĆ”s que mientras el mĆ©todo no haya convergido (i.e., $\vec{r}_n\neq \vec{0}$),
$$||\vec{e}_n||_A \leq ||\vec{e}_{n-1}||_A$$y $\vec{e}_n=\vec{0}$ se logra para algĆŗn $n\leq m$. Todas las propiedades mencionadas son vĆ”lidas suponiendo que no existen errores aritmĆ©ticos en los cĆ”lculos.
CG puede interpretarse como un proceso iterativo para minimizar la funciĆ³n cuadrĆ”tica $E(\vec{x})=\frac{1}{2} \vec{x}^TA\vec{x}-\vec{x}^T\vec{b}$ para $\vec{x} \in \mathbb{R}^m$. La funciĆ³n $E(\vec{x})$ podemos asimilarla a una "energĆa". Vemos que $E(\vec{x})$ pasa por 0 cuando $\vec{x}=\vec{0}$ y tiene su Ćŗnico mĆnimo en $E(\vec{x^*})=-\frac{1}{2}\vec{x^*}^T\vec{b}$.
Esta funciĆ³n es similar a $||\vec{e}_n||_A$, salvo por una constante.
En cada paso, la iteraciĆ³n $\vec{x}_n = \vec{x}_{n-1} + \alpha_n \vec{d}_{n-1}$ minimiza a $E(\vec{x})$ para todos los $\vec{x}$ en el rayo generado por $\vec{x}_{n-1}+\overline{\vec{d}_n}$. AquĆ, $\overline{\vec{d}_n}$ denota el espacio generado por $\vec{d}_n$. Vemos a continuaciĆ³n un ejemplo concreto para fijar el concepto.
Ejemplo. El primer subespacio ($\mathcal{K}_1$) es el rayo generado en la direcciĆ³n $\vec{d_0} = \vec{r_0} = \vec{b}$. La energĆa para los vectores $\vec{x} = \alpha\vec{b}$ es $E(\alpha\vec{b})=\frac{1}{2} \alpha\vec{b}^TA\alpha\vec{b}-\alpha\vec{b}^T\vec{b}$. Minimizando $E(\alpha\vec{b})$ en funciĆ³n de $\alpha$ obtenemos el nĆŗmero $\alpha_1=\frac{\vec{b}^T\vec{b}}{\vec{b}^TA\vec{b}}$. Como vemos, este valor de $\alpha_1$ es el tamaƱo de paso dado en el primer ciclo del algoritmo CG.
Cada ciclo de CG elige $\alpha_k$ para minimizar a $E(\vec{x})$ en la direcciĆ³n de bĆŗsqueda $\vec{x}=\vec{x}_n+\alpha\vec{d}_{n-1}$. El Ćŗltimo ciclo obtiene el minimizador $\vec{x^*}=A^{-1}\vec{b}.$
A modo de comparaciĆ³n, el mĆ©todo de gradient descent es un caso particular del algoritmo para minimizar a $E(\mathbf{x})$, donde se elige $\vec{d}_n=\vec{r}_n=-\boldsymbol\nabla{E}(\mathbf{x})$, en vez de $\vec{d}_n = \vec{r}_n + \beta_n \vec{d}_{n-1}$.
CG: SegĆŗn la tolerancia 1e-15, convergiĆ³ en la iteraciĆ³n 2 [[-5.55111512e-17] [ 1.00000000e+00]]
con minimizador $\mathbf{x}^*=(0,1)^T$.
La matriz $A= \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ es simƩtrica y semidefinida positiva.
Define la siguiente funciĆ³n cuadrĆ”tica a minimizar (convex minimization):
$f(x)= \frac{1}{2}\mathbf{x}^TA\mathbf{x} - \mathbf{x}^T\mathbf{b}$
El mĆ©todo de gradientes conjugados (CG) debe arribar a la soluciĆ³n a lo sumo en dos pasos (nĆŗmero de columnas de $A$).
[[2 1] [1 2]] Energia en la soluciĆ³n: [[-1.]] xnew [[0.35714286 0.71428571]]
Energia en xnew: [[-0.89285714]] xnew [[-5.55111512e-17 1.00000000e+00]]
Energia en xnew: [[-1.]]
En este ejemplo podemos hacer las cuentas en papel y analizar en detalle al algoritmo. Consideramos el caso $A\vec{x}=\vec{b}$ dado por
\begin{equation} \begin{bmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{bmatrix} \begin{bmatrix} \phantom{+}3 \\ -1 \\ -1 \end{bmatrix}= \begin{bmatrix} 4 \\ 0 \\ 0 \end{bmatrix}, \end{equation}claramente el minimizador es $\mathbf{x}^*=(3,-1,-1)^T$.
Partiendo de $\vec{x_0}=\vec{0}, \vec{r_0}=\vec{d_0}=\vec{b}$, en la primer iteraciĆ³n resultan:
\begin{align} \alpha_1 &= \dfrac{1}{2}; \\ \vec{x}_1 &= \vec{0} + \alpha_1 \vec{b} = (2,0,0)^T; \\ \vec{r}_1 &= \vec{b} - A\vec{x}_1 = (4,0,0)^T-(4,2,2)^T=(0,-2,-2)^T ;\\ \beta_1 &= \frac{8}{16};\\ \vec{d}_1 &= \vec{r}_1 + \beta_1\vec{d}_0 = (0,-2,-2)^T + \frac{1}{2}(4,0,0)^T = (2,-2,-2)^T. \end{align}En la segunda iteraciĆ³n:
\begin{align} \alpha_2 &= \dfrac{1}{2}; \\ \vec{x}_2 &= \vec{x}_1 + \alpha_2\vec{d}_1 = (2,0,0)^T + \frac{1}{2}(2,-2,-2)^T = (3,-1,-1)^T\equiv\vec{x^*}; \\ \vec{r}_2 &= \vec{r}_1 - \alpha_2 A\vec{d}_1 = (0,-2,-2)^T-\frac{1}{2}(0,-4,-4)^T=\vec{0}.\\ \end{align}Normalmente, CG tomarĆa 3 iteraciones en llegar a la soluciĆ³n. En este caso alcanzĆ³ con 2 ciclos, donde el residuo es $\vec{r}_2=\vec{0}$ y el algoritmo entonces termina. La razĆ³n de ello es que $A$ es de rango 2, lo que implica que la soluciĆ³n $A^{-1}\vec{b}$ es una combinaciĆ³n de los vectores$\{\vec{b},A\vec{b}\}$. En efecto, un simple cĆ”lculo prueba que
$$\vec{x}^* = (\frac{5}{4})\vec{b}+(-\frac{1}{4})A\vec{b} = \frac{5}{4} (4,0,0)^T -\frac{1}{4}(8,4,4)^T,$$verificando $\vec{x}^* \in \overline{\{\vec{b},A\vec{b}\}}.$
Por Ćŗltimo, observamos en la primer iteraciĆ³n que $\vec{r}_1$ dado por $\vec{b}-A\vec{x}_1$ es el gradiente de $E(\vec{x})$ cambiado de signo. La direcciĆ³n que tomarĆa gradient descent (GD) es justamente esta. La dificultad de GD es que $\vec{r}_1$ puede ser muy similar a $\vec{d}_0$. Esto conduce a poco progreso en la reducciĆ³n de $E(\vec{x})$. Para ello, CG agrega la componente $\beta_1\vec{d}_0$, forma tal que la nueva direcciĆ³n $\vec{d}_1=\vec{r}_1+\beta_1\vec{d}_0$ es A-ortogonal a $\vec{d}_0$ (lo que garantiza un progreso sostenido hacia el mĆnimo). Luego nos movemos en esta direcciĆ³n conjugada $\vec{d}_1$ hacia $\vec{x}_2=\vec{x}_1+\alpha_1\vec{d}_1$.
x0 = [[0 0 0]] ---------------- iteraciĆ³n 1 ---------------- Ī± = [[0.5]] Ī² = [[0.5]] xnew = [[2. 0. 0.]] residuo = [[ 0. -2. -2.]] ---------------- iteraciĆ³n 2 ---------------- Ī± = [[0.5]] Ī² = [[0.]] xnew = [[ 3. -1. -1.]] residuo = [[0. 0. 0.]] ++++++++++++ Cerrar loop.
En este ejemplo mostramos solamente como el residuo tiene a cero a medida que progresan las iteraciones.
Consideramos el caso $A\vec{x}=\vec{b}$ dado por
\begin{equation} \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 1 \\ 1 & 1 & 3 & 1 \\ 1 & 1 & 1 & 4 \end{bmatrix} \begin{bmatrix} 1 \\ 1 \\ 1 \\ 1 \end{bmatrix}= \begin{bmatrix} 4 \\ 5 \\ 6 \\ 7 \end{bmatrix}, \end{equation}claramente el minimizador es $\mathbf{x}^*=(1,1,1,1)^T$.
A [[1 1 1 1] [1 2 1 1] [1 1 3 1] [1 1 1 4]] b= [[4] [5] [6] [7]] x0 [[0 0 0 0]] ---------------- iteraciĆ³n 1 ---------------- Ī± = [[0.17307692]] Ī² = [[0.00280008]] xnew = [[0.69230769 0.86538462 1.03846154 1.21153846]] residuo = [[ 0.19230769 0.32692308 0.11538462 -0.44230769]] ---------------- iteraciĆ³n 2 ---------------- Ī± = [[0.46935747]] Ī² = [[0.07118769]] xnew = [[0.78782571 1.02539961 1.10050361 1.01313773]] residuo = [[ 0.07313335 0.04773374 -0.12787388 0.03372017]] ---------------- iteraciĆ³n 3 ---------------- Ī± = [[0.71154748]] Ī² = [[0.12210046]] xnew = [[0.85017193 1.07663337 1.01621091 1.01571967]] residuo = [[ 0.04126412 -0.03536925 0.00884231 -0.00589487]] ---------------- iteraciĆ³n 4 ---------------- Ī± = [[2.88338055]] Ī² = [[1.52917007e-27]] xnew = [[1. 1. 1. 1.]] residuo = [[-7.49400542e-16 -9.50628465e-16 -1.18828558e-15 -1.34614542e-15]]
En este Ćŗltimo ejemplo, tomamos una matriz definida positiva $500\times 500$. Paulatinamente vamos perturbando la matriz para que vaya perdiendo la propiedad de ser definida positiva. Esto lo hacemos con el parĆ”metro $\tau$. Podemos observar que CG deja de converger cuando la matriz deja de ser definida positiva. Para $\tau=0.2$ la matriz no es definida positiva, entonces el mĆ©todo no converge. Vemos las primeras 20 iteraciones.
Si continuamos las iteraciones en el caso donde $A$ ya no es definida positiva veremos, experimentalmente, que se logra convergencia, pero no tenemos teoremas que garanticen el comportamiento en esta situaciĆ³n.
Eso es todo por hoy.