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 может быть передана как матрица, дескриптор функции или встроенная функция Afun таким образом, что Afun(x) = A * x. Дополнительные параметры к Afun могут быть переданы после x0.

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

• b — вектор правой части.
• tol — требуемая относительная толерантность для погрешности остатка, b - A * x. Итерация останавливается, если norm (b - A * x) ≤ tol * norm (b). Если tol опущена или пуста, то используется толерантность 1e-6.
• maxit — максимальное допустимое количество итераций; если maxit опущена или пуста, то используется значение 20.
• m — матрица предварительного обуславливания ЭПО. Для любого разложения m = p1 * p2 такого, что inv (p1) * A * inv (p2) является ЭПО, метод сопряжённых градиентов формально применяется к линейной системе 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: алгоритм обнаруживает, что входная (предварительно обусловленная) матрица не является ЭПО.
• 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;
Afun = @(x) A * x;
Mfun = @(x) M \ x;
M1fun = @(x) M1 \ x;
M2fun = @(x) M2 \ x;

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

x = pcg (A, b)

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

x = pcg (Afun, b)

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

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

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

x = pcg (Afun, b, 1e-6, 100, Mfun)

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

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

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

x = pcg (Afun, b, 1e-6, 100, M1fun, M2fun)

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

function y = Ap (A, x, p) # compute A^p * x
     y = x;
     for i = 1:p
       y = A * y;
     endfor
  endfunction
Apfun = @(x, p) Ap (A, x, p);
x = pcg (Apfun, 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, Итерационные методы для линейных и нелинейных уравнений, SIAM, 1995. (основной алгоритм PCG)
2. Y. Saad, Итерационные методы для разреженных линейных систем, 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)

Литература:

[1] В. Хакбуш, Итеративное решение больших разреженных систем уравнений, раздел 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] Й. Саад. "Методы предварительного обуславливания." Итерационные методы для разреженных систем линейных уравнений, PWS Publishing Company, 1996.

[2] М. Джонс, П. Плассманн: Улучшенная неполная факторизация Холецкого, 1992.

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

ilu (A)

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"

Версия 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). По умолчанию значение ложно.

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 имеет 7840 ненулевых элементов, полное LU-разложение — 126478 ненулевых элементов, а неполное LU-разложение с 0 уровнем заполнения — 7840 ненулевых элементов, то есть столько же, сколько у 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 как предобуславливатель для случайной матрицы КЭМ, имеющей большой число обусловленности. Без L и U BICG не сойдётся.

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