Estrategias para depurar problemas de estabilidad numérica?

Question

Estoy tratando de escribir una implementación del algoritmo de factorización de densidad espectral de Wilson [1] para Python. El algoritmo iterativamente factoriza una función de matriz [QxQ] en su raíz cuadrada (es una especie de extensión del buscador de raíz cuadrada Newton-Raphson para matrices de densidad espectral).

El problema es que mi implementación solo converge para matrices de tamaño 45x45 y menor. Entonces, después de 20 iteraciones, la diferencia cuadrada sumada entre matrices es de aproximadamente 2.45e-13. Sin embargo, si realizo una entrada de tamaño 46x46, no converge hasta la centésima o más iteraciones. Para 47x47 o más grande, las matrices nunca convergen; el error fluctúa entre 100 y 1000 para aproximadamente 100 iteraciones, y luego comienza a crecer muy rápidamente.

¿Cómo harías para tratar de depurar algo como esto? No parece haber ningún punto específico en el que se vuelva loco, y las matrices son demasiado grandes para que realmente intente hacer el cálculo a mano. ¿Alguien tiene consejos/tutoriales/heurística para encontrar bichos numéricos extraños como este?

Nunca he tratado con algo como esto antes, pero espero que algunos de ustedes tienen ...

Gracias, - Dan

[1] G. T. Wilson. "La factorización de las densidades espectrales matriciales". SIAM J. Appl. Matemáticas (Vol 23, No. 4, diciembre de 1972)

Answer 1

Recomendaría hacer esta pregunta en la lista de correo scipy-user, quizás con un ejemplo de su código. En general, las personas en la lista parecen tener mucha experiencia en computación numérica y son realmente útiles, solo seguir la lista es una educación en sí misma.

De lo contrario, me temo que no tengo ninguna idea ... Si crees que es un problema de precisión numérica/redondeo de coma flotante, lo primero que puedes probar es modificar todos los tipos hasta float128 y ver si hace alguna diferencia

Answer 2

Interval arithmetic puede ayudar, pero no estoy seguro de si el rendimiento será suficiente para realmente permitir una depuración significativa en los tamaños de matriz de su interés (hay que pensar en un par de magnitudes de magnitud de ralentización, entre sustituciones -HW ayudó a las operaciones de punto flotante "escalar" con SW "pesadas" de intervalo, y agregando las verificaciones sobre qué intervalos crecen demasiado, cuándo, dónde y por qué).