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 должна быть эрмитовой и положительно определённой (ЭПО). Если 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;
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, Итерационные методы для линейных и нелинейных уравнений, 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.

: 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"

Версия 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.