Eine große Anzahl von kleinen linearen Systemen lösen

Question

Ich muss eine große Anzahl von 3x3 symmetrischen, postive-definite Systemen mit Python lösen. Bisher habe ich

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

Dies erzeugt das richtige Ergebnis, aber auch ziemlich langsam ist, wenn n groß ist. (Naja ... Ja.) Vielleicht ist 3x3 keine Problemgröße, wo das Aufrufen von solve() Sinn macht.

Irgendwelche Hinweise?

Answer 1

In NumPy 1.8 und später, numpy.linalg.solve actually broadcasts. Für numpy.linalg.solve(a, b), wenn b.ndim == a.ndim - 1, wird eine Broadcast-Matrix-Vektorlösung durchgeführt; Andernfalls wird eine Matrix-Matrix-Lösung gesendet. (Dieses Entscheidungskriterium nicht dokumentiert ist,. Ich hatte an der Quelle sehen)

Wenn Sie effizient einen Stapel von A s und rhs s konstruieren können, können Sie solve einmal anrufen und eine Python-Schleife vermeiden.