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

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

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

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

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

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

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

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

Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to parallel pool with 6 workers.

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

N = 481601

您现在可以定义右侧向量 $b$ 。在这个例子中， $b$ 被定义为 $A$ 的行和，这导致了 $A x = b$ 形式为 $x_{exact} = [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');

当一系列近似解收敛到特定容差或达到最大迭代步数时，迭代计算结束。对于分布式和客户端数组， $p c g$ 使用相同的默认设置：

默认最大公差为 $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:  9.616080e+00 s
Time to solve system with PCG: 1.399194e+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');