수치 안정성 문제를 디버깅하기위한 전략은 무엇입니까?

Question

파이썬 용 윌슨의 스펙트럼 밀도 인수 분해 알고리즘 [1]을 구현하려고합니다. 알고리즘은 [QxQ] 행렬 함수를 그것의 제곱근으로 반복적으로 인수 분해합니다 (이는 스펙트럼 밀도 행렬에 대한 Newton-Raphson 제곱근 파인더의 확장입니다).

문제는 내 구현이 크기가 45x45 이하인 행렬에만 수렴한다는 것입니다. 따라서 20 회 반복 한 결과 행렬의 제곱 된 차이는 약 2.45e-13입니다. 그러나 크기를 46x46으로 입력하면 100 번 정도까지 수렴하지 않습니다. 47x47 이상에서는 행렬이 수렴하지 않습니다. 오류는 약 100 반복에 대해 100과 1000 사이에서 변동하고 매우 빠르게 증가하기 시작합니다.

어떻게 이런 식으로 디버깅하려고합니까? 그것이 미쳐가는 특정 시점이없는 것처럼 보이며 행렬이 실제로 계산을 시도하기에 너무 큽니다. 누구나 이런 기괴한 수치 버그를 찾기위한 팁/자습서/발견 적 방법이 있습니까? 댄

[1] G. T. 윌슨 -

...

감사합니다, 을 내가 전에 이런 일 처리 적이하지만 난 당신의 일부가 바라고 있어요. "Matricial 스펙트럼 밀도의 인수 분해 (Factorization of Matricial Spectral Density)". SIAM J. Appl. Math (Vol 23, No. 4, 1972 년 12 월)

Answer 1

scipy-user 메일 링리스트에이 질문을하는 것이 좋습니다. 예를 들어 코드 예를 들어 보겠습니다. 일반적으로 목록에있는 사람들은 수치 계산에 고도의 경험이있는 것으로 보이며 실제로 도움이됩니다. 목록 바로 다음에 교육 자체가 있습니다.

그렇지 않으면 생각이별로 없습니다 ... 숫자 정밀도/부동 소수점 반올림 문제라고 생각하면 가장 먼저 시도 할 수있는 것은 모든 dtyp을 float128까지 늘려보고 차이가 있다면.

Answer 2

Interval arithmetic도 도움이 될 수 있습니다. 성능이 실제로 관심있는 매트릭스 크기에서 의미있는 디버깅을 허용하기에 충분할 지 확신하지 못합니다. (몇 단계의 속도 저하를 예상해야합니다. - SW- 무거운 "간격"객체로 "스칼라"부동 소수점 연산을 돕고 간격이 너무 넓게, 언제, 어디서, 왜 증가하는지에 대한 검사를 추가했습니다.