Вопрос по python, algorithm, markov-chains, sparse-matrix, math – Лучший способ вычислить фундаментальную матрицу поглощающей цепи Маркова?

11

У меня очень большая поглощающая цепь Маркова (масштабируется до размера задачи - от 10 штатов до миллионов), которая очень разрежена (большинство штатов может реагировать только на 4 или 5 других штатов).

Мне нужно рассчитать одну строку фундаментальной матрицы этой цепочки (средняя частота каждого состояния дана для одного исходного состояния).

Обычно я делаю это путем вычисления(I - Q)^(-1), но я не смог найти хорошую библиотеку, в которой реализован обратный алгоритм разреженной матрицы! Я видел несколько работ, большинство из которых P.h.D. уровень работы.

Большинство моих результатов в Google указывают мне на сообщения о том, что нельзя использовать матрицу, обратную при решении линейных (или нелинейных) систем уравнений ... Я не считаю это особенно полезным. Похоже ли вычисление фундаментальной матрицы на решение системы уравнений, и я просто не знаю, как выразить одно в форме другого?

Итак, я задаю два конкретных вопроса:

What's the best way to calculate a row (or all the rows) of the inverse of a sparse matrix?

ИЛИ ЖЕ

What's the best way to calculate a row of the fundamental matrix of a large absorbing Markov chain?

Решение на Python было бы замечательно (так как мой проект все еще является проверкой концепции), но если бы мне пришлось испачкать руки хорошим старым приложением. Fortran или C, это не проблема.

Редактировать: я только что понял, что обратная B матрицы A может быть определена как AB = I, где I - единичная матрица. Это может позволить мне использовать некоторые стандартные разреженные матричные решатели для вычисления обратного ... Я должен убежать, поэтому не стесняйтесь завершать мой ход мыслей, который, как я начинаю думать, может потребовать только действительно элементарной матрицы имущество...

Я работал над некоторыми вещами на PGM и задавался вопросом, есть ли способ вычислить это вообще - хотя понятия нет для разреженной матрицы, так что удачи! argentage
Если вы хотите решение Python, пожалуйста, пометьте егоpython, Существуют и другие обмены стеков, которые могут быть более или менее полезными. Jared Farrish

Ваш Ответ

2   ответа
2

Предполагая, что вы пытаетесь сделать это,ожидаемое количество шагов до поглощенияуравнение из «конечных цепей Маркова»; (Кемени и Снелл), который воспроизводится в Википедии:

t=N1

Или расширение фундаментальной матрицы

t=(I-Q)^-1 1

Перегруппировка:

(I-Q) t = 1

Который в стандартном формате для использования функций для решения систем линейных уравнений

A x = b

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

import networkx as nx
import numpy

def example(n):
    """Generate a very simple transition matrix from a directed graph
    """
    g = nx.DiGraph()
    for i in xrange(n-1):
        g.add_edge(i+1, i)
        g.add_edge(i, i+1)
    g.add_edge(n-1, n)
    g.add_edge(n, n)
    m = nx.to_numpy_matrix(g)
    # normalize rows to ensure m is a valid right stochastic matrix
    m = m / numpy.sum(m, axis=1)
    return m

Представление двух альтернативных подходов для расчета количества ожидаемых шагов.

def expected_steps_fundamental(Q):
    I = numpy.identity(Q.shape[0])
    N = numpy.linalg.inv(I - Q)
    o = numpy.ones(Q.shape[0])
    numpy.dot(N,o)

def expected_steps_fast(Q):
    I = numpy.identity(Q.shape[0])
    o = numpy.ones(Q.shape[0])
    numpy.linalg.solve(I-Q, o)

Выбираем пример, который является достаточно большим, чтобы продемонстрировать типы проблем, которые возникают при вычислении фундаментальной матрицы:

P = example(2000)
# drop the absorbing state
Q = P[:-1,:-1]

Производит следующие тайминги:

%timeit expected_steps_fundamental(Q)
1 loops, best of 3: 7.27 s per loop

А также:

%timeit expected_steps_fast(Q)
10 loops, best of 3: 83.6 ms per loop

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

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

3

Причина, по которой вы получаете совет не использовать обратные матрицы для решения уравнений, заключается в численной устойчивости. Когда ваша матрица имеет собственные значения, которые равны нулю или близки к нулю, у вас возникают проблемы либо из-за отсутствия обратной (если ноль), либо из-за числовой устойчивости (если она близка к нулю). Таким образом, способ решения проблемы состоит в использовании алгоритма, который не требует наличия обратного. Решение заключается в использованииГауссово исключение, Это не обеспечивает полного обратного, а скорее приводит вас к форме ряд-эшелон, обобщению верхней треугольной формы. Если матрица обратима, то последняя строка матрицы результата содержит строку, обратную. Так что просто сделайте так, чтобы последний ряд, который вы исключаете, был желаемым.

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

Похожие вопросы