Работа с очень большими матрицами в numpy - PullRequest
Новая
       42

Работа с очень большими матрицами в numpy

1 голос
/ 10 марта 2020

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

def steady_state(matrix):
    dim = matrix.shape[0]
    q = (matrix - np.eye(dim))
    ones = np.ones(dim)
    q = np.c_[q, ones]
    qtq = np.dot(q, q.T)
    bqt = np.ones(dim)
    return np.linalg.solve(qtq, bqt)

Однако матрица, с которой я работаю, имеет около 1,5 миллиона строк и столбцов . Это не разреженная матрица; большинство записей маленькие, но ненулевые. Конечно, просто попытка построить эту матрицу приводит к ошибке памяти.

Как я могу изменить приведенный выше код для работы с огромными матрицами? Я слышал о решениях как PyTables, но я не уверен, как их применять, и я не знаю, будут ли они работать для таких задач, как np.linalg.solve.

Будучи очень новым для numpy и очень неопытным с линейным Я очень ценю пример того, что делать в моем случае. Я открыт для использования чего-то другого, кроме numpy, и даже чего-то другого, чем Python, если необходимо.

1 Ответ

0 голосов
/ 14 марта 2020

Вот несколько идей для начала:

Мы можем использовать тот факт, что любой начальный вектор вероятности будет сходиться в стационарном состоянии при эволюции во времени (предполагая, что это ergodi c, aperiodi c, регулярное, et c).

Для небольших матриц мы могли бы использовать 

def steady_state(matrix):
    dim = matrix.shape[0]
    prob = np.ones(dim) / dim
    other = np.zeros(dim)
    while np.linalg.norm(prob - other) > 1e-3:
        other = prob.copy()
        prob = other @ matrix
    return prob

(я думаю, что условием, предполагаемым функцией в вопросе, является распределение go в строках).

Теперь мы можем использовать тот факт, что умножение матриц и norm может быть выполнено чанком по чанку: 

def steady_state_chunk(matrix, block_in=100, block_out=10):
    dim = matrix.shape[0]
    prob = np.ones(dim) / dim
    error = 1.
    while error > 1e-3:
        error = 0.
        other = prob.copy()
        for i in range(0, dim, block_out):
            outs = np.s_[i:i+block_out]
            vec_out = np.zeros(block_out)
            for j in range(0, dim, block_in):
                ins = np.s_[j:j+block_in]
                vec_out += other[ins] @ matrix[ins, outs]
            error += np.linalg.norm(vec_out - prob[outs])**2
            prob[outs] = vec_out
        error = np.sqrt(error)
    return prob

Это должно использовать меньше памяти для временных, хотя вы могли бы добиться большего, используя out параметр np.matmul. Я должен добавить кое-что, чтобы иметь дело с последним срезом в каждом l oop, в случае, если dim не делится на block_*, но я надеюсь, что вы поняли.

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

...