Métodos potenciales de prospección, FCAG, 2023.
Buscamos resolver el siguiente sistema lineal de dos ecuaciones con dos incógnitas:
\begin{equation} A\mathbf{m} = \mathbf{y}, \end{equation}con $A=\begin{bmatrix} 400 & 1 \\ 1 &4\end{bmatrix}$ y $\mathbf{y}=(5,2)^T$. La solución de este sistema es $\mathbf{m}^*=(0.0113,0.4970)^T$.
La función de costo o función objetivo es:
\begin{equation} q = q(\mathbf{m}) := \dfrac{1}{2}|| \mathbf{y} - A\mathbf{m} ||^2. \end{equation}El gradiente de la función de costo respecto de $\mathbf{m}$ es:
\begin{equation} \nabla_{\mathbf{m}}q = A\mathbf{m} - \mathbf{y}. \end{equation}En esta actividad práctica, realizamos el preacondicionamiento del sistema de ecuaciones:
\begin{equation} P\,A\mathbf{m} = P\mathbf{y}, \end{equation}Donde P es una matriz diagonal $P_{j,j}=1/A_{j,j}$.
El preacodicionamiento genera que las lineas de contorno de la función de costo sean circulares y el método de gradient descent no tarde en converger hacia la solución. Cotejamos visualmente esto último.
array([[0.0025, 0. ],
[0. , 0.25 ]])
Programamos el método de gradient descent y las subrutinas necesarias:
Procedemos a realizar la inversión, partiendo de un vector inicial $\mathbf{m}_0=(0,0)^T$ y sin realizar preacondicionamiento:
Gradient descent (sin preacondicionamiento) a = 0.011257035653554502 b = 0.49718573860328286 número de iteraciones: 4772 tamaño de paso: 0.001
Visualizamos el resultado del método gradient descent en función del número de iteraciones $k$:
Necesitamos de varias iteraciones $k$ para obtener el resultado deseado.
Realizamos ahora el preacondicionamiento del sistema:
P = np.diag([1/400.,1/4.])
PA = P @ A
Py = P @ y
Procedemos a realizar la inversión, partiendo de un vector inicial $\mathbf{m}_0=(0,0)^T$ y tomando un tamaño de paso mucho mayor que en el caso anterior.
Gradient descent (con preacondicionamiento) a = 0.011257035811693463 b = 0.49718573715068803 número de iteraciones: 27 tamaño de paso: 0.5
Visualizamos en función del número de iteraciones $k$:
En la notebook del método de la capa equivalente utilizaremos una normalización similar para lograr una optimización aceptable. En ese caso la normalización es respecto a la desviación estándard de cada columna de $A$:
scale = 1/np.std(A,axis=1)
P = np.diag(scale)
PA = P @ A
Py = P @ y
a = 0.011257035651043495 b = 0.49718574158878504 a = 0.011257035651043497 b = 0.49718574158878504
Nota. Para funciones de costo cuadráticas, puede calcularse un tamaño de paso exacto por medio de la expresión:
\begin{equation} \epsilon = \dfrac{\mathbf{g}^T\mathbf{g}}{\mathbf{g}^TA\mathbf{g}}, \end{equation}donde el vector $\mathbf{g}=A\mathbf{m}-\mathbf{y}$ es el gradiente de la función de costo. De haber utilizado este resultado para proponer el tamaño de paso, hubieramos obtenido resultados satisfactorios en una menor cantidad de iteraciones.
Eso es todo por hoy.