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).Estrategias para depurar problemas de estabilidad numérica?
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)
¿Qué quiere decir con "solo converge para matrices de tamaño 45x45?" ¿También fallan también las matrices menores de 45x45? – badp
No, lo siento, editaré la publicación. Converge con éxito para el tamaño 45x45 y más pequeño – Dan