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

Когда вы используете distributed MATLAB автоматически запускает параллельный пул, используя настройки кластера по умолчанию. Этот пример использует матрицу Ватена из MATLAB gallery функция. Эта матрица является разреженной, симметричной и случайной матрицей с общей размерностью $$N=3n^2+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$$, что приводит к точному решению $$Ax=b$$ формы $$x_{\text{exact}}=[1,...,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');

Итеративный расчет заканчивается, когда серия приблизительных решений сходится к определенному допуску или после максимального количества шагов итерации. Как для распределенных, так и для клиентских массивов, $$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

The 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{M}$как дополнительный вход.

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

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');

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

delete(gcp('nocreate'))

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

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