（计算机应用技术专业论文）基于多链mcmc的并行化.pdf

中文摘要 蒙特卡罗马可夫链( m c m c ) 算法是概率推理中的重要方法，在统计学、人 工智能、图形学、统计物理等有重要应用。m c m c 算法对机器的计算速度和存 储空间有很高的要求。并行计算机能够提供更快的计算速度和更大的存贮空间。 现有的m c m c 软件都是基于单处理器的，并行m c m c 算法的研究也还刚刚起 步，离实际应用还有很大距离。 针对这种状况，本文在总结m c m c 算法和并行计算基本技术的基础上，设 计并实现了种基于多链m c m c 的并行化方法，即在并行计算机上应用m c m c 算法的简单多链m c m c 并行模型，具体工作包括： 提出了在并行计算机上应用m c m c 算法的简单多链m c m c 并行模型。该 模型包括输入、数据预处理、并行多链采样、多链收敛性诊断、采样数据处理及 输出六大模块。 以j a g s 为采样引擎，通过c 语言编程和m p i 并行程序设计，实现了该模 型。 开展了实验研究。对几个典型问题，分别应用b u g s 、j a g s 软件及本多链 并行m c m c 模型，对结果作了对比、分析。实验充分说明，简单多链m c m c 并行模型在不减少精度的情况下，能极大地提高计算速度。 关键词：m c m c 并行化j a g sm p i a b s t r a c t m o n t ec a r l om a r k o vc h a i n ( m c m c ) a l o g r i t h m sa r ei m p o r t a n tt op r o b a b i l i s t i c i n f e r e n c e t h e r ea rem a n ya p p l i c a t i o n sf o rm c m ci ns t a t i s t i c s ，a r t i f i c i a li n t e l l i g e n c e ， g r a p h i c s ，s t a t i s t i c a lp h y s i c sa n ds oo n m c m ca l g o r i t h m sh a v eas t r i c tr e q u i r e dt ot h e s p e e da n ds t o r a g eo ft h ec o m p u t e r p a r a l l e lc o m p u t e r sh a v em o r ef a s ts p e e da n dm o r e l a r g es t o r a g e t h ee x i s t i n gs o f t w a r e sf o rm c m c a r ea l m o s tb a s e do i ls i n g l e p r o c e s s t h e r eh a sb e e nf e w e rr e s e a r c h sf o rp a r a l l e lm c m c a l g o r i t h m s ，f a rf r o mp r a c t i c a l a p p l i c a t i o n s i nv i e wo ft h i ss i t u a t i o n , a f t e rc o n c l u d e dm c m ca l g o r i t h m sa n dt h eb a s i c t e c h n o l o g i e so fp a r a l l e lc o m p u t i n g ，t h ep a p e rp r o p o s e sa n di m p l e m e n t sam e d e lf o r a p p l y i n gm c m ca l g o r i t m si np a r a l l e lc o m p u t e r , c a l l e ds i m p l ep a r a l l e l i z a t i o no f m c m cb a s e do nm u l t i c h a i n f i r s tw ep r o p o s ea n da n l y s et h em o d e l t h em o d e lh a s s i xm o d u l e s ，w h i c hi sc a l l e di n p u tm o d u l e ，d a t ap r e p r o c e s s i n gm o d u l e ，p a r a l l e lc h a i n s s a m p l i n gm o d u l e ，m u l t i p l ec h a i n sc o n v e r g e n c eo ft h ed i a g n o s i sm o d u l e ，d a t a p r o c e s s i n gm o d u l ea n de x p o r tm o d u l e t h e nu s i n gj a g ss o f t w a r ea ss a m p l i n ge n g i n e ， w ei m p l e m e n tt h em o d e lw h i t ht h ec p r o g r a m m i n gl a n g u a g ea n dm p il i b r a r yf o r p a r a l l e lp r o g r a m m i n g f i n a l l yw ec o m p a r eo u rp a r a l l e lm o d e lw i t hb u g ss o f t w a r e a n dj a g ss o f t w a r eb ys o m ee x p e r i m e n t s t h ee x p e r i m e n t ss h o wt h a tt h es i m p l e p a r a l l e l i z a t i o no fm c m cb a s e do nm u l t i - c h a i nm o d e lc a ng r e a t l yi m p r o v et h es p e e d o fc o m p u t i n gw i t ht h es a m ea c c u r a c y k e yw o r d s ：m c m cp a r a l l e l i z a t i o nj a g sm p i 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的 研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表 或撰写过的研究成果，也不包含为获得叁盗盘堂或其他教育机构的学位或证 书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中 作了明确的说明并表示了谢意。 学位论文作者签名：走林辱 签字日期：二口p 年占月千日 学位论文版权使用授权书 本学位论文作者完全了解基洼盘鲎有关保留、使用学位论文的规定。 特授权苤鲞盘堂可以将学位论文的全部或部分内容编入有关数据库进行检 索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校 向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名：套稀辱 导师签名： 签字日期：2 。诅年多月哥日签字日期弘占年月厂日 第一章绪论 第一章绪论 本章介绍论文的研究背景和主要工作。 1 1 工作背景 在高维空间上，处理积分或者求和问题主要有三种方法：1 ) 数值方法；2 ) 变分法；3 ) m c m c 算法。 数值方法和变分法有很好的数学理论基础，在精度控制和误差分析上较为方 便；但随着维数增加，传统的数学分析往往无法进行，或者时间和空间上的花费 超过了现有机器的能力。这时，基于马尔可夫链的蒙特卡罗方法( m c m c ) 成为 科学工作者的首选【1 】。 蒙特卡罗方法也称为统计模拟方法，指根据给定概率分布，产生随机数，利 用随机数进行数值模拟的方法。根据特定概率密度7 r ( x ) 产生随机数的方法一般 有四种：1 ) 反函数法，一般是针对一维分布上存在有闭函数形式的概率分布。2 ) 拒绝采样，取m 和易于直接采样的概率密度1 ，( x ) ，使得满足万( x ) m y ( x ) ，依据 _ - ，、 1 ，( x ) 对x 进行采样，并以概率掣接受x 。拒绝采样有两个不足，一是m 可 m v ( x ) 能不存在；二是效率比较低，平均m 个样本点中只接受一个。3 ) 重要性采样， 设f ( x ) 是任意函数，万( x ) 是其概率分布密度，求f ( x ) 的期望；可以根据公式 f f ( x ) ，r ( x ) d x ：f s ( x ) 掣v ( x ) 出，选择易于直接采样的概率密度1 ，( 石) ，依据y ( 工) 对 。 ，l 驯 _ - ，、 x 采样，给每个样本x 赋予权值厂( x ) 掣掣，再运用蒙特卡罗方法计算即可。往往 v 【列 选取接近厂( x ) 万( x ) 的参考分布，这样误差较小。4 ) m c m c 算法，这是一种基于 m a r k o v 链的采样机制，能够自动搜索概率分布较大的区域，所以能把更多的计 算资源用在关键区域。m c m c 算法要求选择合适的转移矩阵，建立以目标分布 万( x ) 为平稳分布的非周期不可约马尔可夫链 x ( f 。非周期不可约马尔可夫链有 ， 一个重要性质：从任意初始状态，运行足够长的时间，最终会达到唯一的平稳状 态，在平稳状态下采样所得到样本的统计性质与直接从目标分布采样得到的样本 的统计性质一致。在高维空间上，如果概率分布比较复杂，m c m c 往往是唯一 有效、可行的计算方法f 2 ，3 1 。 第一章绪论 蒙特卡罗方法最初是由s t a r tu l a m ，n e u m a njv 和e n r i c of e r m i 在1 9 4 6 年研 制e n i a c 时引入的。n i c km e t r o p l i s 被蒙特卡罗方法的简单、强大深深吸引。 m e t r o p l i s 在1 9 4 9 和s t a r tu l a m 合作发表了蒙特卡罗方法的第一个公开文档。随 后，m e t r o p l i s 在1 9 5 3 年发展了蒙特卡罗方法，在物理杂志上发表了著名的论文， 介绍了m e t r o p l i s 算法【4 】。m e t r o p l i s 算法是m c m c 算法的最初形式，在科学研究 和实际应用中发挥了巨大的作用，被评为2 0 世纪十大算法之一p j 。 在m e t r o p l i s 开创性的论文之后，物理学界出现了很多相关研究。从推理的 角度，最重要的是h a s t i n g s 在1 9 7 0 年推广的m e t r o p l i s 算法，即现在的 m e t r o p l i s h a s t i n g s ( m - h ) 算法。h a s t i n g s 和他的博士生p e s k u n 研究发现m e t r o p l i s 算法及其扩展形式m e t r o p l i s h a s t i n g s 算法都是一类算法群的特例，现在我们称 这个算法群为m o n t ec a r l om a r k o vc h a i n ( m c m c ) 。 g e m a ns 和g e m a nd 在1 9 8 4 年在图形学上引进了g i b b s 采样；p e a r lj 把 m c m c 运用到了人工智能领域。g i b b s 采样和贝叶斯推理天然结合，它使得贝叶 斯网络上的推理成为可能，有一个b u g s 的自由软件组织专门负责g i b b s 算法的 开发。但直到1 9 9 0 年m c m c 才第一次在统计学中起来重要作用。在神经网络领 域，n e a l 的工作格外重要。另一个引人关注的是完美采样( p e r f e c ts a m p i n g ) ，一 种m c m c 的特殊采样技术，能够无偏差地从给定分布采集样本。这是p r o p pjg 和w i l s o n db 在1 9 9 5 年提出的【6 1 。 m c m c 算法对机器的计算能力、存储空间有很高要求，但现有的m c m c 算 法和m c m c 软件都是基于单个处理器的，在m c m c 的并行化方面，做的工作 还不是很多 7 q o 。另一方面，并行计算机以其强大的计算能力和日趋便宜的价格 吸引了越来越来的理论研究者和实际用户。把m c m c 算法和并行计算机两者结 合，将是非常有益的工作。 r o s e n t h a ljs 在1 9 9 9 年介绍了并行蒙特卡罗方法，作者探讨了并行m c m c 算 法的可能性及其困难。文中分析了在一个分布式环境中，在每台机器上运行一条 m a r k o v 链的可能性并给出了一些观测结剁1 1 】。g l y n n 和h e i d e l b e r g e r q ! e 先前也做过 这方面的研究。 w i l k i n s o n 在2 0 0 4 年介绍了一种并行方法，把状态空间分解成几个相对独立的 类，在m a r k o v 链的每一次更新中，把对该类的处理放在一个处理器上。w i l s o n 给出了这种方法的_ 个具体实例 1 2 1 。但该方法有很大的局限性，最大的困难是得 根据具体问题划分状态空间，很多时候，这种划分并不存在u 引。 b r o c k w e l la e 在2 0 0 6 年给出了一种“提前预取( p r e f e t c h ) 的并行化方法， 能提供l o g ，n 的并行度，n 是处理器个数。但由于各个处理器处理不同的数据， 计算速度不一样，在实验中，当n 大于四时，该算法的效果已经很差【1 4 1 。 2 第一章绪论 综上所述，m c m c 是概率推理中的重要算法，在科研、生产中发挥着重大作 用；但对并行m c m c 的研究还处在探索阶段，还没有通用的方法和技术。如何在 并行计算机上应用m c m c 就成为一个很有意义的课题，理论上能挖掘m c m c 算 法的内在属性，更重要的是能给实际的生产应用带来帮助。 1 2 本文工作 本文对m c m c 算法和现有m c m c 软件进行了全面调研，总结了m c m c 算 法的核心思想。发现已有算法和软件基本都是基于单处理器的，对并行m c m c 算法的研究至今还没有实质性突破。在这种情况下，本文提出了基于并行计算机 的简单多链m c m c 并行模型并编程实现了这一模型。 简单多链m c m c 并行模型主要分为六部分： 一、输入模块，包括一个对问题描述的数学模型文件和已知观察数据的数据 文件。 二、数据预处理模块，主要是生成在具体m c m c 采样时需要的一些文件。 三、并行多链采样模块，在并行计算机的每个核上，独立进行m c m c 采样， 生成多条马尔可夫链的采样样本，实现简单m c m c 并行。 四、多链收敛性诊断模块，应用收敛诊断程序，对前面生成的采样样本进行 收敛性判断。 五、采样数据处理模块，主要功能是收集各个马尔可夫链的数据，在舍去初 始阶段的数据后，把其余数据看做一个整体，运用蒙特卡罗方法，给出未知参数 的统计值。 六、输出模块。输出未知参数的统计值，一般是未知参数的均值、方差和总 的样本数目。 具体实现主要是应用c 语言和m p i 进行并行程序设计。该实现是独立于具体 的并行计算机的，有良好的可移植性。 最后，对几个实际问题做了实验。实验结果说明，本文提出的并行计算机上 的简单多链m c m c 模型可行、有效，能在相同精度下，加快计算速度。 第二章工作基础 第二章工作基础 本章介绍工作的理论和技术基础。m c m c 算法是本文工作的理论基础，j a g s 软件是本文模型实现中的采样引擎，并行技术是模型实现中的关键技术。 2 1m c m c 通过构造独立同分布随机数，计算积分的蒙特卡罗方法，也称为静态蒙特卡 罗方法。一般有两种方法：1 ) 用频率估计概率来计算积分；2 ) 期望法：用样本 函数的平均值估计的期望来计算积分。 但在高维空间，因为产生独立同分布随机数需要的计算量太大，用静态蒙特 卡罗方法处理速度太慢。所以，在高维空间，一般是通过构造马尔可夫链的极限 不变分布，来模拟计算积分，即m c m c ( m a r k o vc h a i nm o n t ec a r l o ) ，也称动态 蒙特卡罗方法。 用马尔可夫链的样本，来对不变分布、g i b b s 分布、g i b b s 场、高维分布或 样本空间非常大的离散分布等作采样，并用以作随机模拟的方法，统称为 m c m c 。 2 1 1 基本概念 本节介绍m c m c 的基本概念。 蒙特卡罗方法，又称为统计模拟方法，指利用随机数进行数值模拟的方法。 蒙特卡罗方法就是将各种随机事件的概率特征( 概率分布、数学期望) 与随 机事件的模拟联系起来，用试验来确定事件的相应概率与数学期望。蒙特卡罗方 法的突出特点是概率模型的解是用试验得到的，而不是计算出来的。模拟任何一 个实际过程，蒙特卡罗方法都需要用到大量的随机数，计算量很大，需要借助计 算机实现。 蒙特卡罗方法解决的基本问题是数值积分。 蒙特卡罗方法的理论基础是概率论中的弱大数定理。 定理2 1 ( 弱大数定理) 设以为f i d 随机变量，7 = e x i ，s n = x i + 五+ + 瓦， 4 第二章工作基础 女u m e i x , i o ，有 p 以+ 。= jx o = o ，五= ，k = f 。) = 尸 以+ 。= i 以= f 。) 。 公式( 2 - 1 ) 上式刻画了马尔可夫链的特性，称为马尔可夫性( 或“无后效性”) 。马尔可 夫链简称马氏链。 定义2 2 ( 转移概率) 称公式( 2 - 1 ) 中的条件概率p k + 。= j l x 一。) 为马氏链 k ，z 0 ) 的一步转移概率，简称转移概率。 第二章工作基础 定义2 3 ( 时齐马氏链) 当马氏链的转移概率p k + = j l x = f 。) 只与状态f ，歹 有关，$ i - 与n 无光时，称马氏链是时齐的，并记岛= 尸 以+ 。= i 以= f ) 。 本文中除非特别声明，所指的马尔可夫链均为时齐的马尔可夫链。 令p = ( 助) = p o o p o lp 0 2p 0 3 p l op l ip 1 2p 1 3 p 2 0p 2 1p 2 2 p 2 3 ： 3 0p 3 1p 3 2p 3 3 称p 为概率转移矩阵，简称为转移矩阵。乃( 玉s ) 有性质 ( 1 ) p u 0 ，i ，j s ； ( 2 ) 岛= 1 ，v i e s 公式( 2 2 ) y e s 定义2 4 ( 随机矩阵) 若一个矩阵的元素具有公式( 2 2 ) 中两条性质，则称此矩 阵为随机矩阵。 记乃0 ) = p ( 以= f ) ，万( ，z ) = ( 乃( ，1 ) ，死( ，z ) ，乃( 忍) ，) ， 万0 ) 表示以时刻以的概率分布向量。称 乃( o ) ，f s ) 为马尔可夫链的初始分 布。 定理2 2f f ( n ) = 万( o ) p ，其中p ”是尸的刀次幂。 定理表明，一个马尔可夫链的特性完全由它的一步转移概率矩阵p 及初始分 布向量万( o ) 决定。 定义2 5 ( n 步转移概率) 称条件概率矽= 尸( 以棚= f 五= f ，厶s ，朋o ， n l 为马氏链的n 步转移概率。 定义2 6 ( 互通) 状态f 可达状态g j ) ，若存在，z o 使得矽 0 ，记为 i 专。若同时有状态专f ，则称f 与，互通，记为ih 。 。 互通是种等价关系。互通的状态构成一个等价类。 定义2 7 ( 不可约) 称马尔可夫链是不可约的，若马尔可夫链只存在一个等价类。 即所有状态都是互通的。 6 第二章工作基础 定义2 8 ( 周期性) 若集合 耽刀1 ，矽 o ) 非空，则称它的最大公约数d = d ( f ) 为状态i 的周期；若d = l ，则称i 是非周期的。若马尔可夫链中的所有状态都是 非周期的，则称马尔可夫链是非周期的。 定义2 9 ( 常返性) 对于任何状态i ，j ，以帕表示从i 出发经刀步后首次到达的 概率，则有刀= o ，刀帕= p 以_ - j ，鼍j ，k = 1 ，2 ，n - i l x o - - i ，刀1 。 令乃= 刀，若乃= l ，则状态为常返状态。 n = l z ，表示从f 出发在有限步内可以到达的概率。 对于常返状态f ，定义“，= ，耀，表示由f 出发再返回f 所需要的平均步数 n = l ( 时间) 。若u ， 。曲，且乃2 虿1 。 即对于非周期、不可约、正常返的马尔可夫链，极限分布就是该马尔可夫链 的唯一平稳分布。 定义2 1 2 ( 时间可逆) 如果对于所有的f ，r t p o = 乃p 都成立，则称马尔可夫 链为时间可逆的。 定理2 5 ：如果一个马尔可夫链是时间可逆的，则该马尔可夫链是非周期、不可 约的。 定理2 4 和定理2 5 是m c m c 算法的理论基础。一般构造非周期、不可约的、 常返的马尔可夫链，使其极限分布是所求的目标分布；这样在经过充分多的转移 步数后，马尔可夫链到达平稳分布，由定理2 4 可知，此平稳分布就是极限分布， 也即目标分布。而非周期、不可约的马尔可夫链一般是通过构造时间可逆的马尔 可夫链得到的【1 6 】。 定理2 6 ( 马尔可夫链的大数定理) ：设 ) 为不可约的、正常返的马尔可夫链， 无论初始分布“如何，如果对平稳分布万，有易【i 厂( 五) i 】 ，则 一l i m l 窆；。f ( x m ) = 萎乃邝) 也【厂( 驯，w p 1 这里的期望是对概率所取的【1 7 1 。 马尔可夫链的大数定理，是应用蒙特卡罗方法处理马尔可夫链生成的样本的 理论基础。大数定理的本质是随机变量的算术平均依概率收敛于常数一即均值 收敛至期望。 2 1 2m c m c 算法的采样 m c m c 是一种应用马尔可夫链的性质在状态空间上生成满足目标分布p ( x ) 的样本的算法。利用非周期、不可约马尔可夫链的极限分布就是其平稳分布这一 性质，以目标分布为极限分布，构造一条非周期、不可约的马尔可夫链。这样， 从任何初始值出发，马尔可夫链都会达到平稳分布也即目标分布p ( x ) 。依据定 理2 。5 ，在实际中，是通过构造满足时间可逆的马尔可夫链来生成非周期、不可 约的马尔可夫链的【l8 1 。 第二苹工作基础 一_ 一。_ _ _ _ 一 m c m c 是一类算法族。其中，最重要的是m e t r o p o l i s h a s t i n g s ( m - h ) 算法 和g i b b s 算法。 m e t r o p o l i s h a s t i n g s 算法是m c m c 算法族中最重要的算法，大部分m c m c 算法都可以看做是m h 算法的特例或者某种形式的扩展【5 1 9 1 。 m h 算法 已知所求的向量参数x = ( 五，x 2 ，吒) ，其概率密度为p ( 力，一般p ( x ) 较复 杂，很难直接生成符合p ( 曲分布的随机数；可以构造以p ( x ) 为平稳分布的非周 期、不可约马氏链n ，x ( m 。设参考分布为q ( x lz ) ，则算法伪代码为： 1 初始化x ( 0 1 。 2 从i = 0 到一l 一采样“qo i 】。 一采样f q ( x l 工7 ) 。 如果“ 彳c x 。) ，工，= l ， 琶黼) 工( = x 否则 r ( f + i ，= x f n 。 m h 算法简单，但需要认真设计参考分布。一般可以通过次优推理或者和 机器学习得到高效的参考分布。 m h 的转移函数是 其中 k m h ( x ( i x o ，作蚰鄙a = e f 嚣煅鼢 布的马尔可夫链工。，利用遍历定理可以由这个马尔可夫链的一条轨道，得到分 布密度万( x ) 的估计，记为，r o ( x ) 。 二、实现高维积分的数值计算。 用车( 黠) 作为积分，= p ( x 胁的估计。 三、用被估参数的b a y e s 分布的抽样样本，来估计参数。 用m c m c 方法对b a y e s 分布的参数作估计，指对于参数0 的b a y e s 分布， 用m c m c 方法生成以它为极限分布的马尔可夫链见；于是只的分布近似b a y e s 分布，可以取刍：堡正生k 譬业( m n o , n o 充分大) 作为参数秒的估计。 四、求解复杂样本空间上的函数极值( 模拟退火) 2 0 1 。 2 2j a g s m c m c 应用领域非常广泛，但若要自己编程实现，并不容易阱】。本节我们 第二章工作基础 介绍两个现有的m c m c 软件。 一、b u g s b u g s ( b a y e s i a ni n f e r e n c eu s i n gg i b b ss a m p l i n g ) - 一个应用m c m c 算法对 复杂统计模型进行贝叶斯分析的软件【2 3 1 。由英国剑桥公共卫生事务所在1 9 8 9 年 开始研发，最初的版本为c l a s s i cb u g s ，现在有两个版本w i n b u g s 、o p e n b u g s 。 w i n b u g s 是b u g s 在w i n d o w 下的版本，而o p e n b u g s 是在g n u 下发布的开 源软件。 b u g s 是实现m c m c 不确定推理方面最好的软件，有详细的使用文档、详 尽的例子，且基本实现上了对常见分布的采样，用户也能在此基础上扩展新的分 布。b u g s 还定义了一种模型描述语言“b u g s ，可以方便地对数学统计模型用 b u g s 语言描述。b u g s 还实现了与统计软件r 的接v i ，这使得用户可以在r 中调用b u g s 。w m b u g s 和o p e n b u g s 都是菜单界面，很方便用户的使用。b u g s 是用结构化p a s c a l ( c o m p o n e n tp a s c a l ) 语言编写。 二、j a g s j a g s ( j u s ta n o t h e rg i b b ss a m p l e r ) = 采用m c m c 算法对于贝叶斯模型进行 分析，职责在于构建贝叶斯层次模型、选择采样算法及产生随机样本【2 5 1 。j a g s 的数据输入、输出和数据的进一步分析都依赖于r 及r 中的其它软件包。 j a g s 是由m a r t y np l u m m e r 开发的开源软件，这不是b u g s 的版本。j a g s 采用了和b u g s 同样的模型描述语言，但是完全重写了代码。j a g s 是用c + + 编 写的，采用了泛型编程。j a g s 在2 0 0 2 年推出了j a g s 0 1 的最早版本，最新的 版本是j a g s 1 0 。j a g s 在功能上是参照b u g s ，但现在还没能实现b u g s 那么 多功能；它的优势是用c + + 编写的开源代码，且能够运行在u n i x 上。目前，j a g s 能运行大多数经典b u g s 例子，但在运行稳定性及运行效率上都较b u g s 差一 些。 开发j a g s 主要目的有三个： 1 ) 为对贝叶斯模型进行实验提供一个平台； 2 ) 可扩展性，允许用户自己添加函数、分布以及采样方法； 3 ) 提供一个在u n i x l i n u x 环境下运行b u g s 语言描述的模型的采样引擎。 在使用上，j a g s 提供了一个命令行接口；用户可以交互地使用j a g s ，也 可以运行j a g s 脚本。 1 2 第二章工作基础 2 3 并行技术 2 3 1 并行计算简介 并行计算是指在并行计算机上，将一个应用分解成多个子任务，分配给不同 的处理器，各个处理器之间相互协同，并行地执行子任务，从而达到加快求解速 度，或者提高求解问题规模的目的【2 6 】。 对于具体的应用问题，采用并行计算技术的主要目的在于两个方面： 1 ) 加快求解问题的速度， 2 ) 提高求解问题的规模。 并行计算的硬件环境是并行计算机。按内存访问模型、微处理器和互联网络 的不同，目前流行的并行处理器可分为对称多处理器共享并行计算机、分布共享 存储并行计算机、机群、星群和大规模并行计算机五类。 并行计算的软件环境一般指操作系统、并行编程环境。u n i x 操作系统几乎 是目前所有高性能并行机采用的标准操作系统。l i n u x 操作系统也是应用较为广 泛的系统，它能提供很好的机群管理功能。 在并行程序设计方面，比较流行的并行编程环境可以分为三类：消息传递、 共享存储和数据并行。下表给出三种并行编程环境的主要特征 表2 1 三种并行编程环境主要特征 2 3 2m p i 编程 本节介绍消息传递并行计算机模型及消息传递标准界面m p i 。 第二章工作基础 一、消息传递分布式存储并行计算机模型 图2 1 是消息传递进程拓扑结构。其中p 表示m p i 进程，m 表示每个进程 的局部内存空间，多个p m ( 进程内存) 模块通过互联网络相互连接，构成一 个分布式存储的进程拓扑结构。 互联 网络： 消息 传递 函数库 图2 1 消息传递进程拓扑结构 如果每个p m 模块替换成处理器，且规定每个处理器只能分配用户程序的 一个进程，则就是消息传递分布存储并行计算机模型。 特别地，消息传递分布式存储并行计算机模型是对具体并行计算机体系结构 的抽象。无论将该模型映射到何种类型的并行计算机，用户面对的都是该模型隐 含的进程拓扑结构，从而简化并行程序设计，增强程序的可移植性。 二、标准消息传递界面m p i ( m a s s a g ep a s s i n gi n t e r f a c e ) 1 9 9 4 年，全球工业、政府和科研应用部门联合推出消息传递并行编程环境 的标准用户界面，它将消息传递并行环境分解为两部分： 1 ) 消息传递函数的标准接口，即m p i 函数库； ” 2 ) 各并行计算机厂商提供的对这些函数的具体实现。 用户只需学习m p i 库函数的标准接口，设计m p i 并行程序，就可在支持 m p i 并行编程环境的具体并行机上执行该程序。 1 4 第二章工作基础 相对其他并行编程环境，m p i 具有许多优点： 1 ) 很好的可移植性，被当前所有并行环境支持； 2 ) 很好的可扩展性，是目前高效率的大规模并行计算最可信赖的平台； 3 ) 有完备的异步通信功能； 4 ) 有精确定义，为并行软件的发展提供了必要的条件鲫。 三、m p i 并行程序 m p i 并行程序通过对m p i 函数的调用来实现特定的并行算法。一个m p i 并 行程序主要由3 部分组成： 1 ) 进入并行环境：调用m p ii n i t 启动并行计算环境，包括在指定的计算结 点上启动构成并行程序的所有进程及构建初始的m p i 通信环境、通信器 m p i c o m m w o r l d 、m p i c o m m s e l f 。 2 ) 主体并行任务：这是并行程序的实质部分，所有需要并行完成的任务都 在这里进行，实现并行算法在并行计算机上的执行过程。 3 ) 退出并行环境：调用m p if i n a l i z e 退出并行环境。一般，退出并行计算 环境后，程序的运行也马上结束。 m p i 并行程序从程序结构上可分为3 种编程模式，包括主从模式、单程序多 数据模式和多程序多数据模式。 第三章简单多链m c m c 并行模型 第三章简单多链m c m c 并行模型 本章介绍了简单多链m c m c 并行模型的理论基础，重点介绍了模型的架构 及其设计原则、实现方法。 3 1 理论基础 3 1 1 简单蒙特卡罗并行 蒙特卡罗方法通常是计算一个巨大数据集的均值。很直接的想法是，划分数 据集成n 部分，用n 个处理器对每一部分计算均值，然后加权取平均值。这样 处理器间通信很少，可以很容易在并行计算机上实现。 假设要计算未知参数“，并且有程序可以生成独立同分布的样本置，z ， 样本均值为m ，样本方差为1 ，。( 如果m “，则m 一“就是偏差。) 标准的蒙特卡罗方法是采样得到，z 个样本五，五，咒( 刀一般是一个极大 1n 的数) ，然后用e = 三y 置来近似估计“。这个估计的偏差是m u ，方差是刀。 n 百 现在，假设我们在并行计算机上求解上述问题。假设并行计算机有c 个处理 器。我们希望能得到更好的估计值或者较少运行时间。 最简单的方法是在这c 个处理器上运行相同的程序。假设处理器，生成n 个 1h 样本研力，曩n ，计算e ，= 二掣力，并把该值传送给主程序。主程序对这c 个 。 刀百 ， 均值取平均，得到最终结果面：三罗巨。一e 的均值为m ，方差为q c n 。 n 百 上述方法称为简单蒙特卡罗并行方法。简单蒙特卡罗并行方法很容易在并行 计算机上实现，计算结果与单机上结果相比，有相同均值和线性减小的方差【l l ，2 8 | 。 3 1 2 简单m c m c 并行 r 在并行计算机上应用m c m c ，主要有两种策略：1 ) 运行一条很长的单链， 在单链内对每次更新实现并行化；2 ) 运行多条较短的马尔可夫链，综合分析结 果 3 】。 策略1 的优点在于充分了利用了并行计算机的运行能力，马尔可夫链能在较 1 6 第三章简单多链m c m c 并行模型 短时间内达到平稳状态；但实际操作中，对单次更新并行，难度极大，不易实现。 目前，主要有两种方案。第一种是w i l k i n s o n 在2 0 0 4 年介绍的状态空间划分 法。把状态空间分解成几个相对独立的类，在马尔可夫链的每一次更新中，把对 该类的处理放在一个处理器上引。但该方法有很大的局限性，最大的困难是得根 据具体问题划分状态空间，很多时候，这种划分并不存在。第二种是b r o c k w e l l 在2 0 0 6 年介绍的“提前预取( p r e f e t c h ) ”的并行化方法。能提供l o g ，n 的并行 度，n 是处理器个数。但在实际中，当n 大于四时，该算法的效果已经很差【l4 1 。 所以，方法1 ，在实际应用中，目前还没有可行的算法。 策略2 的优点在于实现较简单，但浪费了并行计算机的部分运算能力，因为 在每条马尔可夫链上，都必须舍去初始( b u m - i n ) 阶段的采样样本。 策略2 的理论基础在于：对于以目标分布万( z ) 为平稳分布的、非周期、不 可约马氏链，在其达到平稳分布后，虽然样本点间的关联性依旧很大( 马尔可夫 性) ；但是在马尔可夫链达到平稳分布后，样本数据在统计上服从万( x ) 分布。即 在统计上，样本数据可以看做是从万( x ) 采样的独立同分布的数据。 由此，在并行处理器的每个处理器上运行一条马尔可夫链，在舍去初始阶段 的样本后，其余样本数据从统计上可以看做是服从目标分布的独立同分布的数 据。把这些样本数据综合分析，即实现了m c m c 的并行化。这种方法，我们称 之为简单m c m c 并行方法【 2 9 】。 3 1 3m c m c 收敛性 对m c m c 算法中生成的马尔可夫链进行收敛性判断，是应用m c m c 算法 中的重要一部分。 当m c m c 采样的状态空间是离散的，由第二章的定理可以判断马尔可夫链 的收敛性。两个关键概念是非周期( a p e r i o d i c i t y ) 和不可约( i r r e d u c i b i l i t y ) 。这 两个概念可以扩展到连续空间。t i e m e y 在1 9 9 4 年给出了在连续空间上判断 m c m c 收敛的一个定理。 定理3 1 假设尸是非周期、不可约马尔可夫链的转移矩阵，7 1 是该马尔可夫 链唯一的平稳分布，则对几乎所有的x ，有 l p 加( p ) 一万k 专o ， | | i | 。表示全方差距离嗍。 上述定理表明在平稳分布下，非周期、不可约马尔可夫链收敛。 在实际应用中，对m c m c 收敛性的判断，一般采样收敛诊断程序。 1 7 第三章简单多链m c m c 并行模型 3 2 模型构架 本节具体介绍简单多链m c m c 并行模型的构架及各模块功能。 根据简单m c m c 并行的原理，本文提出了在并行计算机上使用m c m c 算法 的简单多链m c m c 并行模型，即在并行计算机上并行地运行多条独立的马尔可夫 链，进行m c m c 采样，再汇集采样样本，统一用蒙特卡罗方法作数据处理。 简单多链m c m c 并行模型主要分为六部分，见图3 1 。 一、输入模块。 输入模块，包括一个对问题描述的数学模型文件和已知观察数据的数据文 件。 二、数据预处理模块。 数据预处理模块主要是生成在具体进行m c m c 采用时需要的一些文件，包 括为每条马尔可夫链中的未知参数设定的初始值文件。 三、并行多链采样模块。 并行多链采样模块，应用简单m c m c 并行的原理，在并行计算机的每个核 上，独立进行m c m c 采样，生成多条马尔可夫链的采样样本。 四、多链收敛性诊断模块。 多链收敛性模块，基于多条马尔可夫链，应用收敛诊断程序，对前面生成的 采样样本进行收敛性判断；只有在判断收敛后，才可以认为马尔可夫链已经达到 平稳状态、样本服从平稳分布( 即目标分布) 。 五、采样数据处理模块。 采样数据处理模块，主要功能是收集各个马尔可夫链的数据，在舍去初始阶 段的数据后，把其余数据看做一个整体，运用蒙特卡罗方法，给出未知参数的统 计值，一般是未知参数的均值和方差。 六、输出模块。 输出模块，输出未知参数的统计值，一般是未知参数的均值、方差和总的样 本数目。 第三章简单多链m c m c 并行模型 并行处理器 ，、 预处理( 生成初始值等) l 上 厂、 并行多链采样 mmmm cccc mmmm cccc 、， ，、 多链收敛性诊断 l 、 ，r、 采样数据处理 l 、， 黜 图3 1 简单多链m c m c 并行模型 1 9 第三章简单多链m c m c 并行模型 3 3 设计与实现 本节介绍模型的总体设计、各模块的详细设计和本文对该模型的一种具体实 现。 3 3 1 模型设计 软件系统设计由总体设计和详细设计两部分组成，总体设计是确定系统的模 块结构，即将系统划分为模块、确定模块间接口和评价模块划分质量。 m c m c 算法一般应用在高维空间上，进行采样的次数一般都在几千万次以 上，对计算机的计算能力和存贮空间有很高的要求。本文总的设计原则就是在不 降低算法精度的前提下，尽可能地利用并行计算机强大的计算能力和更大的存贮 空