Решение СЛАУ с симметричной разреженной матрицей

В этой статье мы будем рассматривать решения СЛАУ вида Ax = b, где A — симметричная разреженная матрица. Такие матрицы появляются, например, при решении задач методом наименьших квадратов. Для симметричных СЛАУ разработаны специальные методы, такие, как метод Холецкого и LDL^T разложение ( см. [1][2] ). Так как первый из них применим к более узкому классу матриц, чем второй, поэтому далее будем рассматривать только LDL^T разложение, хотя выводы этой статьи применимы к обоим методам.

При LDL^T разложении матрица A представляется в виде произведения LDL^T, где L — нижняя треугольная матрица с единицами на главной диагонали, D — диагональная матрица,
L^T— транспонированная матрица L. Затем последовательно решаются уравнения Ly = b и DL^Tx = y. При этом затрачивается около n³/6 умножений ( где n — размерность матрицы А ),
что в 2 раза меньше, чем в методе Гаусса. Это неудивительно, так как при большой размерности в симметричной матрице содержится приблизительно в 2 раза меньше информации, чем в матрице общего вида.
Далее, чтобы уменьшить количество вычислений, воспользуемся разреженностью матрицы А. Элементы матриц D и L вычисляются по формулам:

Из этих формул следует, что если в какой-то строке матрицы А первые k элементов были нулями, то и в матрице L в той же строке первые k элементов тоже будут нулями. А это значит, что эти элементы вычислять не надо, и, если их много, то получается большая экономия в вычислениях. Это основная идея этой статьи. Назовём эти нули — «лидирующими нулями».

Ниже приведена программа на языке С++, которая решает СЛАУ с учётом лидирующих нулей. Это реальная программа с проверками на переполнение и исчезновение порядка
для чисел с плавающей точкой. Матрица А поступает на вход программы в компактном виде. Во-первых, опускаются лидирующие нули, а во-вторых, указываются элементы только до диагонального, так как матрица симметричная. Новые матрицы L и D записываются на место исходной матрицы А, причём для матрицы D записываются обратные диагональные элементы.
Входные параметры функции slu_LDLt:
n — размерность матрицы А,
m — массив с количеством лидирующих нулей по строкам,
a — матрица А записанная в виде массива по строкам фрагментами от первого ненулевого до диагонального элемента,
b — массив свободных членов,
x — массив неизвестных.

typedef unsigned int nat;  template<class T> inline const T & _max ( const T & i1, const T & i2 ) {     return i1 > i2 ? i1 : i2; }  bool slu_LDLt ( nat n, const nat * m, double * a, const double * b, double * x ) {     nat i, j, k, l, nn = 0; // Заполнение индексов диагональных элементов     std::vector<nat> di ( n );     for ( i = 0; i < n; ++i )     {         if ( m[i] > i )             return false;         di[i] = nn += i - m[i];         ++nn;     } // Построение матриц L и D     for ( i = l = 0; i < n; ++i )     {         const nat mi = m[i];         const nat ni = di[i] - i;         const double * ai = a + ni;         for ( j = mi; j <= i; ++j )         {             double * aj = a + di[j] - j;             double x = ai[j];             if ( i == j )             {                 for ( k = mi; k < j; ++k )                 {                     double & ajk = aj[k];                     const double y = ajk;                     ajk *= a[di[k]];                     x -= y * ajk;                 }                 const double ax = fabs ( x );                 if ( ax < 1e-100 )                 {                     return false;                 }                 a[l] = ax > 1e100 ? 0 : 1 / x;             }             else             {                 for ( k = _max ( mi, m[j] ); k < j; ++k ) x -= ai[k] * aj[k];                 const double ax = fabs ( x );                 if ( ax > 1e100 )                 {                     return false;                 }                 a[l] = ax < 1e-100 ? 0 : x;             }             ++l;         }     } // Решение системы Ly = b     for ( i = j = 0; i < n; ++i, ++j )     {         double t = b[i];         for ( k = m[i]; k < i; ++k ) t -= a[j++] * x[k];         x[i] = t;     } // Решение системы DLtx = y;     for ( i = 0; i < n; ++i ) x[i] *= a[di[i]];     for ( l = nn; --i > 0; )     {         const double t = x[i];         const double * p = a + ( l -= i + 1 );         for ( l += k = m[i]; k < i; ++k ) x[k] -= p[k] * t;     }     return true; } 

В случае отсутствия лидирующих нулей — это будет обычная программа для LDL^T разложения. А если их будет много, то время работы программы резко уменьшится по сравнению с обычной. В следующей статье будет показано, как можно увеличить количество лидирующих нулей, переставляя столбцы и строки исходной матрицы А.

Литература.

1. Уилкинсон, Райнш. «Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра».

2. Голуб, Ван Лоун. «Матричные вычисления».