Ich möchte Ax = b lösen, wobei A eine sehr große quadratische positiv definite symmetrische Blockmatrix ist und x und b Vektoren sind. Wenn ich groß sage, meine ich für eine nxn Matrix eine n so groß wie 300,000.Lösen großer linearer Systeme mit Block-Sparse-Matrizen
Hier ist ein Beispiel für eine viel kleinere, aber repräsentative Matrix, die ich lösen möchte.
Und hier ist die gleiche Matrix zeigen, gezoomt, dass es von Blöcken aus dichtem Matrizes besteht.
ich vorher (siehe here, here und wie here weit zurück) verwendet Cholesky Solver Eigen, die für n<10000 aber mit mit n=300000 der Cholesky-Solver langsam ist zu fein gearbeitet. Allerdings habe ich nicht ausgenutzt, dass ich eine Blockmatrix habe. Offensichtlich existieren Algorithmen zum Lösen von Sparse-Block-Matrizen (z. B. block cholesky factorization).
Ich möchte speziell wissen, ob Eigen Algorithmen mit Faktorisierung oder iterativen Methoden für spärliche dichte Blockmatrizen optimiert hat, die ich verwenden kann?
Können Sie auch andere Algorithmen vorschlagen, die ideal sind, um meine Matrix zu lösen? Ich meine, so weit ich zumindest für Faktorisierung verstehe, eine Permutationsmatrix zu finden, ist NP schwer, so viele verschiedene heuristische Methoden existieren und soweit ich sagen kann, bauen Menschen eine Intuition verschiedener Matrixstrukturen auf (zB banded matrix) und welche Algorithmen am besten funktionieren Sie. Ich habe diese Intuition noch nicht.
Ich bin derzeit in der Suche mit einer conjugate gradient method. Ich habe dies selbst in C++ implementiert, aber ich nutze immer noch nicht die Tatsache, dass die Matrix eine Blockmatrix ist.
//solve A*rk = xk
//Eigen::SparseMatrix<double> A;
//Eigen::VectorXd rk(n);
Eigen::VectorXd pk = rk;
double rsold = rk.dot(rk);
int maxiter = rk.size();
for (int k = 0; k < maxiter; k++) {
Eigen::VectorXd Ap = A*pk;
double ak = rsold /pk.dot(Ap);
xk += ak*pk, rk += -ak*Ap;
double rsnew = rk.dot(rk);
double xlength = sqrt(xk.dot(xk));
//relaxing tolerance when x is large
if (sqrt(rsnew) < std::max(std::min(tolerance * 10000, tolerance * 10 * xlength), tolerance)) {
rsold = rsnew;
break;
}
double bk = rsnew/rsold;
pk *= bk;
pk += rk;
rsold = rsnew;
}
Möchten Sie das wirklich selbst implementieren? Vielleicht solltest du ein paar verfügbare Bibliotheken ausprobieren, die solche Sachen machen? Z.B. [elemental] (http://libelemental.org/). – sascha
@sascha, Ich bin glücklich, Bibliotheken zu verwenden, ich mache das bereits mit Eigen (außer einer benutzerdefinierten konjugierten Gradientenmethode), aber ich würde andere berücksichtigen, wenn Eigen für diesen Fall nicht ausreicht. –
Eigen klingt für mich wie eine einfache Lineare Algebra, ohne (komplexere; strukturverwendende) Optimierungsalgorithmen an sich. Aber ich habe es nie benutzt. Obwohl ich überhaupt keine Erfahrung mit * elemental * habe, denke ich, dass das, was du tust, irgendwie möglich ist (und wahrscheinlich hoch optimiert). Leider gibt es keine Hinweise mehr, die ich dir geben könnte. Vielleicht sehen Sie sich [this] (http://libelemental.org/documentation/dev/optimization/models.html) einen Teil der Dokumente an. Viel Glück! – sascha