Python левое умножение матрицы с обратной разреженной матрицей

Я пытаюсь вычислить выражение вида K = P*C.T*S^-1 (реализация фильтра Калмана)

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

я пытался использовать

import scipy.sparse.linalg as spln

self.K = self.P.dot(spln.spsolve(S.T, C).T)

Проблема в том, что spsolve ожидает, что второй аргумент будет вектором, а не матрицей.

редактировать: Уточнение, проблема может быть решена в Matlab с помощью K = P * (C/S), поэтому я ищу метод, похожий на spsolve, но который может принимать матрицу в качестве второго аргумента. Конечно, это можно сделать, разбив C на несколько векторов-столбцов c1..cn и решив задачу для каждого из них, а затем снова собрав их в матрицу, но я подозреваю, что это будет громоздко и неэффективно.

edit2&3: Размеры матриц обычно составляют около P~10⁶x10^6, S~100x100, C=100x10⁶. P диагональный и S симметричный, а C будет иметь только один элемент в строке. Он будет использоваться для реализации фильтра Калмана с использованием разреженных матриц, см.

http://en.wikipedia.org/wiki/Kalman_filter#The_Kalman_filter


python numpy scipy sparse-matrix linear-equation
person ajn    schedule 14.10.2011    source источник
comment
Не вычислив обратное, вы не сможете вычислить K. Что вы можете сделать без вычисления обратного, так это вычислить Kx для некоторого вектора x, что потребовало бы решения линейной системы.   -  person Sven Marnach    schedule 14.10.2011
comment
Я не согласен, в Matlab решение моего вопроса будет просто: K = P*(S'\C)' или, что то же самое, K = P*(C/S) Тот факт, что C является матрицей, а не вектором, не меняется рассуждение, вы делаете в соответствии с тем, что вы говорите, решая один раз для каждого столбца в C. Мой вопрос о том, что spsolve ограничивает меня тем, что C является вектором, тогда как в Matlab это также можно сделать для матриц. В зависимости от размерности матриц это все еще может быть значительно более эффективным, чем вычисление фактического обратного.   -  person ajn    schedule 14.10.2011
comment
Извините, я неявно предполагал, что все матрицы квадратные. Так почему бы просто не перебрать столбцы C и решить для каждого столбца? Поскольку вам приходится каждый раз решать линейную систему, накладные расходы на цикл будут незначительными.   -  person Sven Marnach    schedule 14.10.2011
comment
Только что увидел вашу правку. Нет, я не думаю, что перебор столбцов C будет неэффективным. Однако в зависимости от размера линейной системы рассмотрите возможность использования итеративного решателя вместо spsolve().   -  person Sven Marnach    schedule 14.10.2011
comment
Если нет лучшего метода, я обязательно попробую решить его, перебирая столбцы, но у меня есть подозрение, что это будет не так эффективно. Размеры матриц обычно составляют около P~10⁶x10^6, S~100x100, C=100x10⁶. P и S будут диагональными, а C будет иметь только один элемент в строке. Я также обновлю свой вопрос с этой информацией.   -  person ajn    schedule 14.10.2011
comment
При заданных габаритах это конечно будет малоэффективно! Но если S диагональная, почему бы просто не инвертировать ее?   -  person Sven Marnach    schedule 14.10.2011
comment
Хорошо, к сожалению, мое утверждение было неверным :( S просто симметрична, а не диагональна. Однако S имеет довольно малые размеры по сравнению с другими, поэтому я думаю, что может быть целесообразно просто инвертировать ее как плотную матрицу. Но я все еще хотел бы получить ответ на общий вопрос, кто-то указал мне, что базовый код для spsolve, похоже, решает общую проблему и что это просто оболочка scipy для ненужного ограничения параметров.   -  person ajn    schedule 14.10.2011
comment
Общее эмпирическое правило должно быть таким: если количество столбцов C меньше размерности S, выполните итерацию по столбцам C, в противном случае инвертируйте S. И если S составляет только 100x100, не беспокойтесь об инвертировании - это займет миллисекунду или около того.   -  person Sven Marnach    schedule 14.10.2011

Ответы (1)


arrow_upward
2
arrow_downward

Как обходной путь может сделать

import numpy as np
from scipy.sparse.linalg import splu

def spsolve2(a, b):
    a_lu = splu(a)
    out = np.empty((A.shape[1], b.shape[1]))
    for j in xrange(b.shape[1]):
        out[:,j] = a_lu.solve(b[:,j])
    return out

self.K = self.P.dot(spsolve2(S.T, C).T)

Но да, это баг, что spsolve не принимает матрицы.

Однако, поскольку ваш S не очень велик, вы также можете использовать плотную инверсию.

person pv.    schedule 18.10.2011