Обработка ошибок округления при умножении матриц в numpy

TL;DR: умножение матриц для матрицы перехода состояний должно сохранять норму, но np.matmul не сохраняет норму. Как я могу это исправить? Есть ли лучший модуль Python для этого?

У меня есть правильная матрица перехода состояния, A, т. Е. s (t) A (tau) = s (t + tau)

где s(t) — матрица-столбец, сумма которой равна 1. Кроме того, мы знаем, что каждая строка матрицы A также дает в сумме 1.

Мы знаем, что A^n также является правой матрицей перехода состояния для любого n в натуральных числах.

Один из способов найти стационарное распределение состоит в том, чтобы вычислить A^n, когда n стремится к бесконечности. Следующий фрагмент вычисляет A^(2^n):

    def p_multiplier(A,n):
        B=copy.deepcopy(A)
        for i in range(n):
            B=numpy.matmul(B,B)
        return B[0]

Это работает нормально, но по какой-то причине суммирование строк начинает не складывать до 1, т.е. numpy начинает пропускать норму A[i].

Причина, вероятно, в ошибках округления. Одним из быстрых решений может быть принудительное сохранение нормы каждой строки на каждой итерации:

def p_multiplier(A,n):
    B=copy.deepcopy(A)
    for i in range(n):
        B=np.matmul(B,B)
        for j in range(len(B)):
            B[j]=B[j]/sum(B[j])
    return B[0]

Насколько научно это исправление? С этим лучше справляются scipy или mpmath?

Редактировать: чтобы полностью выполнить MWE, вы можете использовать следующую матрицу перехода состояний для A:

A=np.asarray([[0.76554539,0.13929202,0,0,0.09516258],[0.04877058,0.76996018,0.18126925,0,0],[0,0.09516258,0.76554539,0.13929202,0],[0,0,0.13929202,0.67943873,0.18126925],[0.09516258,0,0,0.04877058,0.85606684]])

Кроме того, если есть латексная поддержка stackoverflow, это было бы здорово.

C.Koca 09.04.2022 15:09

Пожалуйста, опубликуйте примеры массивов для s и A, а также ожидаемый результат.

AJH 09.04.2022 15:15

@AJH s является произвольным и не встречается в устойчивых распределениях, потому что независимо от вашего начального состояния вы должны достичь одного и того же устойчивого распределения. Я все же добавлю A.

C.Koca 09.04.2022 15:17

Ваша входная матрица A не имеет суммы, близкой к 1 для каждой строки, потому что не хватает цифр. Не могли бы вы поделиться значениями матрицы с большим количеством цифр, например ~ 14 цифр вместо 8 ~ 9?

Jérôme Richard 09.04.2022 16:12

Расчеты с поплавками не являются «точными». Люди часто жалуются на почти нулевые значения в matmul(A, inv(A))

hpaulj 09.04.2022 16:14
Анализ настроения постов в Twitter с помощью Python, Tweepy и Flair
Анализ настроения постов в Twitter с помощью Python, Tweepy и Flair
Анализ настроения текстовых сообщений может быть настолько сложным или простым, насколько вы его сделаете. Как и в любом ML-проекте, вы можете выбрать...
7 лайфхаков для начинающих Python-программистов
7 лайфхаков для начинающих Python-программистов
В этой статье мы расскажем о хитростях и советах по Python, которые должны быть известны разработчику Python.
Установка Apache Cassandra на Mac OS
Установка Apache Cassandra на Mac OS
Это краткое руководство по установке Apache Cassandra.
Сертификатная программа "Кванты Python": Бэктестер ансамблевых методов на основе ООП
Сертификатная программа "Кванты Python": Бэктестер ансамблевых методов на основе ООП
В одном из недавних постов я рассказал о том, как я использую навыки количественных исследований, которые я совершенствую в рамках программы TPQ...
Создание персонального файлового хранилища
Создание персонального файлового хранилища
Вы когда-нибудь хотели поделиться с кем-то файлом, но он содержал конфиденциальную информацию? Многие думают, что электронная почта безопасна, но это...
Создание приборной панели для анализа данных на GCP - часть I
Создание приборной панели для анализа данных на GCP - часть I
Недавно я столкнулся с интересной бизнес-задачей - визуализацией сбоев в цепочке поставок лекарств, которую могут просматривать врачи и...
1
5
30
1
Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Проблема возникает из-за того, что операция не численно стабильный. В результате это быстро расходятся (экспоненциально) до 0 или бесконечности даже при относительно небольших значениях n, таких как 70. Вы можете использовать метод диагонализации на основе собственных значений (см. здесь для получения дополнительной информации), который гораздо более численно стабилен.

def stable_multiplier(A, n):
    d, P = numpy.linalg.eig(A)
    return (P @ numpy.diag(d**n) @ numpy.linalg.inv(P))[0].real

Количество потерянных (десятичных) цифр находится в O(log10(n)), поэтому ошибка линейна относительно n. Обратите внимание, что если n большое, например> 1000, вам нужно быть очень осторожным в отношении численной стабильности всего алгоритма, особенно если это итеративный процесс (с количеством итераций). Для получения дополнительной информации обратитесь к темам это и это.

Спасибо! Это фактически решает мой вопрос. Теперь я помню, что раньше использовал собственное разложение для этой конкретной цели. Извините за плохой вопрос :)

C.Koca 09.04.2022 21:13

Другие вопросы по теме