本页采用了机器翻译。点击此处可查看英文原文。

使用分布式数组通过迭代方法求解线性方程组

打开实时脚本

对于大规模数学计算，迭代方法比直接方法更有效。此示例显示如何使用分布式数组和迭代方法并行求解形式为 Ax=b 的线性方程组。

此示例延续了 使用分布式数组通过直接方法求解线性方程组 中涵盖的主题。mldivide 中实现的直接求解器方法可用于并行求解分布式线性方程组，但对于某些大型稀疏系统可能效率不高。迭代方法根据初始猜测生成一系列解决方案，经过几步后收敛到最终结果。与直接计算解决方案相比，这些步骤的计算量更小。

分布式数组将数据从客户端工作区分发到本地计算机或集群中的并行池。每个工作单元将数组的一部分存储在其内存中，但也可以与其他工作单元通信以访问数组的所有段。分布式数组可以包含不同类型的数据，包括完整矩阵和稀疏矩阵。

此示例使用 pcg 函数演示如何使用共轭梯度法和预条件共轭梯度法求解大型线性方程组。迭代方法可用于密集矩阵和稀疏矩阵，但对于稀疏矩阵系统最有效。

使用稀疏矩阵定义线性方程组

当您使用distributed函数时，软件会使用您的默认集群设置自动启动并行池。此示例使用了 gallery函数中的 Wathen 矩阵。该矩阵是一个稀疏、对称、随机矩阵，总体维度为 N=3n2+4n+1。

n = 400;
A = distributed(gallery("wathen",n,n));
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");

使用共轭梯度法求解线性方程组

pcg MATLAB 函数提供了共轭梯度 (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，若小于 20，则为系数矩阵 A 的阶。

作为第二个输出参量，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 方法少得多。

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:  8.819875e+00 s 
Time to solve system with PCG: 1.307804e+00 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"))
Parallel pool using the 'Processes' profile is shutting down.

本示例中使用的 Wathen 矩阵很好地证明了良好的预处理器如何能够显著提高解决方案的效率。Wathen 矩阵具有相对较小的非对角线分量，因此选择 M=diag(A) 可以提供合适的预条件子。对于任意矩阵 A，找到预条件子可能不是那么简单。

有关如何通过线性系统近似微分方程并使用具有多重网格预处理器的分布式迭代求解器求解的示例，请参阅 使用分布式离散化的多重网格预处理器求解微分方程。

另请参阅

| |

主题