Стратегии для отладки проблем с числовой стабильностью?

Question

Я пытаюсь написать реализацию алгоритма факторизации спектральной плотности Вильсона [1] для Python. Алгоритм итеративно факторизует матричную функцию [QxQ] в ее квадратный корень (это своего рода расширение квадратного корневого искателя Ньютона-Рафсона для спектральных матриц плотности).

Проблема в том, что моя реализация только сходится для матриц размером 45x45 и меньше. Таким образом, после 20 итераций суммарный квадрат разницы между матрицами составляет около 2.45e-13. Однако, если я делаю ввод размера 46x46, он не сходится до 100-й или около того итерации. Для 47x47 или больше матрицы никогда не сходятся; ошибка колеблется от 100 до 1000 для примерно 100 итераций, а затем начинает расти очень быстро.

Как бы вы попытались отладить что-то вроде этого? Кажется, нет какой-либо конкретной точки, в которой она сошла с ума, и матрицы слишком велики для меня, чтобы на самом деле попытаться выполнить расчет вручную. У кого-нибудь есть советы/учебники/эвристики для поиска причудливых числовых ошибок вроде этого?

Я никогда не имел дело с чем-нибудь, как это раньше, но я надеюсь, что некоторые из вас ...

Спасибо, - Dan

[1] Г. Т. Уилсон. «Факторизация матричных спектральных плотностей». SIAM J. Appl. Math (Vol 23, No. 4, Dec. 1972)

Answer 1

Я бы порекомендовал задать этот вопрос в списке рассылки , возможно, на примере вашего кода. Как правило, люди в списке, как представляется, имеют большой опыт работы с численными вычислениями и действительно полезны, просто следуя за списком, это образование само по себе.

В противном случае, я боюсь, что у меня нет идей ... Если вы считаете, что это числовая проблема с округлением с точностью до точки с плавающей точкой, первое, что вы могли бы попробовать, это показать все типы dtypes до float128 и посмотреть если имеет значение.

Answer 2

Interval arithmetic может помочь, но я не уверен, что производительность будет достаточной, чтобы фактически позволить осмысленную отладку в размерах матрицы вашего интереса (вам нужно рассчитывать на пару порядков замедленного действия, что между заменой высоко -HW - помогал «скалярным» операциям с плавающей запятой с SW-тяжелыми «интервалами» и добавлял проверки того, какие интервалы становятся слишком широкими, когда, где и почему).