Regularización de Tikhonov

Question

Estoy revisando algo de código en C que agrega regularización de Tikhonov a una matriz simétrica definida positiva.

Así:

power = (matrix[1][1] + matrix[2][2]);

factor = 0.000136 //Encontrado mediante prueba y error

for(i = 1; i < matrixSize; i++ ) {

matrix[i][i] = matrix[i][i]+( factor * power );

}

Estoy atascado en entender qué significa ese power, ¿por qué se han utilizado las segunda y tercera partes de la diagonal? ¿Cómo está ayudando esto a hacer que los números sean más estables?

Gracias por cualquier ayuda.

Answer 1

Entiendo que mi respuesta puede ser demasiado tardía en cierto sentido, pero me gustaría publicar algunas explicaciones sobre el enfoque de regularización de Tikhonov de todos modos, porque bueno, en primer lugar, tuve una experiencia práctica en la regularización aplicada y en el enfoque de Tikhonov también resolviendo problemas científicos inversos reales (expresados principalmente en forma de ecuaciones integrales) que eran fuertemente mal planteados, y en segundo lugar, cuando comencé me encontré con una notable falta de fuentes que pudieran describir claramente una implementación práctica de la regularización y, lo que fue más importante, de la estimación de la solución.

Hablando en términos simples, el enfoque de Tikhonov no es más que una combinación regularizada del método de mínimos cuadrados de Gauss y el método de matriz seudo-inversa de Moore-Penrose para encontrar una solución estable que existe y está definida de manera única en el espacio de soluciones. Pero a diferencia de los mínimos cuadrados y el método de matriz seudo-inversa, abordamos una solución desconocida real con la posibilidad de una estimación del error de la solución regularizada si se cumplen algunos requisitos relacionados con los datos iniciales y los espacios de soluciones. ¡En caso de una estimación del error punto a punto, es posible definir una región (límites superior e inferior) donde se encuentra la solución exacta!

Para tener la posibilidad de abordar el problema publicado de manera más general, usemos una notación matricial y algo de análisis funcional también (solo por diversión):

Problema inverso inicial:

$A\tilde{y}=\tilde{b}, \tilde{y} \in L_2, \tilde{b} \in L_2\\ \|\bar{b}-\tilde{b}\|_{L_2} \leq \delta,$

donde $A$ es una matriz de coeficientes del SLE que se obtiene después de la discretización numérica del problema inicial, $\tilde{b}$ representa unos datos iniciales, $\delta$ es un error superior de los datos iniciales (si los datos iniciales no tuvieran errores, entonces el problema inicial se volvería correctamente planteado y podría resolverse utilizando métodos comunes) y $\tilde{y}$ es una solución regularizada aproximativa del problema inverso. El error del operador $\|\tilde{A} - \bar{A}\| \leq \xi$ no se tiene en cuenta ya que se asume que es insignificante si la cuadrícula de discretización es lo suficientemente fina.

Aquí se ve que considero deliberadamente que los datos iniciales y la solución regularizada pertenecen al espacio de soluciones de Hilbert $L_2$. Hay algunas razones significativas para ello. En pocas palabras, esta restricción garantiza (junto con algunas afirmaciones adicionales que no proporciono aquí, ya que están relacionadas con la comprensión más profunda de la regularización en general) una existencia de la solución regularizada que puede encontrarse resolviendo el siguiente problema de minimización:

$M_{\delta}^{\alpha}[\tilde{y}]=\|A\tilde{y}-\tilde{b}\|_{L_2}^2 + \alpha\|\tilde{y}\|_{L_2}^2, \\ y_{\alpha}=\text{arg}\,\min M_{\delta}^{\alpha}[\tilde{y}],$

donde $M_{\delta}^{\alpha}[\tilde{y}]$ es un funcional de suavizado, $\alpha$ es un parámetro de regularización y $y_{\alpha}$ es una solución regularizada que proporciona un mínimo global del funcional $M$.

En el caso dado, la minimización resulta en la ecuación de segundo orden de Euler-Tikhonov:

$\alpha y_{\alpha}+A^TAy_{\alpha}=A^T\tilde{b}, \\ y_{\alpha}=(\alpha I + A^TA)^{-1}A^T\tilde{b},$

donde $I$ es una matriz identidad.

Así, como se puede ver, el enfoque numérico simple no es tan malo como se describe en aburridos trabajos matemáticos :) Pero hasta ahora no he presentado otra cosa más que es responsable del éxito de toda la historia. Sí, eso es correcto, ese es el parámetro de regularización $\alpha$. Todo depende de él. Incluso si se demuestra la convergencia del enfoque de regularización, todavía hay una fuerte demanda en el métod...

• Defina un vector de los parámetros de regularización que se utilizarán para calcular un espacio de soluciones regularizadas. Por ejemplo, en este orden:

$\alpha_n = \alpha_0\cdot q^n, n = \overline{1, N},$

donde $N$ denota el número total de las soluciones regularizadas que componen un espacio de soluciones; $\alpha_0$ es un valor inicial, $\alpha_0 > 0$; $q$ es un factor que rige la tasa de convergencia hacia la solución exacta (cuanto menor sea el factor q, más rápida será la tasa), $0 < q < 1$;

• Calcular la ecuación de Euler-Tikhonov para cada parámetro $\alpha_n$;
• Elegir el parámetro más apropiado que corresponda a la solución que está más cerca de la desconocida.

Uno de los métodos más ampliamente utilizados es el principio de la discrepancia. Usando este método, la solución se encuentra como una raíz de la siguiente ecuación:

$\|Ay_{\alpha}-\tilde{b}\|_{L_2}^2 = \delta^2$

Es bastante directo. Pero el diablo se esconde en los detalles. Aquí cabe mencionar que el enfoque de regularización opera con la ecuación de Euler-Tikhonov pero no con la inicial. El análisis detallado del enfoque de regularización dado con atención a la característica señalada muestra que el uso del principio de la discrepancia siempre conducirá a soluciones sobre-regularizadas (el parámetro de regularización $\alpha$ será innecesariamente grande). Tal peculiaridad no es para nada obvia y puede encontrarse muy frecuentemente que muchos investigadores no prestan toda la atención necesaria a estos hechos y tratan las soluciones regularizadas calculadas como si estuvieran lo suficientemente cerca de las soluciones exactas cuando no lo están. Por lo tanto, es probable que se necesite información adicional sobre la solución desconocida (por ejemplo, la suavidad de la solución).

Una palabra más sobre el procedimiento de elección. El requisito mínimo es $\delta$. Eso significa que no hay forma de elegir el parámetro "correcto" $\alpha$ sin el error de los datos iniciales. Sé que, de vez en cuando, surgen intentos de encontrar un método que pueda encontrar una solución regularizada estable y confiable sin error en los datos iniciales. Método $L$-curva, por ejemplo. Pero de todos modos, está teórica y prácticamente demostrado que tales métodos no pertenecen a la familia de algoritmos de regularización y no pueden llevar consecuentemente a la solución exacta.

La estimación del error de la solución regularizada es otra gran historia que no tiene lugar para ser proporcionada aquí...

P.D.: El enfoque de regularización de Tikhonov no es la única forma de resolver problemas mal planteados numéricamente, pero desde mi experiencia, este es muy estable, flexible y relativamente consume no demasiado tiempo para comenzar con el método que debería estar entre los primeros candidatos a probar.