Ich muss eine große, dichte Matrix invertieren, von der ich hoffte, Scipys gmres zu verwenden. Glücklicherweise folgt die dichte Matrix A einem Muster und ich muss die Matrix nicht im Speicher speichern. Die Klasse LinearOperator ermöglicht es uns, ein Objekt zu konstruieren, das als Matrix für GMRES fungiert und direkt das Matrixvektorprodukt A*v berechnen kann. Das heißt, wir schreiben eine Funktion mv(v) , die einen Vektor v als Eingabe verwendet und mv(v) = A*v zurückgibt. Dann können wir die Klasse LinearOperator verwenden, um A_LinOp = LinearOperator(shape = shape, matvec = mv) zu erstellen. Wir können den linearen Operator in den Befehl Scipy gmres setzen, um die Matrix-Vektorprodukte auszuwerten, ohne A vollständig in den Speicher laden zu müssen.
Die Dokumentation für LinearOperator finden Sie hier: LinearOperator Dokumentation .
Hier ist mein Problem: Um die Routine zu schreiben, um das Matrixvektorprodukt mv(v) = A*v zu berechnen, brauche ich einen anderen Eingabevektor C . Die Einträge in A haben die Form A[i,j] = f(C[i] - C[j]) . Was ich wirklich will, ist, dass mv aus zwei Eingaben besteht, einer festen Vektoreingabe C und einer Variableneingabe v , die wir A*v berechnen wollen.
MATLAB hat eine ähnliche Einstellung, wo x = gmres(@(v) mv(v,C),b) geschrieben wird, wobei b die rechte Seite des Problems ist Ax = b , und mv ist die Funktion, die als variable Eingabe v verwendet, die wir wollen A*v und C zu berechnen ist der feste, bekannte Vektor, den wir für die Assemblierung von A benötigen.
Mein Problem ist, dass ich nicht herausfinden kann, wie ich der LinearOperator -Klasse erlauben kann, zwei Eingaben zu akzeptieren, eine Variable und eine "feste", wie ich es in MATLAB kann.
Gibt es eine Möglichkeit, die analoge Operation in SciPy durchzuführen? Wenn jemand eine bessere Möglichkeit kennt, eine große, dichte Matrix (50000, 50000) zu invertieren, wo die Einträge einem Muster folgen, würde ich mich über Vorschläge freuen.
Danke!
EDIT: Ich hätte diese Information eigentlich angeben sollen. Die Matrix ist tatsächlich (in Blockform) [A C; C^T 0] , wobei A ist N x N ( N groß) und C ist N x 3 und 0 ist 3 x 3 und C^T ist die Transponierung von C . Dieses Array C ist dasselbe Array wie das oben erwähnte Array. Die Einträge von A folgen einem Muster A[i,j] = f(C[i] - C[j]) .
Ich schrieb mv(v,C) , um Zeile für Zeile zu gehen, konstruiere A*v[i] für i=0,N , indem ich sum% f(C[i]-C[j)*v[j] (eigentlich mache ich numpy.dot(FC,v) wo FC[j] = f(C[i]-C[j]) , was gut funktioniert). Am Ende werden die Berechnungen für die C^T Zeilen ausgeführt. Ich hatte gehofft, die große for-Schleife mit einem multiprocessing -Aufruf zu ersetzen, um die for-Schleife zu parallelisieren, aber das ist eine zukünftige Sache, die man beachten sollte. Ich werde mich auch mit Cython befassen, um die Berechnungen zu beschleunigen.
Das ist sehr spät, aber wenn Sie immer noch interessiert sind ...
Ihre A-Matrix muss einen sehr niedrigen Rang haben, da es sich um eine nicht linear transformierte Version einer Rang-2-Matrix handelt. Plus es ist symmetrisch. Das bedeutet, dass es invers ist, die abgeschnittene Eigenwertzerlegung mit etwa 5 Eigenwerten zu berechnen: A = U * S * U ', dann invertiere das: A ^ -1 = U * S ^ -1 * U'. S ist diagonal, so dass dies kostengünstig ist. Sie können die abgeschnittene Eigenwertzerlegung mit ach erhalten.
Das kümmert sich um A. Dann für den Rest: Verwenden Sie die Blockmatrix-Inversionsformel . Sieht böse aus, aber ich wette 100,000,000 Preussische Francs, dass es 50x schneller ist als die direkte Methode, die Sie benutzten.
Tags und Links python matlab numpy matrix matrix-multiplication