Лучший способ вычислить фундаментальную матрицу поглощающей цепи Маркова?

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

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

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

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

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

Каков наилучший способ вычисления строки (или всех строк) для инверсии разреженной матрицы?

ИЛИ

Какой лучший способ вычислить строку фундаментальной матрицы большой поглощающей цепи Маркова?

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

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

Ответ 1

Предполагая, что то, что вы пытаетесь сделать, - это ожидаемое количество шагов перед поглощением, уравнение от "Конечный Марков" Цепи "(Kemeny and Snell), который воспроизводится в Википедии:

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

Переупорядочивание:

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

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

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

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

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

Ответ 2

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

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