Определение системы линейных уравнений с помощью разреженной матрицы

При использовании distributed функция MATLAB автоматически запускает параллельный пул с использованием настроек кластера по умолчанию. В этом примере используется матрица Wathen из MATLAB gallery функция. Эта матрица является разреженной, симметричной и случайной матрицей с общей размерностью $$N\mathrm{=}{}^{\mathrm{}}\mathrm{3n2}+\mathrm{}$$4n + 1.

n = 400;
A = distributed(gallery('wathen',n,n));

Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).

N = 3*n^2+4*n+1

N = 481601

Теперь можно определить правый вектор $$b$$. В этом примере $$b$$ определяется как сумма строк $$A$$, что приводит к точному решению $$\mathrm{Ax}=$$ b $${\mathrm{вида}}_{}xexact\mathrm{}=[1,\mathrm{.}{.}^{.}$$, 1] T.

b = sum(A,2);

С тех пор sum действует в распределенном массиве, b также распределяется, и его данные хранятся в памяти работников вашего параллельного пула. Наконец, можно определить точное решение для сравнения с решениями, полученными итеративными методами.

xExact = ones(N,1,'distributed');

Решите систему линейных уравнений с помощью метода сопряженного градиента

Функция MATLAB pcg обеспечивает метод сопряженного градиента (CG), который итеративно генерирует ряд приближенных решений для $\mathit{x}$, улучшая раствор с каждой стадией.

[xCG_1,flagCG_1,relres_CG1,iterCG_1,resvecCG_1] = pcg(A,b);

Когда система решена, можно проверить ошибку между каждым элементом полученного результата. xCG_1 и ожидаемые значения xExact. Ошибка в вычисленном результате относительно высока.

errCG_1 = abs(xExact-xCG_1);

figure(1)
hold off
semilogy(errCG_1,'o');
title('System of Linear Equations with Sparse Matrix');
ylabel('Absolute Error');
xlabel('Element in x');

Итеративное вычисление заканчивается, когда ряд приближенных решений сходится к определенному допуску или после максимального числа шагов итерации. Как для распределенных, так и для внутриклиентных массивов $$\mathrm{pcg}$$ использует одни и те же настройки по умолчанию:

В качестве второго выходного аргумента pcg функция также возвращает флаг сходимости, который предоставляет дополнительную информацию о полученном результате, в том числе о том, сходилось ли вычисленное решение к требуемому допуску. Например, значение 0 указывает на правильное совпадение решения.

flagCG_1

flagCG_1 = 1

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

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

tolerance = 1e-12;
maxit = N;

tCG = tic;
    [xCG_2,flagCG_2,relresCG_2,iterCG_2,resvecCG_2] = pcg(A,b,tolerance,maxit);
tCG = toc(tCG);

flagCG_2

flagCG_2 = 0

С помощью пользовательских настроек решение сходится. Это решение имеет улучшенную абсолютную ошибку по сравнению с предыдущим решением.

errCG_2 = abs(xExact-xCG_2);
figure(2)
hold off
semilogy(errCG_1,'o');
hold on
semilogy(errCG_2,'d');
title('Comparison of Absolute Error');
ylabel('Absolute Error');
xlabel('Element in x');
legend('Default tolerance and iterations','Improved tolerance and iterations');
hold off

pcg метод также возвращает вектор остаточной нормы на каждом шаге итерации, norm(b-A*x)/norm(b). Относительная остаточная норма показывает отношение точности между последовательными шагами итерации. Эволюция остатков в течение итерационного процесса может помочь понять, почему решение не сходится без пользовательских настроек.

figure(3)
f=semilogy(resvecCG_2./resvecCG_2(1));
hold on
semilogy(f.Parent.XLim,[1e-6 1e-6],'--')
semilogy([20 20], f.Parent.YLim,'--')
semilogy(f.Parent.XLim,[1e-12 1e-12],'--')
title('Evolution of Relative Residual');
ylabel('Relative Residual');
xlabel('Iteration Step');
legend('Residuals of CG','Default Tolerance','Default Number of Steps','Custom Tolerance')
hold off

Очевидно, что количество шагов по умолчанию недостаточно для достижения хорошего решения для этой системы.

Решите свою систему линейных уравнений с помощью метода предварительно кондиционированного сопряженного градиента

Можно повысить эффективность решения системы с помощью метода предварительно кондиционированного сопряженного градиента (PCG). Во-первых, подготовьте систему линейных уравнений с помощью матрицы предварительного кондиционирования $\mathit{}$M. Затем решите систему предварительного кондиционирования с помощью метода CG. Метод PCG может принимать гораздо меньше итераций, чем метод CG.

Функция MATLAB pcg также используется для метода PCG. Можно ввести подходящую матрицу предварительного кондиционирования $\mathit{}$Mas в качестве дополнительного ввода.

Идеальной матрицей предварительного условия является матрица, обратная ${\mathit{}}^{\mathrm{M-1\; которой}}$является близкой аппроксимацией к обратной матрице коэффициентов, ${\mathit{}}^{\mathrm{A-1}}$, но её легче вычислить. В этом примере диагональ $$\mathrm{А\; используется}$$для предварительного кондиционирования системы линейных уравнений.

M = spdiags(spdiags(A,0),0,N,N);
tPCG = tic;
    [xPCG,flagPCG,relresPCG,iterPCG,resvecPCG]=pcg(A,b,tolerance,maxit,M);
tPCG = toc(tPCG);

figure(4)
hold off;
semilogy(resvecCG_2./resvecCG_2(1))
hold on;
semilogy(resvecPCG./resvecPCG(1))
title('Evolution of Relative Residual');
ylabel('Relative Residual');
xlabel('Iteration Step');
legend('Residuals of CG','Residuals of PCG with M \approx diag(A)')

На предыдущем рисунке показано, что метод PCG требует значительно меньшего количества шагов для схождения по сравнению с непринятой системой. Этот результат также отражается во времени выполнения.

fprintf([...
    '\nTime to solve system with CG:  %d s', ...
    '\nTime to solve system with PCG: %d s'],tCG,tPCG);

Time to solve system with CG:  1.244593e+01 s
Time to solve system with PCG: 7.657432e-01 s

Помимо решения этой примерной системы за меньшее количество шагов итерации, метод PCG также возвращает более точное решение.

errPCG = abs(xExact-xPCG);
figure(5)
hold off
semilogy(errCG_1,'o');
hold on
semilogy(errCG_2,'d');
semilogy(errPCG,'x');
title('Comparison of absolute error');
ylabel('Absolute error');
xlabel('Element in x');
legend('CG default','CG custom','PCG');

После завершения вычислений можно удалить параллельный пул. gcp функция возвращает текущий объект параллельного пула, чтобы можно было удалить текущий пул.

delete(gcp('nocreate'))

Матрица Ватена, используемая в этом примере, является хорошей демонстрацией того, как хороший специалист по предварительному кондиционированию может значительно повысить эффективность решения. У матрицы Wathen есть относительно маленькие недиагональные компоненты, таким образом$$\text{, choosingM=diag}$$(A) дает подходящий предварительный кондиционер. Для произвольной матрицы $\mathit{}$А поиск предварительного условия может быть не таким простым.

Пример того, как аппроксимировать дифференциальное уравнение линейной системой и решить его с помощью распределенного итерационного решателя с помощью мультигридного предварительного условия, см. в разделе Решение дифференциального уравнения с помощью мультигридного предварительного условия о распределенной дискретизации.