22.3 Итерационные методы для разреженных матриц

Операторы левого деления \ и правого деления /, рассмотренные в предыдущем разделе, используют прямые методы решения для решения линейного уравнения вида x = A \ b или x = b / A. Octave также включает ряд функций для решения разреженных линейных уравнений с помощью итерационных методов.

: x = pcg (A, b, tol, maxit, m1, m2, x0, …)

: x = pcg (A, b, tol, maxit, M, [], x0, …)

: [x, flag, relres, iter, resvec, eigest] = pcg (A, b, …)

Решает систему линейных уравнений A * x = b с помощью итерационного метода сопряжённых градиентов с предобуславливанием.

Входные аргументы:

• A — матрица линейной системы; она должна быть квадратной. В качестве A можно передать матрицу, дескриптор функции или встроенную функцию Afcn такую, что Afcn(x) = A * x. Дополнительные параметры для Afcn можно передать после x0.

  A должна быть эрмитовой и положительно определённой (HPD). Если pcg обнаруживает, что A не является положительно определённой, выводится предупреждение и устанавливается выходной параметр flag.

• b — вектор правой части.
• tol — требуемая относительная точность для ошибки невязки, b - A * x. Итерации прекращаются, если norm (b - A * x) ≤ tol * norm (b). Если tol не указан или пуст, используется точность 1e-6.
• maxit — максимально допустимое число итераций; если maxit не указан или пуст, используется значение 20.
• m — HPD-матрица предобуславливания. Для любого разложения m = p1 * p2 такого, что inv (p1) * A * inv (p2) является HPD, метод сопряжённых градиентов формально применяется к линейной системе inv (p1) * A * inv (p2) * y = inv (p1) * b, где x = inv (p2) * y (раздельное предобуславливание). На практике на каждой итерации метода сопряжённых градиентов решается линейная система с матрицей m с помощью mldivide. Если доступна определённая факторизация m = m1 * m2 (например, неполная факторизация Холецкого матрицы a), можно передать две матрицы m1 и m2, а соответствующие линейные системы будут решаться с помощью оператора mldivide. Обратите внимание, что правильный выбор предобуславливателя может значительно повысить общую производительность метода. Вместо матриц m1 и m2 пользователь может передать две функции, возвращающие результат применения обратных матриц m1 и m2 к вектору. Если m1 не указана или пуста [], предобуславливание не применяется. Если факторизация m недоступна, m2 можно не указывать или оставить равной [], а входную переменную m1 можно использовать для передачи предобуславливателя m.
• x0 — начальное приближение. Если x0 не указано или пусто, функция по умолчанию устанавливает x0 равным нулевому вектору.

Аргументы, следующие за x0, рассматриваются как параметры и соответствующим образом передаются любой из функций (A, m1 или m2), переданных в pcg. Дополнительные сведения см. в примерах ниже.

Выходные аргументы:

• x — вычисленное приближение к решению A * x = b. Если алгоритм не сошёлся, x соответствует итерации с минимальной невязкой.
• flag сообщает о сходимости:
  • 0: Алгоритм сошёлся с заданной точностью.
  • 1: Алгоритм не сошёлся и достиг максимального числа итераций.
  • 2: Матрица предобуславливателя вырождена.
  • 3: Алгоритм застопорился, то есть абсолютное значение разности между текущей итерацией x и предыдущей меньше eps * norm (x,2).
  • 4: Алгоритм обнаружил, что входная (предобусловленная) матрица не является HPD.
• relres — отношение конечной невязки к её начальному значению, измеренное в евклидовой норме.
• iter указывает номер итерации, на которой было вычислено x. Поскольку выходное значение x соответствует решению с минимальной невязкой, общее число выполненных методом итераций задаётся выражением length(resvec) - 1.
• resvec описывает историю сходимости метода. resvec (i, 1) — евклидова норма невязки, а resvec (i, 2) — норма предобусловленной невязки после (i-1)-й итерации, i = 1, 2, …, iter+1. Норма предобусловленной невязки определяется как r' * (m \ r), где r = b - A * x; см. также описание m. Если eigest не требуется, возвращается только resvec (:, 1).
• eigest возвращает оценки наименьшего eigest(1) и наибольшего eigest(2) собственных значений предобусловленной матрицы P = m \ A. В частности, если предобуславливание не используется, возвращаются оценки крайних собственных значений A. eigest(1) является завышенной оценкой, а eigest(2) — заниженной, поэтому eigest(2) / eigest(1) служит нижней границей для cond (P, 2), которая, тем не менее, теоретически должна стремиться к фактическому значению числа обусловленности.

Рассмотрим простую задачу с трёхдиагональной матрицей.

n = 10;
A = toeplitz (sparse ([1, 1], [1, 2], [2, 1], 1, n));
b = A * ones (n, 1);
M1 = ichol (A); # in this tridiagonal case it corresponds to chol (A)'
M2 = M1';
M = M1 * M2;
Afcn = @(x) A * x;
Mfcn = @(x) M \ x;
M1fcn = @(x) M1 \ x;
M2fcn = @(x) M2 \ x;

ПРИМЕР 1: Простейшее использование pcg

x = pcg (A, b)

ПРИМЕР 2: pcg с функцией, вычисляющей A * x

x = pcg (Afcn, b)

ПРИМЕР 3: pcg с матрицей предобуславливателя M

x = pcg (A, b, 1e-06, 100, M)

ПРИМЕР 4: pcg с функцией в качестве предобуславливателя

x = pcg (Afcn, b, 1e-6, 100, Mfcn)

ПРИМЕР 5: pcg с матрицами предобуславливателя M1 и M2

x = pcg (A, b, 1e-6, 100, M1, M2)

ПРИМЕР 6: pcg с функциями в качестве предобуславливателей

x = pcg (Afcn, b, 1e-6, 100, M1fcn, M2fcn)

ПРИМЕР 7: pcg с функцией, которой требуется аргумент, в качестве входных данных

function y = Ap (A, x, p) # compute A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfcn = @(x, p) Ap (A, x, p);
x = pcg (Apfcn, b, [], [], [], [], [], 2);

ПРИМЕР 8: явный пример, демонстрирующий, что pcg использует раздельный предобуславливатель

M1 = ichol (A + 0.1 * eye (n)); # factorization of A perturbed
M2 = M1';
M = M1 * M2;

## reference solution computed by pcg after two iterations
[x_ref, fl] = pcg (A, b, [], 2, M)

## split preconditioning
[y, fl] = pcg ((M1 \ A) / M2, M1 \ b, [], 2)
x = M2 \ y # compare x and x_ref

Литература:

1. C.T. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, 1995. (основной алгоритм PCG)
2. Y. Saad, Iterative Methods for Sparse Linear Systems, PWS 1996. (оценка числа обусловленности по PCG) Переработанная версия этой книги доступна в интернете по адресу https://www-users.cs.umn.edu/~saad/books.html

См. также: sparse, pcr, gmres, bicg, bicgstab, cgs.

: x = pcr (A, b, tol, maxit, m, x0, …)

: [x, flag, relres, iter, resvec] = pcr (…)

Решает систему линейных уравнений A * x = b с помощью итерационного метода сопряжённых невязок с предобуславливанием.

Входные аргументы:

• A может быть квадратной (предпочтительно разреженной) матрицей, дескриптором функции, встроенной функцией или строкой с именем функции, вычисляющей A * x. В принципе, A должна быть симметричной и невырожденной; если pcr обнаружит, что A численно вырождена, будет выведено предупреждение и установлен выходной параметр flag.
• b — вектор правой части.
• tol — требуемая относительная точность для ошибки невязки, b - A * x. Итерации прекращаются, если norm (b - A * x) <= tol * norm (b - A * x0). Если tol пуст или не указан, функция по умолчанию устанавливает tol = 1e-6.
• maxit — максимально допустимое число итераций; если для maxit указано [] или pcr содержит меньше аргументов, используется значение по умолчанию, равное 20.
• m — матрица (левого) предобуславливания, так что итерации (теоретически) эквивалентны решению методом pcr P * x = m \ b, где P = m \ A. Обратите внимание, что правильный выбор предобуславливателя может значительно повысить общую производительность метода. Вместо матрицы m пользователь может передать функцию, возвращающую результат применения обратной матрицы m к вектору (обычно это предпочтительный способ использования предобуславливателя). Если для m указано [] или m не указана, предобуславливание не применяется.
• x0 — начальное приближение. Если x0 пуст или не указан, функция по умолчанию устанавливает x0 равным нулевому вектору.

Аргументы, следующие за x0, рассматриваются как параметры и надлежащим образом передаются любой из функций (A или m), переданных в pcr. Дополнительные сведения см. в примерах ниже.

Выходные аргументы:

• x — вычисленное приближение к решению A * x = b.
• flag сообщает о сходимости. flag = 0 означает, что решение сошлось и выполнен критерий точности, заданный tol. flag = 1 означает, что достигнуто ограничение maxit на число итераций. flag = 3 сообщает о сбое pcr; подробности см. в [1].
• relres — отношение конечной невязки к её начальному значению, измеренное в евклидовой норме.
• iter — фактическое число выполненных итераций.
• resvec описывает историю сходимости метода: resvec (i) содержит евклидовы нормы невязки после (i-1)-й итерации, i = 1,2, …, iter+1.

Рассмотрим простую задачу с диагональной матрицей (используем разреженность A).

n = 10;
A = sparse (diag (1:n));
b = rand (N, 1);

ПРИМЕР 1: Простейшее использование pcr

x = pcr (A, b)

ПРИМЕР 2: pcr с функцией, вычисляющей A * x.

function y = apply_a (x)
  y = [1:10]' .* x;
endfunction

x = pcr ("apply_a", b)

ПРИМЕР 3: Итерационный метод с предобуславливанием и полной диагностикой. Предобуславливатель (довольно необычный, поскольку исходная матрица A тоже тривиальна) задан в виде функции.

function y = apply_m (x)
  k = floor (length (x) - 2);
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], "apply_m")
semilogy ([1:iter+1], resvec);

ПРИМЕР 4: Наконец, предобуславливатель, зависящий от параметра k.

function y = apply_m (x, varargin)
  k = varargin{1};
  y = x;
  y(1:k) = x(1:k) ./ [1:k]';
endfunction

[x, flag, relres, iter, resvec] = ...
                   pcr (A, b, [], [], "apply_m"', [], 3)

Литература:

W. Hackbusch, Iterative Solution of Large Sparse Systems of Equations, раздел 9.5.4; Springer, 1994

См. также: sparse, pcg.

Сходимость итерационного решателя к решению можно ускорить с помощью матрицы предобуславливания M. В этом случае вместо этого решается линейное уравнение M^-1 * x = M^-1 * A \ b. Типичные матрицы предобуславливания представляют собой частичные факторизации исходной матрицы.

: L = ichol (A)

: L = ichol (A, opts)

Вычисляет неполную факторизацию Холецкого разреженной квадратной матрицы A.

По умолчанию ichol использует только нижний треугольник A и формирует нижнюю треугольную матрицу L, такую что L*L' приближает A.

Факторизация, вычисленная этой функцией, может быть полезна в качестве предобуславливателя для системы линейных уравнений, решаемой итерационными методами, такими как PCG (метод сопряжённых градиентов с предобуславливанием).

Факторизацию можно изменить, передав параметры в структуре opts. Имя параметра является полем структуры, а его настройка — значением этого поля. Имена и спецификаторы чувствительны к регистру.

type

Тип факторизации.

"nofill" (по умолчанию)

Неполная факторизация Холецкого без заполнения (IC(0)).

"ict"

Неполная факторизация Холецкого с отбрасыванием по порогу (ICT).

diagcomp

Неотрицательный скаляр alpha для неполной факторизации Холецкого матрицы A + alpha * diag (diag (A)) вместо A. Это может быть полезно, если A не является положительно определённой. Значение по умолчанию — 0.

droptol

Неотрицательный скаляр, задающий порог отбрасывания при факторизации ICT. Значение по умолчанию — 0, что приводит к полной факторизации Холецкого.

Внедиагональные элементы L устанавливаются равными 0, если не выполняется условие

abs (L(i,j)) >= droptol * norm (A(j:end, j), 1).

michol

Модифицированная неполная факторизация Холецкого:

"off" (по умолчанию)

Суммы строк и столбцов не обязательно сохраняются.

"on"

Диагональ L изменяется так, чтобы суммы строк (и столбцов) сохранялись даже при отбрасывании элементов во время факторизации. Сохраняется соотношение: A * e = L * L' * e, где e — вектор из единиц.

shape

"lower" (по умолчанию)

Используется только нижний треугольник A и возвращается нижняя треугольная матрица L, такая что L*L' приближает A.

"upper"

Используется только верхний треугольник A и возвращается верхняя треугольная матрица U, такая что U'*U приближает A.

ПРИМЕРЫ

Следующая задача показывает, как выполнить полное и неполное разложение Холецкого для симметричной положительно определённой матрицы.

A = [ 0.37, -0.05,  -0.05,  -0.07;
     -0.05,  0.116,  0.0,   -0.05;
     -0.05,  0.0,    0.116, -0.05;
     -0.07, -0.05,  -0.05,   0.202];
A = sparse (A);
nnz (tril (A))
ans =  9
L = chol (A, "lower");
nnz (L)
ans =  10
norm (A - L * L', "fro") / norm (A, "fro")
ans =  1.1993e-16
opts.type = "nofill";
L = ichol (A, opts);
nnz (L)
ans =  9
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.019736

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

nx = 400; ny = 200;
hx = 1 / (nx + 1); hy = 1 / (ny + 1);
Dxx = spdiags ([ones(nx, 1), -2*ones(nx, 1), ones(nx, 1)],
               [-1 0 1 ], nx, nx) / (hx ^ 2);
Dyy = spdiags ([ones(ny, 1), -2*ones(ny, 1), ones(ny, 1)],
               [-1 0 1 ], ny, ny) / (hy ^ 2);
A = -kron (Dxx, speye (ny)) - kron (speye (nx), Dyy);
nnz (tril (A))
ans =  239400
opts.type = "nofill";
L = ichol (A, opts);
nnz (tril (A))
ans =  239400
norm (A - L * L', "fro") / norm (A, "fro")
ans =  0.062327

Литература по реализованным алгоритмам:

[1] Y. Saad. «Методы предобуславливания». Iterative Methods for Sparse Linear Systems, PWS Publishing Company, 1996.

[2] M. Jones, P. Plassmann: An Improved Incomplete Cholesky Factorization, 1992.

См. также: chol, ilu, pcg.

: LUA = ilu (A)

: LUA = ilu (A, opts)

: [L, U] = ilu (…)

: [L, U, P] = ilu (…)

Вычисляет неполную LU-факторизацию разреженной квадратной матрицы A.

ilu возвращает унитреугольную нижнюю треугольную матрицу L, верхнюю треугольную матрицу U и, при необходимости, матрицу перестановок P, такие что L*U приближает P*A.

Факторы, вычисленные этой функцией, могут быть полезны в качестве предобуславливателей для систем линейных уравнений, решаемых итерационными методами, такими как BICG (метод бисопряжённых градиентов) или GMRES (обобщённый метод минимальной невязки).

Факторизацию можно изменить, передав параметры в структуре opts. Имя параметра является полем структуры, а его настройка — значением этого поля. Имена и спецификаторы чувствительны к регистру.

type

Тип факторизации.

"nofill" (по умолчанию)

Факторизация ILU без заполнения (ILU(0)).

Дополнительные поддерживаемые параметры: milu.

"crout"

Вариант факторизации ILU по Крауту (ILUC).

Дополнительные поддерживаемые параметры: milu, droptol.

"ilutp"

Факторизация ILU с порогом и выбором главного элемента.

Дополнительные поддерживаемые параметры: milu, droptol, udiag, thresh.

droptol

Неотрицательный скаляр, задающий порог отбрасывания при факторизации. Значение по умолчанию — 0, что приводит к полной LU-факторизации.

Внедиагональные элементы U устанавливаются равными 0, если не выполняется условие

abs (U(i,j)) >= droptol * norm (A(:,j)).

Внедиагональные элементы L устанавливаются равными 0, если не выполняется условие

abs (L(i,j)) >= droptol * norm (A(:,j))/U(j,j).

milu

Модифицированная неполная LU-факторизация:

"row"

Модифицированная неполная LU-факторизация с сохранением сумм строк. Факторизация сохраняет суммы строк: A * e = L * U * e, где e — вектор из единиц.

"col"

Модифицированная неполная LU-факторизация с сохранением сумм столбцов. Факторизация сохраняет суммы столбцов: e' * A = e' * L * U.

"off" (по умолчанию)

Суммы строк и столбцов не обязательно сохраняются.

udiag

Если значение истинно, все нули на диагонали верхнего треугольного фактора заменяются локальным порогом отбрасывания droptol * norm (A(:,j))/U(j,j). Значение по умолчанию — false.

thresh

Порог выбора главного элемента при факторизации. Он может принимать значения от 0 (выбор диагонального элемента) до 1 (значение по умолчанию), при котором в качестве главного элемента выбирается элемент столбца с наибольшим модулем.

Если ilu вызывается только с одним выходным аргументом, возвращаемая матрица равна L + U - speye (size (A)), где L — унитреугольная нижняя треугольная матрица, а U — верхняя треугольная матрица.

При двух выходных аргументах ilu возвращает унитреугольную нижнюю треугольную матрицу L и верхнюю треугольную матрицу U. Для opts.type == "ilutp" один из факторов переставляется в соответствии со значением opts.milu. Если opts.milu == "row", U является верхним треугольным фактором с переставленными столбцами. В противном случае L является унитреугольным нижним треугольным фактором с переставленными строками.

Если заданы три именованных выходных аргумента и opts.milu != "row", возвращается P, такая что L и U являются неполными факторами P*A. Если opts.milu == "row", возвращается P, такая что L и U являются неполными факторами A*P.

ПРИМЕРЫ

A = gallery ("neumann", 1600) + speye (1600);
opts.type = "nofill";
nnz (A)
ans = 7840

nnz (lu (A))
ans = 126478

nnz (ilu (A, opts))
ans = 7840

Здесь показано, что A содержит 7,840 ненулевых элементов, полная LU-факторизация содержит 126,478 ненулевых элементов, а неполная LU-факторизация с нулевым уровнем заполнения содержит 7,840 ненулевых элементов — столько же, сколько A. Заимствовано из: https://www.mathworks.com/help/matlab/ref/ilu.html

A = gallery ("wathen", 10, 10);
b = sum (A, 2);
tol = 1e-8;
maxit = 50;
opts.type = "crout";
opts.droptol = 1e-4;
[L, U] = ilu (A, opts);
x = bicg (A, b, tol, maxit, L, U);
norm (A * x - b, inf)

В этом примере ILU используется как предобуславливатель для случайной FEM-матрицы с большим числом обусловленности. Без L и U метод BICG не сошёлся бы.

См. также: lu, ichol, bicg, gmres.