快速投影Hessian矩阵算法
基于Hessian矩阵的冠状动脉中心线的跟踪算法许燕
ISSN 1000-0054CN 11-2223/N 清华大学学报(自然科学版)J T singh ua Un iv (Sci &Tech ),2007年第47卷第6期2007,V o l.47,N o.634/36889-892基于Hessian 矩阵的冠状动脉中心线的跟踪算法许 燕1, 胡广书1, 商丽华2, 耿进朝2(1.清华大学生物医学工程系,北京100084;2.清华大学第一附属医院,北京100016)收稿日期:2006-06-12基金项目:裕元医学科学研究基金资助项目作者简介:许燕(1980—),女(汉),浙江,博士研究生。
通讯联系人:胡广书,教授,E -mail :hgs -d ea @ts inghu a .edu .cn摘 要:冠状动脉血管造影对医生临床诊断心血管疾病非常有帮助。
在冠状动脉造影图像上对冠脉中心线的正确提取是冠脉边缘定位和三维重建的基础。
该文提出的中心线提取方法是对经典Sun 算法的改进,此方法结合了Hessian 矩阵特征向量和Canny 算子来进行准确提取。
实验结果表明:结合Hessian 矩阵特征向量的方向,较好地解决了冠脉造影图像中冠脉曲率变化剧烈而跟踪不准确的情况;结合Can-ny 算子的半径计算,较好地解决了冠脉造影图像中冠脉重叠和交叉出现时错误跟踪的情况。
与经典Sun 算法相比,本方法有很好的鲁棒性和较高的准确性。
关键词:图像分析;冠状动脉;中心线;Hessian 矩阵;特征向量中图分类号:T N 911.73文献标识码:A文章编号:1000-0054(2007)06-0889-04Adaptive tracking extraction of vessel centerlines in coronary arteriogramsusing Hessian matrixXU Yan 1,HU G uangshu 1,SHA NG Lih ua 2,G ENG Jinzh ao 2(1.Department of Biomedical Engineering ,T s inghua University ,Beij ing 100084,China ;2.The No .1Hospital Attached to Tsinghua University ,Beij ing 100016,China )Abstract :A method for accurate ex tr action of the coronar y arter ial centerline was pres ented for au tom ated positioning of th e coronary ves sel border s and 3-D reconstruction to aid clin ical diagnoses of cardiovas cular d iseases.T he meth od is an improvement of the Sun ar ith metic.Th e m ethod implemen ts an accu rate ex traction using Hes sian matrix eigenvectors and the C ann y operator.The tes t results s how that the determination of the centerline point directions by th e Hess ian matrix eigenvector resolves tr acking inaccur acies resulting from ab rupt changes of the arterial curvatu re.T hecalculation of th e artery radiu s by th e Canny operator r esolves inter ruptions of the tr acking res ulting from the s uper pos ition and overlappin g of the coronary artery.T he m ethod has betterrobus tn es s and accu racy than the Sun algorithm.Key words :image analysis ;coronar y artery;center line;Hess ianmatrix;eigenvector冠状动脉血管造影可以辅助医生准确诊断心血管疾病。
改进SURF算法在图像汉字识别中的应用
改进SURF算法在图像汉字识别中的应用孟伟;钟娜【摘要】针对复杂背景下汉字匹配准确率较低的问题，提出一种改进的SURF算法。
该算法利用灰度分级的字符分割方法，先进行灰度分割增强图像的对比度，采用灰度分级树将图像中的所有像素处理为树的模式进行计算，根据灰度分级确定主节点，根据主节点的级别所对应的灰度值对图像进行分割。
同时，根据汉字结构的特殊性，取消了SURF算法的旋转不变性。
实验结果表明，与未使用改进的SURF 算法相比，对图像质量较差的文本图像，改进的SURF算法能有效地提高其匹配的准确率。
%Aiming at the low matching accuracy of Chinese characters, an improved algorithm of SURF is presented. The algorithm is based on gradation character segmentation. Contrast of image is enhanced by using gray level segmentation, and then with the gray level classification tree, all pixels in the image are processed to the tree model. According to the gray level classification, the main node is determined. Grey level corresponding to the main node level is used in image segmentation. According to the particularity of Chinese characters, the rotation invariance of SURF algorithm is cancelled. Experimental results show that the improved algorithm can improve the matching accuracy effectively, especially for text image of poor quality.【期刊名称】《计算机工程与应用》【年(卷),期】2015(000)012【总页数】5页(P156-160)【关键词】复杂背景;汉字匹配;快速鲁棒特征(SURF)算法;灰度分级;字符分割【作者】孟伟;钟娜【作者单位】北京林业大学信息学院，北京 100083;北京首钢自动化信息技术有限公司，北京 100043【正文语种】中文【中图分类】TP3911 引言近年来，汉字识别一直是模式识别等相关领域内长期的研究热点[1]。
Hessian矩阵【转】
Hessian矩阵【转】在数学中，海塞矩阵是⼀个⾃变量为向量的实值函数的⼆阶偏导数组成的⽅块矩阵，⼀元函数就是⼆阶导，多元函数就是⼆阶偏导组成的矩阵。
求向量函数最⼩值时可以使⽤，矩阵正定是最⼩值存在的充分条件。
经济学中常常遇到求最优的问题，⽬标函数是多元⾮线性函数的极值问题，尚⽆⼀般的求解⽅法，但判定局部极⼩值的⽅法就是⽤hessian矩阵：在x0点上，hessian矩阵是负定的，且各分量的⼀阶偏导数为0，则x0为极⼤值点。
在x0点上，hessian矩阵式正定的，且各分量的⼀阶偏导数为0，则x0为极⼩值点。
矩阵是负定的充要条件是各个特征值均为负数。
矩阵是正定的充要条件是各个特征值均为正数。
函数如下：如果f所有的⼆阶导数都存在，那么f的海塞矩阵即为：H(f)ij(x) = D i D j f(x)，即（也有⼈把海⾊定义为以上矩阵的）海赛矩阵被应⽤于⽜顿法解决的⼤规模优化问题。
性质对称性：如果函数f在D区域内⼆阶连续可导，那么f海塞矩阵H（f）在D内为对称矩阵。
原因是：如果函数f连续，则⼆阶偏导数的求导顺序没有区别，即：则对于海塞矩阵H（f）,有，所以为对称矩阵。
多元函数极值的判定如果实值多元函数⼆阶连续可导，并且在临界点M(x i)(其中i=1,2,...,n，并且X i已知)处梯度（⼀阶导数）等于0，即，则M为驻点。
仅通过⼀阶导数⽆法判断在临界点M处是极⼤值还是极⼩值。
记f在M点处的⿊塞矩阵为H(M)。
由于f在M点处连续，所以H(M)是⼀个的对称矩阵，对于H(M),由如下结论：如果H(M)是，则临界点M处是⼀个局部的极⼩值。
如果H(M)是，则临界点M处是⼀个局部的极⼤值。
如果H(M)是，则临界点M处不是极值。
hessian矩阵计算公式
海森矩阵（Hessian matrix）是一个二阶导数矩阵，用于描述函数在某一点的局部曲率。
对于一个二元函数f(x, y)，其海森矩阵计算公式为：
H(f) = ∂²f/∂x² * ∂²f/∂y² - ∂²f/∂x∂y * ∂²f/∂y∂x
其中，∂²f/∂x² 和∂²f/∂y² 是函数在某一点的二阶偏导数，∂²f/∂x∂y 是函数在某一点的混合偏导数。
对于多元函数，海森矩阵是一个对称矩阵，行和列都对应于函数的每个自变量。
在二元函数的情况下，海森矩阵是一个2x2矩阵。
需要注意的是，海森矩阵在某些情况下可能不是唯一的，对于一个函数在不同的点可能会有不同的海森矩阵。
此外，海森矩阵也可以通过数值方法进行计算，例如有限差分法或有限元法等。
牛顿法的hessian矩阵
牛顿法的hessian矩阵牛顿法是一种常用的最优化算法，可以用于求解函数的极小值或者最小值。
其基本思想是通过利用函数的二阶导数信息来更新搜索方向，从而提高迭代的收敛速度。
在牛顿法中，关键的一步就是求解函数的Hessian矩阵，本文将介绍如何计算Hessian矩阵。
1. Hessian矩阵的定义Hessian矩阵是一个二阶偏导数组成的方阵，用$\mathbf{H}$表示，其元素如下所示：$$H_{i,j}=\frac{\partial^2 f}{\partial x_i\partial x_j}$$其中，$f$是要进行优化的目标函数，$x_i$和$x_j$是函数的自变量。
2. 牛顿法中的Hessian矩阵在牛顿法中，我们希望通过Hessian矩阵来计算函数的搜索方向。
具体来说，我们会用Hessian矩阵的逆矩阵来计算搜索方向，从而使得每一次迭代的变化方向更加准确。
3. 求解Hessian矩阵的方法求解Hessian矩阵的方法有很多种，这里介绍两种常用的方法。
（1）符号计算法符号计算法是一种基于数学计算的方法，可以精确地计算函数的Hessian矩阵。
这种方法的缺点是计算复杂度比较高，适用于简单的函数，对于复杂的函数而言，其计算时间会非常长。
（2）数值计算法数值计算法是一种基于数值计算的方法，通过对函数进行数值求导，得到数值近似的Hessian矩阵。
这种方法的优点是计算速度比较快，适用于复杂的函数，缺点是精度相对较低。
4. 总结Hessian矩阵是牛顿法中非常重要的一步，其可以帮助我们更加准确地计算函数的搜索方向。
在实际应用过程中，我们可以选择符号计算法或者数值计算法来求解Hessian矩阵，具体方法要根据实际情况来选择。
找特征点的算法SIFT和SURF算法
找特征点的算法SIFT和SURF算法SIFT算法和SURF算法是用于图像特征点的检测与描述的两种经典算法。
它们在图像处理、计算机视觉和模式识别等领域得到广泛应用。
下面将分别介绍SIFT算法和SURF算法，并对其原理和应用进行详细阐述。
一、SIFT算法（Scale-Invariant Feature Transform）SIFT算法是由Lowe于1999年提出的一种用于图像特征点检测与描述的算法。
它通过分析图像的局部特征来提取与尺度无关的特征点，具有尺度不变性、旋转不变性和仿射不变性等优点。
1.特征点检测SIFT算法首先通过高斯差分金字塔来检测图像中的特征点。
高斯差分金字塔是由一系列模糊后再进行差分操作得到的，通过不同尺度的高斯核函数对图像进行卷积，然后对结果进行差分运算，得到图像的拉普拉斯金字塔。
在拉普拉斯金字塔上，通过寻找局部最大值和最小值来确定特征点的位置。
2.特征点描述在确定特征点的位置后，SIFT算法使用梯度直方图表示特征点的局部特征。
首先，计算特征点周围邻域内每个像素点的梯度幅值和方向，然后将邻域分为若干个子区域，并统计每个子区域内的梯度幅值和方向的分布，最后将这些统计结果组合成一个向量作为特征点的描述子。
3.特征点匹配SIFT算法通过计算特征点描述子之间的欧式距离来进行特征点的匹配。
欧式距离越小表示两个特征点越相似，因此选择距离最近的两个特征点作为匹配对。
二、SURF算法（Speeded Up Robust Features）SURF算法是由Bay等人于2024年提出的一种在SIFT算法的基础上进行改进的图像特征点检测与描述算法。
它通过加速特征点的计算速度和增强特征点的稳定性来提高算法的实时性和鲁棒性。
1.特征点检测SURF算法使用Hessian矩阵来检测图像中的特征点。
Hessian矩阵是图像的二阶导数矩阵，通过计算Hessian矩阵的行列式和迹来确定图像的局部最大值和最小值，从而找到特征点的位置。
Hessian矩阵
Hessian矩阵1. Jacobian在向量分析中, 雅可⽐矩阵是⼀阶偏导数以⼀定⽅式排列成的矩阵, 其⾏列式称为雅可⽐⾏列式. 还有, 在代数⼏何中, 代数曲线的雅可⽐量表⽰雅可⽐簇：伴随该曲线的⼀个代数群, 曲线可以嵌⼊其中. 它们全部都以数学家卡尔·雅可⽐(Carl Jacob, 1804年10⽉4⽇－1851年2⽉18⽇)命名；英⽂雅可⽐量”Jacobian”可以发⾳为[ja ˈko bi ən]或者[ʤəˈko bi ən].雅可⽐矩阵雅可⽐矩阵的重要性在于它体现了⼀个可微⽅程与给出点的最优线性逼近. 因此, 雅可⽐矩阵类似于多元函数的导数.雅可⽐⾏列式如果m = n, 那么FF是从n维空间到n维空间的函数, 且它的雅可⽐矩阵是⼀个⽅块矩阵. 于是我们可以取它的⾏列式, 称为雅可⽐⾏列式.在某个给定点的雅可⽐⾏列式提供了在接近该点时的表现的重要信息. 例如, 如果连续可微函数FF在pp点的雅可⽐⾏列式不是零, 那么它在该点附近具有反函数. 这称为反函数定理. 更进⼀步, 如果pp点的雅可⽐⾏列式是正数, 则FF在pp点的取向不变；如果是负数, 则FF的取向相反.⽽从雅可⽐⾏列式的绝对值, 就可以知道函数FF在pp点的缩放因⼦；这就是为什么它出现在换元积分法中.对于取向问题可以这么理解, 例如⼀个物体在平⾯上匀速运动, 如果施加⼀个正⽅向的⼒FF, 即取向相同, 则加速运动, 类⽐于速度的导数加速度为正；如果施加⼀个反⽅向的⼒FF, 即取向相反, 则减速运动, 类⽐于速度的导数加速度为负.2. 海森Hessian矩阵在数学中, 海森矩阵(Hessian matrix或Hessian)是⼀个⾃变量为向量的实值函数的⼆阶偏导数组成的⽅块矩阵, 此函数如下：2), 最优化在最优化的问题中, 线性最优化⾄少可以使⽤单纯形法(或称不动点算法)求解, 但对于⾮线性优化问题, ⽜顿法提供了⼀种求解的办法. 假设任务是优化⼀个⽬标函数ff, 求函数ff的极⼤极⼩问题, 可以转化为求解函数ff的导数f′=0f′=0的问题, 这样求可以把优化问题看成⽅程求解问题(f′=0f′=0). 剩下的问题就和第⼀部分提到的⽜顿法求解很相似了.这次为了求解f′=0f′=0的根, 把f(x)f(x)的泰勒展开, 展开到2阶形式：。
matlab求解机器人的hessian矩阵算法
matlab求解机器人的hessian矩阵算法如何用MATLAB求解机器人的Hessian矩阵算法引言：机器人在各个领域扮演着重要的角色，特别是在工业自动化和服务行业中。
为了使机器人能够更加智能和高效地完成任务，我们需要对其控制算法进行优化。
其中，Hessian矩阵是一种常用的工具，可以用于机器人的姿态控制和目标优化等问题。
本文将详细介绍如何用MATLAB求解机器人的Hessian矩阵算法，以帮助读者深入理解并应用该算法。
第一步：了解Hessian矩阵的基本概念和作用Hessian矩阵是多元函数的二阶偏导数构成的方阵，用于描述函数局部极值点的性质。
在机器人姿态控制中，Hessian矩阵可以用来衡量机器人在给定位置周围的姿态变化情况，从而优化控制算法，使机器人更稳定地工作。
第二步：构建机器人姿态模型在实际应用中，机器人的姿态模型通常使用欧拉角表示，即通过三个旋转（roll、pitch和yaw）来定义机器人的姿态。
假设机器人姿态表示为[x, y, theta]，其中x和y表示机器人的位置坐标，theta表示机器人的旋转角度。
我们可以通过编写MATLAB代码来构建机器人姿态模型，具体内容如下：MATLABfunction robot_pose = robot_model(x, y, theta)robot_pose = [x; y; theta];end第三步：计算Hessian矩阵的偏导数根据Hessian矩阵的定义，我们需要计算机器人姿态模型的二阶偏导数。
在MATLAB中，我们可以使用符号计算工具箱（Symbolic Math Toolbox）来计算函数的二阶偏导数。
具体步骤如下：1. 导入符号计算工具箱MATLABsyms x y theta2. 计算机器人姿态模型MATLABrobot_model = robot_model(x, y, theta);3. 计算姿态模型对x的一阶偏导数MATLABd_model_dx = diff(robot_model, x);4. 计算一阶偏导数对x的一阶偏导数MATLABd2_model_dx2 = diff(d_model_dx, x);5. 以此类推，计算其他变量（y和theta）的一阶和二阶偏导数第四步：将计算结果转换为数值进行计算由于计算的结果是符号表达式，我们需要将其转换为数值进行实际计算。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
文章编号:1671 1114(2009)03 0018 04快速投影Hessian 矩阵算法收稿日期:2008 03 10基金项目:天津市高校发展基金项目(20060402)作 者:汤大林(1965 ),男,高级工程师,主要从事数学建模及应用方面的研究.汤大林(天津理工大学理学院,天津300191)摘 要:分析了求解等式约束非线性规划问题的投影H essian 矩阵算法,找出了算法两步Q 超线性收敛的原因,并用BY RD 的例子说明此算法的收敛效果较差,即甚至不是线性收敛;对算法进行了合理的改进,并用改进后的算法求解BY RD 问题,得到了满意的收敛效果,即Q 超线性收敛.借助数值试验验证了改进算法的快速收敛性.关键词:等式约束非线性规划;投影H essian 矩阵算法;超线性收敛中图分类号:O 221.2 文献标识码:AQ uick projection method with H essian matrixT AN G Dalin(School of Science,T ianjin University of Techn ology,T ian jin 300191,China)Abstract:T he project ion method wit h H essian mat rix used to so lve nonlinear prog ramming w ith equality co nstr aint is analyzed and the r easo n w hy the method is superlinear conver gent by two steps is found o ut.Its bad converg ent effect at linear ity is illuminated by BYRD's example.T he method is impr ov ed and quickly super linear conver gence o f the impr ov ed metho d is illuminated using BY RD's ex ample.T he quickly co nv erg ent effect o f t he impro ved method is verified by a numerical experiment.Key words:nonlinear pr og ramming w ith equality constra int ;project ion method wit h H essian mat rix ;super linearconver gence1 投影Hessian 矩阵算法的缺点考察等式约束非线性规划问题: m in x Rnf (x ),约束c(x)=0,(1)其中,目标函数f (x ):R n !R,约束c(x):R n !R m 是二次可微函数,且m ∀n,即m 个等式约束.为叙述方便,引入如下记号:x =(x (1),x (2),#,x (n)),c(x)=(c (1)(x ),c (2)(x ),#,c (m)(x ))T, g(x)= f (x )=( f x (1), f x (2),#, f x (n))T,A(x)= c(x)=c (1)x (1)c (2)x (1)# c (m) x (1) c (1) x (2) c (2) x (2)# c (m) x (2) ! c (1) x (n)c (2) x (n)#c (m) x (n ),L (x, )=f (x )-∃mi=1(i)c (i)(x ),其中, (i)为拉格朗日乘子,i =1,2,#,m.将A(x)QR 分解为A(x)=(y(x),z(x))(R(x)O ),其中y (x),z(x)均为n 阶正交矩阵,R(x)为m 阶上三角矩阵.V ol.29N o.3Jul.2009第29卷 第3期2009年7月 天津师范大学学报(自然科学版)Jour nal of T ianjin N orma l U niver sity (N atural Science Edit ion)定理1[1] 1)x*是问题(1)的解的一阶必要条件是存在拉格朗日乘子 ,满足q(x, )= x, L(x, )=(g(x)-A(x)c(x))=0,(2)记W(x, )= 2x L(x, )= 2f(x)-∃m i=1 (i) c2(i)(x).2)x*是问题(1)的解的二阶必要条件是投影H essian阵Z T(x)W(x, )Z(x)半正定.3)x*是问题(1)的解的二阶充分条件是式(2)成立,且Z T(x)W(x, )Z(x)正定.假设x*是问题(1)的局部解,且A*=A(x*)列满秩.为求解问题(1),T.F.Colemen等设计了投影H essian矩阵修正算法[2].下面由该算法的简化推导过程说明该算法的缺点.由定理1的1)知,求解问题(1)即求式(2)的零点,可用牛顿法.为了叙述方便,在点(x k, k)处,记g k=g(x k),c k=c(x k),W k=W(x k, k),其他类似.W k-A k A T k Ox kk=-g k-A k kc k,x k+1 k+1=x kk+x kk,∀W k-A k A T k O p kk=-g kc k,x k+1=x k+x k,∀Q T k W k-A kA T k OQ k Q T kx kk+1=-Q T kg kc k,其中,x k=y k p yk +z k p zk,Q T k=(y k,z k)OO I,∀y T k W k y k y T k W k z k-R kz T k W k y k z T k W k z k OR T k O Op ykp zkk+1=-y T k g kz T k g kc k,(3)式(3)是m+n个未知数的线性方程组,由此可知R k p yk =-c k∀得出p yk,(4a)z T k W k y k p yk +z T k W k z k p zk=-z T k g k∀得出p zk,(4b)x k=y k p yk +z k p zk,x k+1=x k+x k,(4c)R k k+1=y T k g k+y T k W k x k,(4d)式(4d)为 k+1的二阶估计,也可由其一阶估计即最小二乘估计给出,即由m in%A k+1 -g k+1%2给出,k+1=R-1k+1y T k+1g k+1,(5)将式(5)代入式(4d)得z T k W k y k z T k W k z kR k Op ykp zk=-z T k g kc k,(6)由定理1的二阶必要条件知z T k W k z k为半正定阵,其计算量很大.若用一个正定矩阵B k近似z T k W k z k似乎是可行的方法,但是这样又会导致z T k W k y k无法处理,解决该问题的最简单的办法就是直接略去z T k W k y k,如此便得到投影H essian矩阵算法.算法1 投影H essian矩阵算法.Step1 1)选择精度!,置k=0.2)置初始点x0,B0,其中B0为正定近似投影阵.3)分别计算f0,g0,c0,A0,分解A0=(y0 z0)R0O.Step2 解方程组R k p yk=-c k(7)得出p yk,解方程组B k p zk=-z T k g k(8)得出p zk,置x k+1=x k+y k p yk+z k p zk.Step3 计算f k+1,g k+1,c k+1,A k+1,分解A k+1=(y k+1 z k+1)R k+1O,计算%z T k+1c k+1%2.Step4 如果%z T k+1c k+1%2<!,则结束.否则,用DFT或BFGS修正公式计算B k+1,并置k=k+1,返回Step2.本算法被证明[3 4]在适当的假设条件下两步Q 超线性收敛于问题(1)的一个局部极小值点x*,即%x k+1-x*%2%x k-x*%2!0,k!&.关于此算法能否一步Q 超线性收敛,下面通过BYRD的一个例子给出了否定的回答.例min f(x)=12x2(1)-∀x(1)x(2)+12x2(2)-(x(1)-∀)33∀,∀为常数,约束c(x)=12-x(2)-1=0.(9)由于约束相当于x(2)=1,不难看出本例的最优解为x*=(∀,1)T.∋19∋第29卷 第3期 汤大林:快速投影H essian矩阵算法现用算法1求解,取x0=(∀,1+#)T,B0为投影H essian矩阵,于是x1=(∀+∀#,1+#2)T,x2=(∀,1+#4)T,x3=(∀+∀#4,1+#8)T,#,当k为偶数时,x k=(∀,1+#2k)T,x k+1=(∀+∀#2k,1+#2k+1)T,易知%x k-x*%2=#2k,%x k+1-x*%2=#2k∀2+#2k,%x k+1-x*%2%x k-x*%2=∀2+#2k,(10)由式(10)知,只有当|#|∀1时算法才收敛,但只是局部收敛.另外当∀(0时,算法即使收敛也不是Q 超线性收敛,当∀>1时甚至不是线性收敛.2 算法的改进从上例看出算法1的收敛效果不能令人满意,这是因为算法1在式(6)中简单地去掉了z T k W k y k,这样势必导致不精确.重新考虑z T k W k y k,在式(6)中将其移到等式右边,并做近似,则式(6)变为R k p yk=-c k,z T k W k z k-p zk =-z T k g k-z T k W k y k p yk)-z T k x L(x k+y k p yk , k),(11)如此处理比直接从式(6)中简单地去掉z T k W k y k提高了一阶精度,从而得到算法1的改进算法.算法2 改进的投影H essian矩阵算法Step1 同算法1.Step2 解方程组R k p yk =-c k,得出p yk,置x∗k=x k+y k p yk.解方程组B k p zk =-z T k x L(x∗k, k),得出p zk,置x k+1=x k+y k p yk +z k p zk.Step3 同算法1.Step4 同算法1.下面用算法2来求解BYRD的例子,说明算法2的优点.易知A k=( 01(2-x(2)k)2)=0110(1(2-x(2)k)20),y k=01,z k=1.取B k=2fx2(1)k=1-2(x(1)k-∀)∀,g k=x(1)k-∀x(2)k-(x(1)k-∀)2∀-∀x(1)k+x(2)k.第一次迭代:x0=(∀,1+#)T,B0=1,R0=1(1-#)2,p y=-(R T0)-1c0=-(1-#)2(11-#-1)=#(#-1),x∗0=(∀,1+#)T+(0,1)T#(#-1)=(∀,1+#2)T,g∗0=(-∀#2,1+#2-∀2)T,p z=-B-10(1,0)(-∀#2,1+#2-∀2)T=∀#2,x1=x∗0+z0p z=(∀+∀#2,1+#2)T.第二次迭代:B1=1-2#2,R1=1(1-#2)2,c1=#21-#2,p y1=-(R T1)-1c1=-(1-#2)#2,x∗1=(∀+∀#2,1+#2)T+(0,1)T#2(#2-1)=(∀+∀#2,1+#4)T,g∗1=∀#2-2∀#4-∀(1+#2)+1+#4,p z1=-B-11z1g∗1=-∀#2,x2=x∗1+z1p z1=(∀,1+#4)T.第三、四次迭代:x3=(∀+∀#8,1+#8)T,x4=(∀,1+#16)T,#.由归纳法知:x k=(∀,1+#2k)T, k为偶数,(∀+∀#2k,1+#2k)T,k为奇数.故%x k-x*%=#2k, k为偶数,#2k∀+1,k为奇数.%x k+1-x*%2%x k-x*%2=#2k∀2+1, k为偶数,#2k∀2+1, k为奇数,k!&0.(12)式(12)表明,无论∀取何值,只要|#|<1,算法2即为一步Q 超线性收敛.3 数值试验利用BYRD的例子,应用数学软件M athematica6.0编程,分两种情形对两种算法进行数值对比试验.∋20∋天津师范大学学报(自然科学版)2009年7月情形1:取迭代初值x0=(0,1)T,∀=0.2+0.2+3,收敛精度为%x k-x*%∀max(10-9,10-8%x k%)且3f(x k)∀10-9.试验结果见表1.情形2:取迭代初值x0=(∀,1)T,∀=0.2+0.2+3,收敛精度为%x k-x*%∀max(10-9,10-8%x k%)且3f(x k)∀10-9.试验结果见表2.两种情形的对比试验结果均表明:1)同等条件下,算法2比算法1收敛速度快;2)不论∀<1或>1,算法2均收敛,而算法1只在∀<1时收敛.表1 情形1∀算法2(x0=(0,1)T)算法1(x0=(0,1)T)迭代次数x*最优值迭代次数x* 最优值0.22x(1)!0,x(2)!10.513333 5x(1)!0.07779,x(2)!0.01556 0.00594670.42x(1)!0,x(2)!10.553333 5x(1)!0.16473,x(2)!0.06589 0.02224930.62x(1)!0,x(2)!10.620000 5x(1)!0.27502,x(2)!0.16502 0.04327120.82x(1)!0,x(2)!10.713333 5x(1)!0.44287,x(2)!0.35429 0.05428281.02x(1)!0,x(2)!10.833333 23x(1)!1.00000,x(2)!1.00000 01.22x(1)!0,x(2)!10.980000 2x(1)!2.71726,10103,x(2)!8.09698,1018-5.5730120,10309 1.42x(1)!0,x(2)!11.153330 101x(1)!2.28077,10102,x(2)!1.37916,1020-2.8248700,10306 1.62x(1)!0,x(2)!11.353330 101x(1)!1.69299,10103,x(2)!9.95408,1019-1.0109318,103091.82x(1)!0,x(2)!11.580000 101x(1)!1.88499,10103,x(2)!1.10902,1020-1.2403130,103092.02x(1)!0,x(2)!11.833330 101x(1)!1.30476,10103,x(2)!-2.22676,1020-3.7019995,10308 2.22x(1)!0,x(2)!12.113330 101x(1)!1.48876,10103,x(2)!2.12286,1019-4.9995037,10308 2.42x(1)!0,x(2)!12.420000 101x(1)!1.81431,10103,x(2)!1.23266,1020-8.2947708,10308 2.62x(1)!0,x(2)!12.753330 101x(1)!5.28182,10103,x(2)!3.32167,1019-1.8891057,103102.82x(1)!0,x(2)!13.113330 101x(1)!2.77578,10103,x(2)!4.03269,1019-2.5460911,103093.02x(1)!0,x(2)!13.500000 101x(1)!3.20881,10103,x(2)!5.90223,1019-3.6710456,10309表2 情形2∀算法2(x0=(∀,1)T)算法1(x0=(∀,1)T)迭代次数x*最优值迭代次数x* 最优值0.22x(1)!0.2,x(2)!1 0.48 6x(1)!0.07779,x(2)!0.01556 0.00594670.42x(1)!0.4,x(2)!1 0.42 6x(1)!0.16473,x(2)!0.06589 0.02224930.62x(1)!0.6,x(2)!1 0.32 6x(1)!0.27502,x(2)!0.16502 0.04327120.82x(1)!0.8,x(2)!1 0.18 5x(1)!0.44287,x(2)!0.35429 0.05428281.02x(1)!1.0,x(2)!1 0 2x(1)!1.00000,x(2)!1.00000 01.22x(1)!1.2,x(2)!1-0.22 2x(1)!4.99756,10103,x(2)!-3.29194,1018-3.4671348,10310 1.42x(1)!1.4,x(2)!1-0.48 101x(1)!3.81463,10103,x(2)!-2.10787,1018-1.3216242,10310 1.62x(1)!1.6,x(2)!1-0.78 101x(1)!3.58878,10103,x(2)!3.54846,1019-9.6293836,103091.82x(1)!1.8,x(2)!1-1.12 101x(1)!3.64849,10103,x(2)!6.18958,1019-8.9938505,103092.02x(1)!2.0,x(2)!1-1.50 101x(1)!4.09807,10103,x(2)!3.84194,1019-1.1470649,10310 2.22x(1)!2.2,x(2)!1-1.92 101x(1)!3.28328,10103,x(2)!1.57167,1020-5.3626603,10309 2.42x(1)!2.4,x(2)!1-2.38 101x(1)!3.62176,10103,x(2)!1.47605,1020-6.5982074,10309 2.62x(1)!2.6,x(2)!1-2.88 101x(1)!5.14778,10103,x(2)!1.68313,1019-1.7488987,103102.82x(1)!2.8,x(2)!1-3.42 101x(1)!6.50098,10103,x(2)!1.70756,1019-3.2708260,103103.02x(1)!3.0,x(2)!1-4.00 101x(1)!2.19182,10104,x(2)!5.10395,1019-1.1699615,10312参考文献:[1] Broyden C G.T he con vergen ce of a clas s of double rank m inimiz ation algorith ms[J].J In st M ath Appl,1970,6:76 90.[2] C oleman T F,Conn A R.On the local convergence of a qu asiN ew ton m ethod for the n onlinear programm ing problem[J].S IAM J Num er Anal,1984,21:769 775.[3] Fletcher R.Practical methods of optimization[M].NewYork:John Wiley&Sons Inc,1980.[4] Shannon D F.Conditioning of quasi New ton methods for functionmini mization[J].M ath ematics of Com puting,1970,24: 647 656.(责任编校 马新光)∋21∋第29卷 第3期 汤大林:快速投影H essian矩阵算法。