Решите квадратную линейную систему, используя gmres с настройками по умолчанию, а затем настроить допуск и количество итераций, используемых в процессе решения.

Создайте случайную разреженную матрицу A с плотностью 50% и ненулями на основной диагонали. Также создайте случайный вектор b для правой стороны $\mathrm{Ax}=\mathit{b}$.

rng default
A = sprandn(400,400,0.5) + 12*speye(400);
b = rand(400,1);

Решить $\mathrm{Ax}=\mathit{b}$ использование gmres. Выход отображения включает в себя значение относительной остаточной ошибки $\frac{\Vert \mathit{b}-\mathrm{Ax}\Vert }{\Vert \mathit{b}\Vert }$.

x = gmres(A,b);

gmres stopped at iteration 10 without converging to the desired tolerance 1e-06
because the maximum number of iterations was reached.
The iterate returned (number 10) has relative residual 0.35.

По умолчанию gmres использует 10 итераций и допуск 1e-6и алгоритм не может сходиться в этих 10 итерациях для этой матрицы. Поскольку невязка все еще велика, это хороший показатель того, что необходимо больше итераций (или матрица предкондиционера). Можно также использовать больший допуск, чтобы упростить сходимость алгоритма.

Снова решить систему с помощью допуска 1e-4 и 100 итераций.

tol = 1e-4;
maxit = 100;
x = gmres(A,b,[],tol,maxit);

gmres stopped at iteration 100 without converging to the desired tolerance 0.0001
because the maximum number of iterations was reached.
The iterate returned (number 100) has relative residual 0.0045.

Даже с более свободным допуском и большим количеством итераций остаточная ошибка недостаточно улучшается для сходимости. Когда итерационный алгоритм застопорится таким образом, это хорошее указание на то, что необходима матрица предкондиционера. Однако gmres также имеет вход, который управляет количеством внутренних итераций. Путем определения значения для внутренних итераций, gmres делает больше работы на каждую внешнюю итерацию.

Снова решить систему с помощью restart значение 100 и a maxit значение 20. Вместо того, чтобы делать 100 итераций один раз, gmres выполняет 100 итераций между перезапусками и повторяет это 20 раз.

restart = 100;
maxit = 20;
x = gmres(A,b,restart,tol,maxit);

gmres(100) converged at outer iteration 2 (inner iteration 75) to a solution with relative residual 9.3e-05.

В этом случае задается большое значение перезапуска для gmres позволяет ему сходиться к решению в пределах допустимого количества итераций. Однако большие значения перезапуска могут потреблять много памяти при A также является большим.

Исследуйте эффект использования матрицы предкондиционера с не перезапущенными gmres для решения линейной системы.

Загрузка west0479, действительная 479 на 479 несимметричная разреженная матрица.

load west0479
A = west0479;

Определите b чтобы истинное решение $\mathrm{Ax}=\mathit{b}$ является вектором всех таковых.

b = sum(A,2);

Установите допуск и максимальное количество итераций.

tol = 1e-12;
maxit = 20;

Использование gmres найти решение с требуемым допуском и количеством итераций. Укажите пять выходы, чтобы вернуть информацию о процессе решения:

[x,fl0,rr0,it0,rv0] = gmres(A,b,[],tol,maxit);
fl0

fl0 = 1

rr0

rr0 = 0.7603

it0

it0 = 1×2

     1    20

fl0 равен 1, потому что gmres не сходится к запрошенному допуску 1e-12 в пределах запрошенных 20 итераций. Лучшее приблизительное решение, которое gmres возвраты - последний (как обозначено it0(2) = 20). MATLAB хранит историю остатков в rv0.

Для помощи с медленной сходимостью можно задать матрицу предкондиционера. Начиная с A несимметричен, используйте ilu чтобы сгенерировать предварительный кондиционер $\mathit{M}=\mathit{L}\mathit{U}$. Задайте допуск удаления, чтобы игнорировать недиагональные значения со значениями меньше 1e-6. Решение предварительно обусловленной системы ${\mathit{M}}^{-1}\mathit{A}\mathit{x}={\mathit{M}}^{-1}\mathit{b}$ путем определения L и U в качестве входов в gmres. Обратите внимание, что когда вы задаете предварительное условие, gmres вычисляет норму невязки предварительно обусловленной системы для выходов rr1 и rv1.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
[x1,fl1,rr1,it1,rv1] = gmres(A,b,[],tol,maxit,L,U);
fl1

fl1 = 0

rr1

rr1 = 7.2303e-14

it1

it1 = 1×2

     1     6

Использование ilu предварительное кондиционер создает относительную невязку меньше, чем предписанный допуск 1e-12 на шестой итерации. Первый остаточный rv1(1) является norm(U\(L\b)), где M = L*U. Последняя невязка rv1(end) является norm(U\(L\(b-A*x1))).

Можно следить за прогрессом gmres путем построения графика относительных невязок при каждой итерации. Постройте график истории остатков каждого решения с линией для заданного допуска.

semilogy(0:length(rv0)-1,rv0/norm(b),'-o')
hold on
semilogy(0:length(rv1)-1,rv1/norm(U\(L\b)),'-o')
yline(tol,'r--');
legend('No preconditioner','ILU preconditioner','Tolerance','Location','East')
xlabel('Iteration number')
ylabel('Relative residual')

Использование предкондиционера с перезапущенными gmres.

Загрузка west0479, действительная 479 на 479 несимметричная разреженная матрица.

load west0479
A = west0479;

Определите b чтобы истинное решение $\mathrm{Ax}=\mathit{b}$ является вектором всех таковых.

b = sum(A,2);

Создайте неполный предварительный кондиционер LU с допуском на падение 1e-6.

[L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));

Преимущество использования перезапущенных gmres ограничивает объем памяти, требуемый для выполнения метода. Без перезапуска, gmres требует maxit векторы памяти для хранения базиса подпространства Крылова. Кроме того, gmres должна быть ортогональна всем предыдущим векторам на каждом шаге. Перезапуск ограничивает количество используемой рабочей области и объем работы, проделанной за внешнюю итерацию.

Выполните gmres(3), gmres(4), и gmres(5) использование неполных факторов LU в качестве предварительных кондиционеров. Используйте допуск 1e-12 и максимум 20 внешних итераций.

tol = 1e-12;
maxit = 20;
[x3,fl3,rr3,it3,rv3] = gmres(A,b,3,tol,maxit,L,U);
[x4,fl4,rr4,it4,rv4] = gmres(A,b,4,tol,maxit,L,U);
[x5,fl5,rr5,it5,rv5] = gmres(A,b,5,tol,maxit,L,U);
fl3

fl3 = 0

fl4

fl4 = 0

fl5

fl5 = 0

fl3, fl4, и fl5 все 0, потому что в каждом случае перезапускаются gmres приводит относительную невязку к меньшему, чем предписанный допуск 1e-12.

Следующий график показывает историю сходимости каждого перезапущенного gmres способ. gmres(3) сходится во внешней итерации 5, внутренней итерации 3 (it3 = [5, 3]), которая будет такой же, как внешняя итерация 6, внутренняя итерация 0, следовательно, маркировка 6 на конечной отметке деления.

semilogy(1:1/3:6,rv3/norm(U\(L\b)),'-o');
h1 = gca;
h1.XTick = (1:6);
title('gmres(N) for N = 3, 4, 5')
xlabel('Outer iteration number');
ylabel('Relative residual');
hold on
semilogy(1:1/4:3,rv4/norm(U\(L\b)),'-o');
semilogy(1:1/5:2.8,rv5/norm(U\(L\b)),'-o');
yline(tol,'r--');
hold off
legend('gmres(3)','gmres(4)','gmres(5)','Tolerance')
grid on

В целом, чем больше количество внутренних итераций, тем больше работы gmres делает по внешней итерации и тем быстрее она может сходиться.

Исследуйте эффект подачи gmres с начальным угадыванием решения.

Создайте трехугольную разреженную матрицу. Используйте сумму каждой строки в качестве вектора для правой оси графика $\mathrm{Ax}=\mathit{b}$ так, чтобы ожидаемое решение для $\mathit{x}$ является вектором таковых.

n = 900;
e = ones(n,1);
A = spdiags([e 2*e e],-1:1,n,n);
b = sum(A,2);

Использование gmres решить $\mathrm{Ax}=\mathit{b}$ дважды: один раз с начальным предположением по умолчанию и один раз с хорошим начальным предположением решения. Используйте 200 итераций и допуск по умолчанию для обоих решений. Задайте начальное предположение во втором решении как вектор со всеми элементами, равными 0.99.

maxit = 200;
x1 = gmres(A,b,[],[],maxit);

gmres converged at iteration 27 to a solution with relative residual 9.5e-07.

x0 = 0.99*e;
x2 = gmres(A,b,[],[],maxit,[],[],x0);

gmres converged at iteration 7 to a solution with relative residual 6.7e-07.

В этом случае подача начального предположения включает gmres чтобы сходиться быстрее.

x0 = zeros(size(A,2),1);
tol = 1e-8;
maxit = 100;
for k = 1:4
    [x,flag,relres] = gmres(A,b,[],tol,maxit,[],[],x0);
    X(:,k) = x;
    R(k) = relres;
    x0 = x;
end

Решить линейную систему путем предоставления gmres с указателем на функцию, который вычисляет A*x вместо матрицы коэффициентов A.

Одна из тестовых матриц Уилкинсона, сгенерированных gallery является тридиагональной матрицей 21 на 21. Предварительный просмотр матрицы.

A = gallery('wilk',21)

A = 21×21

    10     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     1     9     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     1     8     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     1     7     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     1     6     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     1     5     1     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     1     4     1     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     1     3     1     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     1     2     1     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     1     1     1     0     0     0     0     0     0     0     0     0     0
      ⋮

Матрица Уилкинсона имеет специальную структуру, поэтому можно представлять операцию A*x с указателем на функцию. Когда A умножает вектор, большинство элементов получившегося вектора являются нулями. Ненулевые элементы в результате соответствуют ненулевым триидиагональным элементам A. Кроме того, только главная диагональ имеет не ненули, которые не равны 1.

Выражение $\mathrm{Ax}$ становится:

$\mathrm{Ax}=\left[\begin{array}{cccccccccc}10& 1& 0& \cdots & & \cdots & & \cdots & 0& 0\\ 1& 9& 1& 0& & & & & & 0\\ 0& 1& 8& 1& 0& & & & & ⋮\\ ⋮& 0& 1& 7& 1& 0& & & & \\ & & 0& 1& 6& 1& 0& & & ⋮\\ ⋮& & & 0& 1& 5& 1& 0& & \\ & & & & 0& 1& 4& 1& 0& ⋮\\ ⋮& & & & & 0& 1& 3& \ddots & 0\\ 0& & & & & & 0& \ddots & \ddots & 1\\ 0& 0& \cdots & & \cdots & & \cdots & 0& 1& 10\end{array}\right]\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\\ {\mathit{x}}_3\\ {\mathit{x}}_4\\ {\mathit{x}}_5\\ ⋮\\ \\ ⋮\\ \\ {\mathit{x}}_{21}\end{array}\right]=\left[\begin{array}{c}10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}\end{array}\right]$.

Получившийся вектор может быть записан как сумма трех векторов:

$\mathrm{Ax}=\left[\begin{array}{c}0+10{\mathit{x}}_1+{\mathit{x}}_2\\ {\mathit{x}}_1+9{\mathit{x}}_2+{\mathit{x}}_3\\ {\mathit{x}}_2+8{\mathit{x}}_3+{\mathit{x}}_4\\ ⋮\\ {\mathit{x}}_{19}+9{\mathit{x}}_{20}+{\mathit{x}}_{21}\\ {\mathit{x}}_{20}+10{\mathit{x}}_{21}+0\end{array}\right]$=$\left[\begin{array}{c}0\\ {\mathit{x}}_1\\ ⋮\\ {\mathit{x}}_{20}\end{array}\right]+\left[\begin{array}{c}10{\mathit{x}}_1\\ 9{\mathit{x}}_2\\ ⋮\\ 10{\mathit{x}}_{21}\end{array}\right]+\left[\begin{array}{c}{\mathit{x}}_2\\ ⋮\\ {\mathit{x}}_{21}\\ 0\end{array}\right]$.

В MATLAB ® запишите функцию, которая создает эти векторы и добавляет их вместе, таким образом придавая значение A*x:

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end

(Эта функция сохранена как локальная функция в конце примера.)

Теперь решите линейную систему $\mathrm{Ax}=\mathit{b}$ путем предоставления gmres с указателем на функцию, который вычисляет A*x. Используйте допуск 1e-12, 15 внешних итераций и 10 внутренних итераций перед перезапуском.

b = ones(21,1);
tol = 1e-12;  
maxit = 15;
restart = 10;
x1 = gmres(@afun,b,restart,tol,maxit)

gmres(10) converged at outer iteration 5 (inner iteration 10) to a solution with relative residual 5.3e-13.

x1 = 21×1

    0.0910
    0.0899
    0.0999
    0.1109
    0.1241
    0.1443
    0.1544
    0.2383
    0.1309
    0.5000
      ⋮

Проверяйте это afun(x1) создает вектор таковых.

afun(x1)

ans = 21×1

    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
    1.0000
      ⋮

function y = afun(x)
y = [0; x(1:20)] + ...
    [(10:-1:0)'; (1:10)'].*x + ...
    [x(2:21); 0];
end