Решение больших линейных систем с блочными разреженными матрицами

Question

Я хочу решить Ax = b где A - очень большая квадратная положительно определенная симметричная матричная матрица и x и b - векторы. Когда я говорю «большой», я имею в виду матрицу nxnn размером 300,000.

Вот пример гораздо меньшей, но представительной матрицы, которую я хочу решить.

А вот та же матрица увеличено, показывая, что она состоит из блоков плотных матриц,.

Я ранее (см here, here, и еще в here) используется Чолески решатель Эйген, который работал отлично для n<10000 но с n=300000 решатель Cholesky слишком медленно. Однако я не воспользовался тем, что у меня есть блок-матрица. По-видимому, существуют алгоритмы для решения разреженных блочных матриц (например, block cholesky factorization).

Я хотел бы знать конкретно, если Eigen оптимизировал алгоритмы, используя факторизацию или итеративные методы, для разреженных плотных матричных матриц, которые я могу использовать?

Также вы можете предложить другие алгоритмы, которые могут быть идеальными для решения моей матрицы? Я имею в виду, насколько я понимаю, по крайней мере, для факторизации поиск матрицы перестановок является NP трудным, так что существует множество различных эвристических методов, и насколько я могу судить, что люди создают интуицию различных матричных структур (например, banded matrix) и какие алгоритмы работают лучше всего с их. У меня еще нет этой интуиции.

В настоящее время я изучаю conjugate gradient method. Я сам реализовал это на C++, но все еще не воспользовался тем, что матрица является блочной матрицей.

//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; 
} 

Answer 1

Я думаю, что ARPACK был разработан, чтобы быть эффективным для задач, как найти решение системы уравнений Ах = Ь, когда А sparce, как в вашем случае. Вы можете найти исходный код для ARPACK here.

Answer 2

Отъезд SuiteSparse, http://faculty.cse.tamu.edu/davis/suitesparse.html Автор, Тим Дэвис, известен в области редких матриц. Он также получил награду от Google за качество кода. Выполнение его кода также выдающееся.

Cheers