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)

Ссылка:

В. Хаукбуш, Итеративное решение больших разреженных систем уравнений, раздел 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"

Версия факторизации 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

Если True, любые нули на диагонали верхней треугольной матрицы заменяются локальным порогом отбрасывания 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 имеет 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 используется как предобусловитель для случайной матрицы FEM, которая имеет большое число обусловленности. Без L и U BICG не сойдётся.

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