利用 GPU 加速求解线性方程

本示例演示了如何利用 GPU 加快线性方程组的求解速度。

例如，在此示例中，您将使用 CPU 和 GPU 求解一组线性方程组，并比较 CPU 和 GPU 的执行速度。在 GPU 上，您还可以同时使用全数组和稀疏数组的 gpuArray 输入来求解该系统。使用稀疏数组来存储包含大量零的矩阵，可以减少内存要求。

线性方程组

技术计算中一个重要问题是如何高效地求解线性方程组。在矩阵表示法中，可以将一组线性方程组写成如下形式：

$Ax = b$ .

尽管这并非标准的数学符号，但 MATLAB ® 采用了标量情况下常见的除法术语来描述一般联立方程组的解。反斜杠符号 \ 对应于 MATLAB 函数 mldivide。该运算符用于表示通过 mldivide 得到的矩阵方程 $Ax = b$ 的解：

$x = A \ b$ .

大多数线性方程组包含一个方阵系数矩阵 A 和一个右边列向量 b。

定义问题

该图描绘了一个平面桁架，该桁架由 13 根构件（带编号的线段）连接 8 个等间距的节点（带编号的圆圈）[1]。分别在接头 2、5 和 6 处施加标明的 10、15 和 20 吨载荷，目的是确定桁架各构件所受的内力。

Plane truss with 13 members connected by 8 joints.

为了使桁架处于静力平衡状态，任何节点处都不应存在水平或垂直方向的合力。因此，可以通过将每个节点的左右水平力相等来确定构件内力，同样地，也可以通过将每个节点的上下垂直力相等来确定构件内力。对于这八个关节，这将产生 16 个方程，而待求的未知量只有 13 个。为了使桁架成为静定结构，即存在唯一解，假设接头 1 在水平和垂直方向上均被刚性固定，且接头 8 在垂直方向上被固定。通过将构件内力分解为水平组件和垂直组件，并定义 $α = \frac{1}{\sqrt{2}}$ ，可以得到关于构件内力 $f_{i}$ 的如下方程组：

关节	力量
2	$f_{2} = f_{6}$ $f_{3} = 10$
3	$α f_{1} = f_{4} + α f_{5}$ $α f_{1} + f_{3} + α f_{5} = 0$
4	$f_{4} = f_{8}$ $f_{7} = 0$
5	$α f_{5} + f_{6} = α f_{9} + f_{10}$ $α f_{5} + f_{7} + α f_{9} = 15$
6	$f_{10} = f_{13}$ $f_{11} = 20$
7	$f_{8} + α f_{9} = α f_{12}$ $α f_{9} + f_{11} + α f_{12} = 0$
8	$f_{13} + α f_{12} = 0$

为了使您能够在 MATLAB 中求解 x，其中 x 是由力 $f_{2}$ 至 $f_{12}$ 组成的列向量，首先利用受力方程定义矩阵 A 和载荷列向量 b。A 和 b 这两行分别对应上表中的受力方程。

a = 1/sqrt(2);

A = [0 1 0 0 0 -1 0 0 0 0 0 0 0; ...
    0 0 1 0 0 0 0 0 0 0 0 0 0; ...
    a 0 0 -1 -a 0 0 0 0 0 0 0 0; ...
    a 0 1 0 a 0 0 0 0 0 0 0 0; ...
    0 0 0 1 0 0 0 -1 0 0 0 0 0; ...
    0 0 0 0 0 0 1 0 0 0 0 0 0; ...
    0 0 0 0 a 1 0 0 -a -1 0 0 0; ...
    0 0 0 0 a 0 1 0 a 0 0 0 0; ...
    0 0 0 0 0 0 0 0 0 1 0 0 -1; ...
    0 0 0 0 0 0 0 0 0 0 1 0 0; ...
    0 0 0 0 0 0 0 1 a 0 0 -a 0; ...
    0 0 0 0 0 0 0 0 a 0 1 a 0; ...
    0 0 0 0 0 0 0 0 0 0 0 a 1];

b = [0 10 0 0 0 0 0 15 0 20 0 0 0]';

解方程

使用函数 mldivide（运算符 \）求解 x。为了尽可能缩短计算时间，MATLAB 会根据矩阵 A 的属性采用不同的算法。

x = A\b

在使用 GPU 计算 x 之前，请先确认您使用的 GPU 是否受支持。

gpu = gpuDevice;
disp(gpu.Name + " GPU selected.")

NVIDIA RTX A5000 GPU selected.

将矩阵 A 转换为 gpuArray。这会将数据移动到 GPU 内存中。如果输入数据是 gpuArray，那么许多 MATLAB 线性代数函数（包括 mldivide）都会在 GPU 上运行。如果某个函数支持 gpuArray 输入，则会在该函数的参考页面“扩展功能”部分中予以说明。

A = gpuArray(A);

在 GPU 上求解 x。输出也是一个 gpuArray。

x = A\b;

验证计算得到的解能否准确重构向量 b。

A*x

使用 plotTrussForces 函数将计算出的力叠加到桁架构件上。此函数作为支持文件包含在本示例中。

plotTrussForces(x)

分别使用 timeit 和 gputimeit 函数比较 CPU 和 GPU 的执行时间。对于使用 GPU 的函数，gputimeit 比 timeit 更可取，因为它能确保在记录时间之前 GPU 上的所有操作都已完成，并能补偿由此产生的开销。不过，对于不使用 GPU 的运算，timeit 能提供更高的精度。

A = gather(A);
tCPU = timeit(@() A\b,1)

tCPU = 
3.2756e-06

A = gpuArray(A);
tGPU = gputimeit(@() A\b,1)

tGPU = 
0.0013

disp("Speedup of CPU over GPU for a small problem: " + round(tGPU/tCPU,1) + "x")

Speedup of CPU over GPU for a small problem: 387.6x

扩大规模

对于一个包含八个关节且具有 13×13 矩阵 A 的系统，CPU 求解该系统的速度明显快于 GPU。这是意料之中的，因为 GPU 拥有众多核，因此只有在处理大量数据时才能发挥其优势。

为了展示 GPU 计算能够发挥效力的规模，请定义一个更大的问题。使用 trussMatrix 和 randomLoad 函数生成矩阵 A 和列向量 b，以表示作用于一个具有 1000 个节点的平面桁架上的内力和外载荷。这些函数作为辅助文件附在示例中。

A = trussMatrix(1000);
b = randomLoad(1000,30);

比较 CPU 和 GPU 的执行时间。

tCPU = timeit(@() A\b,1)

tCPU = 
0.0638

A = gpuArray(A);
x = A\b;
tGPU = gputimeit(@() A\b,1)

tGPU = 
0.0246

disp("Speedup of GPU over CPU for a large problem: " + round(tCPU/tGPU,1) + "x")

Speedup of GPU over CPU for a large problem: 2.6x

验证计算得到的解能否准确重构向量 b。由于向量较大，因此仅比较前 10 个元素。

bCheck = A*x;
[b(1:10) bCheck(1:10)]

ans =

         0         0
   19.0000   19.0000
         0         0
         0         0
         0   -0.0000
         0         0
         0    0.0000
   18.0000   18.0000
         0         0
   17.0000   17.0000

利用矩阵的稀疏性

随着关节数量的增加，矩阵 A 的稀疏性也会增加；也就是说，元素为零的比例会增加。计算矩阵 A 的密度。只有大约 0.01% 的元素值不为零。

density = nnz(A)/numel(A)

density = 
0.0013

绘制矩阵 A 的稀疏性分布图。大多数非零元素都位于主对角线附近。

figure
spy(A)

MATLAB 中的稀疏数组可高效存储零值占比很高的矩阵。虽然全（或稠密）矩阵会将每个元素都存储在内存中（无论其值为何），但稀疏矩阵仅存储非零元素及其位置。因此，使用稀疏矩阵可以显著减少数据存储所需的内存量。有关详细信息，请参阅稀疏矩阵的计算优点和在 GPU 上使用稀疏数组。

将矩阵 A 转换为稀疏矩阵，并比较以全矩阵形式和稀疏矩阵形式存储时该矩阵的大小。所需内存的减少程度取决于矩阵的大小和结构。

sparseA = sparse(A);
whos A sparseA

  Name            Size                 Bytes  Class       Attributes

  A            1997x1997            31904072  gpuArray              
  sparseA      1997x1997               67872  gpuArray    sparse

使用稀疏 gpuArray 输入求解该方程组所需的时间。在 MATLAB 中，许多能够求解线性方程组的函数都支持稀疏数组，包括 mldivide。

tGPUSparse = gputimeit(@() sparseA\b,1)

tGPUSparse = 
0.0278

参考资料

[1] 莫勒，克利夫·B. 《使用 MATLAB 进行数值计算》。工业与应用数学学会，2004 年。https://doi.org/10.1137/1.9780898717952。

另请参阅

mldivide | gpuArray | sparse | spy