一种基于快速矩阵分解法的近海气候预测方法和装置

1.本发明涉及一种基于快速矩阵分解法的近海气候预测方法和装置，该方法运行在分布式系统中，依赖于mpi以及openmp任务模型，可用于近海地区地理统计应用中如气温、降水量、风速、气压、湿度等物理数据的快速预测。


背景技术：

2.气候以及环境预测应用是当今各国超级计算机上运行的主要物理应用之一，该类应用对于科学家进行地理空间统计相关研究不可或缺。近海区域气候对于海洋渔业资源具有重大影响。受限于采集设备投放地点，通常需要根据现有观测数据以及模型对其余分布不均匀地点进行预测来获取完整气候物理数据。在预测前，科学家通常会定期或不定期地在指定范围的地理地区上进行大量物理数据(物理数据可以是气温，降水量，湿度，风速，气压等)测量，然后对这些数据使用高斯空间随机场模型进行建模。通过对基于大型协方差矩阵的高斯对数似然函数进行模型参数估计，来最终确定数据模型。借助已经获得的观测数值，该模型可以用于该地区其他地理点位物理数据的预测。这类预测问题的主要步骤可以归纳为一个大型对称正定线性系统求解过程，对称正定线性系统求解过程又可以进一步分解为矩阵cholesky预分解以及基于正反向传播的求解。由于求解过程计算量过于庞大，预测程序的主要运行时间都会集中在该过程，因此优化预测程序的求解过程是十分必要的。
3.一种优化方式是使用更高性能的求解算法。共享内存架构上基于openmp任务依赖图的并行模型，已被充分证明非常适合具有非规则特性的分解算法。此模型允许程序在运行时动态生成任务依赖图，并在不违反任务之间依赖的前提下，动态调度所有可用线程到分解任务上进行执行。与传统并行模型相比，这种模型消除了算法中的显式同步屏障，最大限度地利用了动态调度策略来实现性能和资源利用率的提升。然而在分布式环境下设计基于任务依赖图的分解算法十分困难，因为mpi标准本身没有提供任何与任务模型相关的功能，同时也与openmp等方案提供的任务模型之间没有任何互操作性，这使得算法难以被移植。目前存在的一些分布式环境下任务模型解决方案，都是绕过mpi标准以及openmp标准，由不同团队自己实现的工具。对于其他开发者而言，需要从头学习每种工具的api，大大增加了分解算法并行化学习成本。
4.另一种优化方式是使用低秩矩阵分解，预测程序中的协方差系数矩阵通常在整体上是稠密矩阵，但在局部上具有稀疏特征，可以使用低秩逼近算法对分解过程中计算强度过大的部分计算核心进行低秩化处理，从而以降低部分精确度为代价，有效地提升算法效率。
5.因此针对预测程序中的大规模低秩矩阵分解以及求解过程，有必要使用目前最新标准下的mpi openmp任务模型对其重新设计，并尝试在新的混合架构下对低秩分解的运行效率，内存占用等方面进行分析和优化。


技术实现要素：

6.本发明的目的在于针对现有技术的不足，提出一种基于快速矩阵分解法的近海气候预测方法，本发明提出一种新型mpi x范式(即mpi openmp任务模型)下稳定且可以轻松复用的并行化策略，以此实现依赖于新范式下低秩cholesky分解的近海气候快速预测方法。
7.本发明解决其技术问题采用的技术方案步骤如下：
8.一种基于快速矩阵分解法的近海气候预测方法，是使用已有基于matern核函数高斯随机场模型以及收集到近海气候物理数据，对该区域的任意空间点位进行物理值预测；使用openmp来生成任务依赖图并按依赖顺序调度其中的分解任务；使用mpi来控制矩阵分块后，然后进行不同进程间的非阻塞通信。该方法采用以下步骤实现：
9.步骤(1)：各进程对数据进行划分，并计算任意位置之间的距离；
10.1-1获取待预测位置的坐标数据矩阵l1、观测位置的坐标数据矩阵l2以及对应的物理值观测向量z2；
11.所述物理值观测向量z2包括气温、降水量、风速、气压、湿度等；
12.所述坐标数据矩阵包括各坐标点的经度和纬度值；
13.1-2配置多个基于mpi通信的mpi进程，各进程结合改进的二维数据布局，根据进程数对步骤1-1数据分布进行划分，各进程会分到若干个矩阵块；
14.1-3根据l1和l2，结合公式(1)各进程分别计算出任意两个观测位置之间的距离矩阵d1以及所有观测位置和待预测位置之间的距离矩阵d2；
15.步骤(2)：利用现有基于matern核函数高斯随机场模型，各进程获取任意位置之间观测值的相关性；
16.步骤(3)：各进程结合两组协方差矩阵∑1和∑2，求解物理值预测向量z1；具体是：
17.3-1根据克里金法，获取物理值预测向量z1计算公式：
[0018][0019]
3-2取z1*∑
2-1
＝x，然后对上述公式(4)建立对称正定线性系统∑1*x＝z2，求解x；具体是：
[0020]
1)分布式低秩cholesky分解
[0021]
2)基于cholesky分解结果的正反向传播
[0022]
步骤(4)：主进程利用完整解向量，根据标准矩阵乘法计算物理值预测向量z1＝∑2*x，即为近海气候预测值。
[0023]
本发明的另一个目的是提供一种近海气候预测装置，包括：
[0024]
数据获取模块，用于获取待预测位置的坐标数据矩阵l1、观测位置的坐标数据矩阵l2以及对应的物理值观测向量z2；
[0025]
并行数据确定模块，用于根据改进的二维数据布局确定openmp的线程数和mpi的进程数；
[0026]
第一计算模块，用于各mpi进程计算出任意两个观测位置之间的距离矩阵d1以及所有观测位置和待预测位置之间的距离矩阵d2；
[0027]
第二计算模块，利用现有基于matern核函数高斯随机场模型，各进程获取任意位置之间观测值的相关性；
[0028]
第三计算模块，包括对称正定线性系统构建模块、分布式低秩cholesky分解模块，正反向传播模块；对称正定线性系统构建模块，用于基于物理值预测向量构建对称正定线性系统方程；分布式低秩cholesky分解模块，用于利用openmp多线程执行分解任务、压缩任务、三角求解任务、尾更新任务、通信任务，获得cholesky分解结果；正反向传播模块，用于利用openmp多线程执行三角求解任务、向量更新任务、通信任务，实现对称正定线性系统的求解；
[0029]
第四计算模块，用于利用对称正定线性系统的解根据标准矩阵乘法计算得到物理值预测向量。
[0030]
本发明具有的有益效果是：
[0031]
1、本发明提出一种全新的mpi x范式(即mpi openmp任务模型)，既可以运行在分布式环境中，也可以使用原生的openmp任务模型进行编程，多线程执行通信可以有效提升大规模分解算法运行效率。
[0032]
2、本发明提出一种分布式环境下全新的二维数据布局，平衡分布式分解算法中各进程上负载；该布局能分摊计算量较高的任务到不同进程上，从而降低分解算法的通信开销，并提升效率。
[0033]
3、本发明将新范式应用于面向近海气候预测方法中分布式低秩cholesky方法中，并进一步提出一种高效分解和压缩策略来进行计算和通信的重叠；该策略具有非常高的运行效率，并在内存占用方面具有优势。
[0034]
4、本发明将新范式应用于面向近海气候预测方法中基于正反向传播的求解过程，使该过程也可以在分布式环境中使用openmp任务模型。
附图说明
[0035]
图1为本发明方法流程图。
[0036]
图2(a)为本发明中改进的分布式二维数据布局，图2(b)为传统分布式二维数据布局。
[0037]
图3为线程在通信任务中切换时的工作原理图。
[0038]
图4为更新任务中低秩更新过程图。
[0039]
图5为分解任务中任务执行区域图。
具体实施方式
[0040]
下面结合附图，对本发明的具体实施方案作进一步详细描述。其中：
[0041]
一种基于快速矩阵分解法的近海气候预测方法，如图1具体是：
[0042]
步骤(1)：各进程对数据进行划分，并计算任意位置之间的距离；具体是：
[0043]
1-1获取待预测位置的坐标数据矩阵l1、观测位置的坐标数据矩阵l2以及对应的物理值观测向量z2；
[0044]
所述物理值观测向量z2包括气温、降水量、风速、气压、湿度等；
[0045]
所述坐标数据矩阵包括各坐标点的经度和纬度值；
[0046]
1-2配置多个基于mpi通信的mpi进程，各进程结合改进的二维数据布局，根据进程数对步骤1-1数据分布进行划分，各进程会分到若干个矩阵块；
[0047]
所述改进的二维数据布局具体是在一个二维循环布局中存在多个单次循环区域，其中每个单次循环区域的边长为vp，则单次循环区域内每列都包含vp个不同进程序号的矩阵块；将整个二维循环布局中对角线上的进程序号依次替换为它所在单次循环区域对应列上的序号；若是对角线上的进程序号所在的单次循环区域位于第i行，则替换的序号为所在列的第i行进程序号；
[0048]
如图2(a)-(b)所示，其二者区别在于对角线上的进程顺序。以图2中9个进程的布局为例，传统布局仅仅是进行粗略划分，对角线上数据块只有0、4、8号进程持有。改进型布局则是对对角线上数据块进一步地划分，如图2(a)虚线块所示。这种划分使得所有进程都拥有对角线数据块，分担了计算量较大的对角线矩阵块任务，另外此改动也在最大程度上保留了传统布局对于数据划分的平衡性。
[0049]
1-3根据l1和l2，结合公式(1)各进程分别计算出任意两个观测位置之间的距离矩阵d1以及所有观测位置和待预测位置之间的距离矩阵d2；
[0050]
其中距离矩阵计算方式如下所示：
[0051][0052]
其中e函数为标准欧氏距离计算函数，(x1,y1)、(x2,y2)为两个任意位置的坐标经纬度；
[0053]
步骤(2)：利用现有基于matern核函数高斯随机场模型，各进程获取任意位置之间观测值的相关性；具体是：
[0054]
所述基于matern核函数高斯随机场模型中matern核函数如下：
[0055][0056]
其中θ1代表方差，θ2代表空间范围，θ3代表平滑度，γ代表标准gamma函数，k代表根据平滑度θ3修正后的贝塞尔曲线函数；
[0057]
基于matern核函数的协方差矩阵计算方式如下所示：
[0058][0059]
其中σ矩阵中对角线位置元素为θ1，非对角线位置元素为任意位置之间matern核函数的计算结果；
[0060]
对步骤(1)获得的距离矩阵d1和d2，根据公式(2)-(3)分别得到两组协方差矩阵∑1和∑2；
[0061]
步骤(3)：各进程结合两组协方差矩阵∑1和∑2，求解物理值预测向量z1；具体是：
[0062]
3-1根据克里金法，获取物理值预测向量z1计算公式：
[0063][0064]
3-2取z1*∑
2-1
＝x，然后对上述公式(4)建立对称正定线性系统∑1*x＝z2，求解x；具体是：
[0065]
3)分布式低秩cholesky分解；本发明中低秩cholesky算法主要步骤的伪代码如下所示：
[0066][0067][0068]
s11：初始化mpi，并指定mpi启动参数为mpi_thread_multiple，使得每个进程中同时启用openmp多线程进行独立的非阻塞通信；上述多线程包括一个主线程和多个从线程；主线程负责创建所有任务，并生成任务依赖图；在主线程工作任务的同时，其余空闲从线程会自动根据任务依赖图到openmp任务队列中获取可执行任务进行绑定并执行；其中openmp任务队列中包括分解任务、压缩任务、三角求解任务、尾更新任务、通信任务；
[0069]
所述任务依赖图中的图节点为openmp任务队列中的各任务，两图节点的边为两任务的前后依赖关系；
[0070]
s12：某一进程选取某线程执行分解任务，具体在对角线块上采用标准型cholesky分解；图5为分解任务中任务执行区域图。
[0071]
s13：当前从线程执行完分解任务后，由所有进程的线程进行通信任务同步，接着释放通信任务的依赖；然后执行分解任务后的进程选取任意一个从线程执行通信任务，将cholesky分解结果发送给对应进程；同步任务由一个单独的依赖变量完成，在此任务中输
出该变量的依赖，同时在其他所有通信任务中接收该变量为任务依赖，这样的同步方式可以保证上述伪代码中第一个通信区域和第二个通信区域中通信任务不会发生区域间的重叠；在不存在全局依赖的情况下，该同步屏障可以避免全局依赖不正确导致的死锁；
[0072]
图3为线程在通信任务中切换时的工作原理图。在通信任务中，为非阻塞通信任务设置显式线程调度点；当线程执行发送或者接收操作后遇到调度点时，会暂停当前任务的执行，并切换到openmp任务队列中其他可执行的任务；等当前线程再次可调度时，回到调度点处继续执行原来的任务；调度前通过mpi进行通信情况的测试，线程会轮询通信情况并判断需要调度或者继续完成任务；该通信任务的设计方式解决了mpi openmp任务模型线程阻塞导致的效率问题；线程切换的工作原理如图3所示。以图3为例，t1时刻，该线程进入一个通信任务开始执行。到t2时刻，由于mpi_test的结果为假，通信任务未完成，所以进行切换其他任务上(可以是计算任务或者其他通信任务)。在t3时刻，非阻塞通信完成，该通信任务可以继续执行下去，但是其他任务中可能不存在同样的切换点，所以必须等待该线程在t4时刻完成那个任务并再次进入调度状态。在t4时刻，当线程再次调度并执行原通信任务时，发现可以跳出while循环，则继续执行至任务末尾，在t5时刻彻底结束该任务。
[0073]
本发明中的发送类型通信任务的伪代码如下(接收类型通信任务与该伪代码类似)：
[0074][0075]
s14：不采用现有方法中预压缩处理，而是在每次循环的步骤s12后，由各进程选取某从线程依次执行压缩任务和三角求解任务；所述压缩任务具体使用自适应交叉算法adaptive cross approximation对矩阵块进行低秩压缩，得到两个压缩结果矩阵x和y；所述三角求解任务，具体对压缩过后矩阵块进行低秩版本的三角求解，获得求解结果过程为：
[0076]
ar＝x*(y*ac-1
)＝x*y'
ꢀꢀꢀ
公式(5)
[0077]
其中ar为每次三角求解后的结果，该结果可由压缩结果x和y，以及在通信中接收到的分解结果逆矩阵ac-1
获得；该步骤使用的低秩三角求解算法，不会完全计算出最终结果，而是继续保留低秩结果表示形式，如结合公式(5)的后两项为y'；这一方法虽然会增加第二个通信区域的任务执行时间(相比全秩通信，需要以全秩尺寸传输两个低秩矩阵)，但
却是步骤5执行更高效低秩更新的必要前提；
[0078]
s15：当前从线程执行完三角求解任务后，由所有进程的线程进行通信任务同步，然后释放通信任务的依赖；最后执行过低秩三角求解的进程会将两个低秩矩阵和对应矩阵秩，发送给对应进程；该通信任务的实现方式同步骤s13；
[0079]
s16：各进程选取某从线程对所有尾矩阵块进行尾更新任务。图4为更新任务中低秩更新过程图。
[0080]
所述尾更新任务具体求解过程如下：
[0081]
s16-1：构建尾更新公式：
[0082]
br＝b-x1*((y'1*x2)*y'2)
ꢀꢀꢀꢀ
公式(6)
[0083]
其中，br为每次尾更新后的结果，可由进程所在矩阵块b，以及在步骤s15中，通过通信获得两个矩阵的低秩结果表示形式(x1，y'1)和(x2，y'2)；
[0084]
s16-2：首先计算公式(6)中间结果w1＝y'1*x2，再基于w1求解下一中间结果w2＝w1*y'2；由于低秩矩阵具有较小秩，这两步仅需要较低计算量；最后计算完整的全秩矩阵w3＝x1*w2，并将这个全秩矩阵用于最后的全秩矩阵更新br＝b-w3；图4显示了低秩更新的过程，其中中间计算结果较小的边长揭示了这种低秩更新效率更高的原因。
[0085]
s17：更新之后，结束本次循环，返回步骤s12处进行下一循环；下一循环只会对上一循环中尾矩阵更新区域进行重复操作；直到进入最后一个循环并对最后一个矩阵块进行分解任务；每一次循环中，任务执行区域的划分如图5所示。由于cholesky分解一般执行在对称矩阵上，只需要存储完整矩阵下三角部分。需要上三角矩阵块的场合，可从对应下三角矩阵块转置结果中获得。
[0086]
4)基于cholesky分解结果的正反向传播；本发明中基于cholesky分解结果的正反向传播求解步骤伪代码如下所示：
[0087][0088]
s21：获得上述分布式低秩cholesky分解结果；mpi初始化方式和openmp工作方式同上述分解过程；其中openmp任务队列中包括三角求解任务、向量更新任务、通信任务；
[0089]
s22：某一进程选取某从线程对对角线块进行三角求解任务；当前从线程执行完三角求解任务后，由所有进程的线程进行通信任务同步，然后释放通信任务的依赖；最后执行过三角求解的进程会将求解结果，发送给对应进程；该通信任务的实现方式同步骤s13；
[0090]
s23：各进程选取某一线程对矩阵块进行向量更新任务；该更新使用标准的全秩矩阵乘法即可完成；
[0091]
s24：更新之后，结束本次循环，返回步骤s22处，进行下一循环；下一循环只会对上一循环中尾矩阵区域进行重复操作；直到进入最后一个循环并求解出最后一部分解向量；全部完成后，主进程开始收集其他进程上的解向量，组合成一个完整解向量；这些数据分布在改进的二维数据布局中，对角线数据块所在进程；
[0092]
步骤(4)：主进程利用完整解向量，根据标准矩阵乘法计算物理值预测向量z1＝∑2*x，即为近海气候预测值；预测值向量可以和观测值向量合并为新的完整数据向量，并应用于其他研究。