Я пытаюсь написать реализацию алгоритма факторизации спектральной плотности Уилсона [1] для Python. Алгоритм итеративно разлагает матричную функцию [QxQ] на ее квадратный корень (это своего рода расширение искателя квадратного корня Ньютона-Рафсона для матриц спектральной плотности).
Проблема в том, что моя реализация сходится только для матриц размером 45x45 и меньше. Таким образом, после 20 итераций суммарная квадратичная разница между матрицами составляет около 2,45e-13. Однако, если я сделаю ввод размером 46x46, он не сходится до 100-й итерации или около того. Для размеров 47x47 и больше матрицы никогда не сходятся; ошибка колеблется между 100 и 1000 примерно на 100 итерациях, а затем начинает очень быстро расти.
Как бы вы попытались отладить что-то подобное? Кажется, нет какой-то конкретной точки, в которой он сходит с ума, и матрицы слишком велики, чтобы я мог попытаться выполнить расчет вручную. Есть ли у кого-нибудь советы/учебники/эвристика для поиска странных числовых ошибок, подобных этой?
Я никогда не сталкивался с чем-то подобным раньше, но я надеюсь, что некоторые из вас...
Спасибо, - Дэн
[1] Г. Т. Уилсон. «Факторизация матричных спектральных плотностей». СИАМ Дж. Заявл. Математика (Том 23, № 4, декабрь 1972 г.)
