Решение большого числа небольших линейных систем

Question

Мне нужно решить большое количество симметричных, поведенчески определенных систем 3x3 с Python. До сих пор я сделал

res = numpy.zeros(n) 
for k, obj in enumerate(data_array): 
    # construct A, rhs, idx from obj 
    res[idx] += numpy.linalg.solve(A, rhs) 

Это производит правильный результат, однако также довольно медленно, если n велико. (Ну ... да.) Возможно, 3x3 не является проблемой, где вызов solve() имеет смысл.

Любые подсказки?

Answer 1

В NumPy 1.8 и более поздних версиях, numpy.linalg.solve actually broadcasts. Для numpy.linalg.solve(a, b), если b.ndim == a.ndim - 1, он будет выполнять широковещательный матричный вектор; в противном случае это приведет к разрешению широковещательной матричной матрицы. (Этот критерий решение не документированы;. Я должен был смотреть на источник)

Если вы можете эффективно построить стек A с и rhs с, вы можете вызвать solve раз и избежать цикла Python.