Используйте распределенные массивы, чтобы решить системы линейных уравнений с итерационными методами

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

Этот пример продолжает темы, затронутые в Использовании Распределенные Массивы, чтобы Решить Системы Линейных уравнений с Прямыми методами. Прямые методы решателя реализованы в mldivide может использоваться, чтобы решить распределенные системы линейных уравнений параллельно, но может не быть эффективным для определенных больших и разреженных систем. Итерационные методы генерируют серию решений от исходного предположения, сходясь к конечному результату после нескольких шагов. Эти шаги могут быть менее в вычислительном отношении интенсивными, чем вычисление решения непосредственно.

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

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

Задайте свою Систему Линейных уравнений с помощью Разреженной матрицы

Когда вы используете distributed функция, MATLAB автоматически запускает параллельный пул с помощью кластерных настроек по умолчанию. Этот пример использует матрицу Wathen от gallery MATLAB функция. Эта матрица является разреженной, симметричной, и случайной матрицей с габаритным размером N=3n2+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 из формы xexact=[1,...,1]T.

b = sum(A,2);

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

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

Решите свою Систему Линейных уравнений с Методом сопряженных градиентов

Функция MATLAB pcg предоставляет метод метода сопряженных градиентов (CG), который итеративно генерирует серию приближенных решений для 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 использование те же настройки по умолчанию:

  • Максимальный допуск по умолчанию 10-6.

  • Максимальное количество по умолчанию шагов итерации равняется 20 или порядку матрицы коэффициентов A если меньше чем 20.

Как второй выходной аргумент, 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). Во-первых, предусловие ваша система линейных уравнений с помощью матрицы перед формирователем M. Затем решите свою предобусловленную систему с помощью метода CG. Метод PCG может взять намного меньше итераций, чем метод CG.

Функция MATLAB pcg также используется в методе PCG. Можно предоставить подходящую матрицу перед формирователем Mкак дополнительный вход.

Идеальная матрица перед формирователем является матрицей чья инверсия M-1 близкое приближение к инверсии матрицы коэффициентов, 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');

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

delete(gcp('nocreate'))

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

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

Смотрите также

| |

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте