有效的数值算法课件

几种有效的数值算法,报告人：王 武 中科院超级计算中心 Email: 20010年5月,2,Fast Algorithm,n,n,If n=64, this table implies an overall reduction in flops of 160 million, which meets the Moores Law! (doubling in 1.5 year),SciDAC Initiative, DOE, CSGF, 2005,n,3,1. 快速多极子方法,快速多极子方法克服了多粒子模拟中最大的瓶颈：精确计算N个粒子之间通过万有引力或静电力的相互作用（比如星系中的星体，或蛋白质中的分子）需要O(N2)的量级。而FMM达到了O(N)的量级。FMM显著的优点之一是它可以任意调整精度 这种算法通过多极展开（空间的粒子或质点、偶极子，四重极子等等）来近似远处的粒子组对近端的局部粒子组的作用 一个递归划分的空间用来描述随距离增大的更大的组,Fast Multipole Method (FMM) , 1987 L Greengard V Rokhlin, Yale University,4,静电场和引力场 位势的多极子展开 矩阵向量乘积形式,N 体问题,5,FMM的应用综述1,Electromagnetic Scattering,Stellar Clusters,Protein Folding,8,多极子展开,Helmholtz 型 Laplace 型 Gauss 型,9,2. Monte Carlo方法,基于随机模拟的计算方法 确定性问题。建立概率模型，再进行随机抽样观察，其算术平均即为所求解的近似估计。精度可用估计值的标准误差表示。 例如计算积分（多重积分，与维数无关）、计算面积（圆周率） 随机性问题。根据物理情况的概率法则，用计算机抽样试验。 例如中子在介质中的扩散、随机服务系统中的排队、动物的生态竞争（MCNP是Los Alamos实验室开发的大型蒙特卡罗程序，可计算中子、光子和电子的联合输运问题以及临界问题）,/4=S(round)/S(square) =4N(round)/N(square) N= no. of random points,The Metropolis Algorithm for Monte Carlo John von Neumann, et al, Los Alamos Lab, 1946,10,收敛性 误差 对于独立同分布随机序列，由中心极限定理可知,Monte Carlo方法求积分,任何一个积分，都可看作某个随机变量的期望值， 因此，可以用这个随机变量的平均值来近似它,f(x): 概率密度函数,正态分布概率误差,11,3. 单纯形方法,具有约束条件的线性规划问题如何求最优解？ 单纯形方法的基本思想是：从可行域的某一个极点出发，迭代到另一个极点，并使目标函数的值有所改善，直到找出有无最优解时为止 该方法用到了单纯形的概念，单纯形是指N维中的N + 1个顶点的凸包，是一个多胞体(比如直线上的一个线段，平面上的一个三角形，三维空间中的一个四面体） 单纯形法尽管理论上讲效果是指数衰减的，但在实践中却是高度有效的,在线性空间中极大化Z,极大化 并满足约束,Simplex Method for Linear Programming George Dantzig, RAND Corporation, 1947,其中,12,4. Krylov子空间迭代法,Krylov子空间迭代法是用来求解形如Ax=b 的方程组，A是一个NxN 的矩阵，当N充分大时，直接计算x=A-1b变得非常困难。 Krylov方法则巧妙地将其变为如下迭代形式求解。 Kx(i+1)=Kx(i)+b-Ax(i) 这里的K是一个构造出来的接近于A的矩阵，而迭代形式的算法的妙处在于，它将复杂问题化简为逐步的易于计算的子步骤。 当 A是对称矩阵时，Lanczos找到了生成子空间K的正交基的方法。Hestenes 和Stiefel提出了共轭梯度法。后来的GMRES、BiCGStab等改进并扩展了这些算法。 Krylov子空间: Km=spanr0,Ar0,Am-1r0, rm=b-Ax(m) 伴随迭代法的是预条件子：构造M，用迭代法求解M Ax=M b,Krylov Subspace Iteration Method (KSP), 1950 Hestenes, Stiefel, Lanczos; National Bureau of Standards,13,5. QR算法,把一个方阵变换为一个“几乎是”上三角的矩阵（在主对角线下面的一斜列上可能有非零元素）是相对容易的，但要想不产生大量的误差就把这些非零元素消去，就不是平凡的事了。QR 算法正好能达到目的。 基于Householder变换，A可以写成正交矩阵Q 和一个上Hessenberge矩阵H 的乘积:H0=Q0-1AQ0，用下面的QR迭代计算A的本征值，可使迭代次数大大减少。 A0=H0, Am-1=QmRm , Am=RmQm 实际上稠密矩阵的特征问题复杂度都是O(N3),QR Algorithm for Computing Eigenvalues J.G.F. Francis, Ferranti Ltd., London, 1959,14,6. 快速排序方法,它利用古老的分而治之的递归策略来解决排序的问题：挑一个元素作为“主元”、把其余的元素分成“大的”和“小的”两堆(当和主元比较时)、再在每一堆中重复这一过程。 尽管主元选的最差时要做N(N-1)/2 次比较，但对随机分布的数组，快速排序还是具有O( Nlog(N) )的平均速度，其优美的简洁性使之成为计算复杂性的著名的例子。,void QuickSort(SeqList R，int low，int high) /对Rlowhigh快速排序 int pivotpos； /划分后的基准记录的位置， lowpivotposhigh。 if(lowhigh) /仅当区间长度大于1时才须排序 pivotpos=Partition(R，low，high)； /对Rlowhigh做划分 QuickSort(R，low，pivotpos-1)； /对左区间递归排序 QuickSort(R，pivotpos+1，high)； /对右区间递归排序 /QuickSort,Quicksort Algorithms for Sorting, 1962 Tony Hoare, Elliott Brothers, Ltd., London,15,7. 快速Fourier变换,应用数学中意义最深远的算法，无疑是使数字信号处理实现突破性进展的FFT。FFT依赖于分而治之策略，把DFT的复杂度由O(N2) 降到O(Nlog(N),Fast Fourier Transform(FFT), 1965 J Cooley, IBM; J Tukey, AT&T Bell Lab,N=2m,16,连续小波变换 为小波母函数 二进制离散小波 母函数通过伸缩平移变换生成一组正交的离散小波基函数， L2(R)中的函数 f(t) 可通过级数展开获得时域或频域的多尺度信息 离散小波变换 常用的小波正交基 Harr小波，Battle-Lemarie小波 等 WT用于信号分析，图像处理，地震勘探数据处理 等,8. 小波变换,Wavelet Transform (WT), J. Morlet Geophysicist, France,1974,17,9. 遗传算法,遗传算法是一种解全局(局部)优化问题的模拟进化搜索方法，与传统优化方法相比, 它采用概率转移准则进行群体搜索,而且不需要计算目标函数的导数 基本思想:将待优化问题的目标函数理解为生物种群对环境的适应性;将优化变量与生物种群的个体对应;在编码方案、适应度函数及遗传算子确定后，将寻找最优解的过程与生物种群的进化过程类比 该方法在每次迭代中都保留一组候选解，并按某种指标从解群中选取较优的个体，利用遗传算子(选择、交叉和变异)对这些个体进行组合，产生新一代的候选解群，重复此过程，直到满足某种收敛指标为止 应用领域：机器学习、模式识别、神经网络、工业优化控制等方面。解决的问题：旅行商问题、通信网络链接长度的优化问题、铁路运输计划优化、VLSI版面设计、药物分子设计等,Genetic Algorithms (GA), 1975 John Holland, psychologist, John Hopkins,18,GA的基本算子: 选择、杂交、变异 Pseudo-code of GA 1. Choose initial population 2. Evaluate the fitness of each individual in the population 3. Repeat Select best-ranking individuals to reproduce b. Breed new generation through crossover and mutation (genetic operations) and give birth to offspring c. Evaluate the individual fitnesses of the offspring d. Replace worst ranked part of population with offspring 4. Until termination,19,10. 多重网格方法,A hierarchy of discretisations (grids) is considered first The important steps of MG: Smoothing reducing high frequency errors, for example using a few iterations of the GaussSeidel method Restriction downsampling the residual error to a coarser grid Prolongation interpolating a correction computed on a coarser grid into a finer grid 多重网格方法的三要素：光滑，限制，延拓 松弛迭代(如GS、SSOR) 作为光滑 相邻点的加权平均作为限制 分片线性插值作为延拓,Multigrid Method (MG), A. Brandt, 1972,S,E,p,p,p,Level=0,1,2,R,S,R,R,S,3,Ax=b的预条件迭代: M(xk+1-xk)=b-Axk , ek=xk+1-xk , rk=b-Axk MG可作为一种迭代法，也可作为PCG等迭代法的预条件子 Algorithm: MG (Ak, bk, uk, k) %V-循环结构的多重网格方法 Step1. If k is the coarsest level, uk= Inv(Ak) bk, return Step2. uk = Sk (Ak, bk, uk ) % 前光滑化(pre-smoothing), % Sk为光滑(smoothing)算子 Step3. rk = bk - Ak uk % 计算余量 Step4. bk-1 = Rk rk % Rk为限制(restriction)算子 Step5. vk-1 = 0 Step6. MG (Ak-1, b k-1, vk-1, k-1) % 粗网格上的多重网格迭代求解 Step7. uk = uk + Pk vk-1 % 粗网格校正, Pk为延拓(prolongation)算子 Step8. uk = Sk (Ak, bk, uk )