（计算机应用技术专业论文）生物数据库搜索和可视化的研究.pdf

浙江大学硕士学位论文 摘要 随着计算机技术的发展，计算机的应用对生物合成制药领域的研究产生 了重大影响。许多新药的研制不再是传统经验型的一个一个地实验合成化合 物，现在新药的合成往往离不开对海量的生化成分数据进行分析和比较。因 此如何在浩瀚的化学数据库中快速地搜索到研究人员所关心的分子结构并 以图象的形式展现出来，成为一个当今的研究热点。 本文在前人的基础上，提出了一种更高效方法来实现大型生化数据库的 搜索和可视化。该算法通过奇异值分解定理( s v d ) 实现分子数据集从高维 到低维映射( 2 d 或者3 d ) ，并采用了截断式牛顿法来最小化误差，从而确保 维持了各个分子间差异的“大小”。在算法中我们也采用了稀疏矩阵的存储 结构来降低存储需求和计算开销，并结合计算机图形技术，将结果可视化地 展现在研究者面前。 本文首先对问题进行了数学建模，描述了问题的数学模型，接着给出了 详细算法描述，实现细节及实验结果。最后介绍了整个软件的图形人机界面。 实验证明整个算法是比较高效的。 关键词： 药物合成，奇异值分解，截断式牛顿，稀疏矩阵，可视化 塑坚奎兰堡主兰垡笙塞 a b s t r a c t a st h ed e v e l o d m e n to ft h ec o m p u t e rs c i e n c et e c h n o l o g y ， 1 t 1 m p a c t s a l m o s ta 1 1 t h es c i e n t i f i cf i e l d s ， i n c l u d i n gt h es y n t h e t i cd r u ga r e a n a 释a d a v s ， m o s to ft h en e wm e d i c i n e sa r en ol o n g e rd e v e l o p e db y t h e t r a d i t i o n a lw a y ，w h i c ht h ec o l p o u n d sa r em i x e dt o g e t h e ro n eb yo n e1 nt n e 1 a b t h e ya r ec r e a t e db yu s i n gv i r t u a le x p e r i m e n t o nt h ec o m p u t e r l t u s u a l l yr e q u i r e st oc o m p a r ea n da n a l y s i sal a r g ea m o u n to fd a t a l nt h e c h e m i c a ld a t a b a s e ，t h e r e f o r e ，h o wt os e a r c ht h ei n f o r m a t i o nt h ed e v e l o p e r c o n c e r n e di nt h el a r g ec h e m i c a ld a t a b a s er a p i d l y ， a n ds h o wi t i ng r a p h b e c o m eah o tr e s e a r c ha r e a i nt h i sp a p e r ，b a s e do nt h eo t h e rp e o p l e sw o r k ， ar a p i da l g o r l t h m f o rv i s u a l i z i n g1 a r g ec h e m i c a ld a t a b a s ei na1 0 w d i m e n s i o n a ls p a c e ( 2 d o r3 d )i sp r e s e n t e da saf i r s ts t e pi nt h ed a t a b a s ea n a l y s i sa n dd e s l g n a p p l i c a t i o n s t h ea l g o r i t h mi sb a s e do nt h es i n g u l a rv a l u ed e c o m p o s i t i o n ( s v d )c o m b i n e dw i t ham i n i m i z a t i o np r o c e d u r eb yu s i n gt r u n c a t e d n e w t o n m e t h o d i tm a i n t a i n st h eo r i g i n a lr e l a t i o n s h i pb e t w e e ne a c hp a i ro ft h e c o m p o u n d sv e r yw e l l ， a n di nt h ea l g o r i t h mt h es p a r s em a t r i xc o n c e p t1 s u s e di no r d e rt or e d u c et h em e m o r ya n dc o m p u t i n gt i m ec o s t i l h ec o m p u t e r v is u a l i z a t i o nt e c h n i q u ei sa l s oi n t r o d u c e dt ot h ep a c k a g et og l v e t h e d e v e l o p e ran i c eg u i f i r s t ，t h ep a p e re s t a b l i s h e st h em a t h e m a t i c a lm o d e l ，a n dt h e ni tw i l l s h o wt h ed e t a i lo ft h ew h 0 1 ea l g o r i t h m ， h c 岍i tp r o g r a m e da n d t h e e x d e r i m e n t a lr e s u lt s i nt h ee n d ，i td e s c r i b e st h eg u io ft h ep a c k a g e t h ee x d e r i m e n t a lr e s u l t sw i l ls h o wt h ee f f i c i e n c yo ft h ea l g o r i t h m k e y l r o r d ：s y n t h e t i cd r u g ， s v d ， t r u n c a t e d n e w t o n ， s p a r s e衄t r i x ， v i s u a li z a t i o n 浙江大学硕士学位论文 第一章绪论 1 1 问题的产生背景及其应用 近几十年来，随着生物技术的兴起，生物信息学成为一门计算机科学和 数学应用于分子生物学而形成的新兴的交叉学科，它涉及生物学、数学、计 算机科学和工程学，依赖于计算机科学、工程学和应用数学的基础，也依赖 于生物实验和衍生数据的大量储存“1 。生物信息学不只是一门为了建立、更 新生物数据库及获取生物数据而联合使用多项计算机科学技术的应用性学 科，也不仅仅是只限于生物信息学这一概念的理论性学科。事实上，它是一 门理论概念与实践应用并重的学科，随着生物技术的进步，它正在越来越显 示出它的重要性。在该领域中，由生物学家和计算机科学家共同研究生物分 子信息的获取、管理、分析和利用。这些生物分子数据具有丰富的内涵，其 背后隐藏着人类目前尚不知道的生物学知识。充分利用这些数据，通过数据 分析、处理，揭示这些数据的内涵，得到对人类有用的信息，是生物信息学 创建的目的。生物信息学接要凭借这些数据，以计算机、网络为工具，用数 学和信息科学的理论、方法和技术去研究生物大分子，发现生物分子信息组 织的规律。该学科一个重要应用就是用于生物制药的药物合成领域。 现今，药物合成技术已经被大大小小的药物发明公司广泛采用。药物合 成的基本概念就是，建立巨大的分子“库”，不是象传统做法那样逐个的合 成化合物；然后对这个库进行高小的筛选，鉴别出最有希望的“先导”制药 化合物。这个概念象大潮一样席卷了制药工业和生物技术工业，几乎没有一 家公司对此无动于衷。现在，几乎每一个从事制药研究开发的机构都在用或 者宣称自己在用某种组合技术，许多有创新精神的公司已赶忙向制药行业提 供化合物库、高效的分类筛选方法、合成方法、硬件和软件、分析工具等。1 。 浙江大学硕士学位论文 然而，至今还没一个通用的高效的筛选平台得到广泛的应用。 但是，随着分子生物学技术的不断进步和基因组研究的不断深入，以及 生物数据的本身一些特点，例如，数据量巨大，其中既有生物分子序列的信 息，又有结构和功能的信息；既有生命本质信息，又有生命表象信息，并且 数据之间存在着密切的联系，因此，这样的分子库随着人类产生与获取数据 的能力增长而成数量级地增大。例如，常常存在如下的情况，发现某几种成 分的组合对某种疾病有很好的疗效，但是这些成分中有特定的成分对人体有 害或者有副作用，因此需要在分子库中找到一种新成分来替代该有害成分 成，即如何在海量的化学数据库中快速地搜索到专该种有害成分的结构和属 性相同的化学物质并可视化地展现出来。面对这浩如烟海的数据，想通过人 工地分析这些数据进而得到各个数据之间的关系几乎是不可能的了。为了最 大限度地减少冗余，提高筛选效率，众多的科学家们都开发新的算法帮助设 计出针对具体受体需要和可口服药物的物理化学特性的化合物库。同时，由 于计算机速度与存储器容量的持续不断提高，这也为各学科采用数据模拟进 行研究开辟了通途，因此，一些公司也正在开发达到类似目的的计算机程序。 由于图象的在分析数据方面的直观性等优点，数据可视化技术也逐渐被 应用到这个领域中来。 数据可视化技术凭借计算机强大的处理能力和计算机图形图像学的基 本算法以及可视化算法把巨大数量的数据转换为静态或动态图象或者图形 呈现在人们面前，并允许通过交互的手段控制数据的抽取或画面的显示，使 隐含在数据之中不可见的现象成为可见，为人们分析数据，理解数据形成概 念，找出规律提供了强有力的手段。在生物分子合成领域，可视化技术可使 研究者通过直观交互的方法增添或者删除某类分子或者分子个数来控制合 成物质的性质从而缩短设计新物质的周期。高技术研究需要实验，但是实验 往往需要花费大量的人力物力和财力，越是复杂，高端的实验越是开销巨大， 而且很多时候还受到很多客观因素的制约，所以大量的实验用计算机来通过 浙江大学硕士学位论文 数据模拟进行是个颇受欢迎的方法。少量的、必要的实际实验只在为校正数 值模拟使用参数和检验数据摸的方法而进行。这正是众多制药公司所期望 的。 本文正是在这样一个背景下，提出了一种新的数据分析算法，并结合可 视化技术，为制药研究者提供了一种较为高效，直观的在分子数据库中快速 搜索查看分子数据集的方法。 1 2 研究现状及理论背景 在许多分子数据库分析的应用中，研究数据集中各个数据之间的距离关 系一直是很重要的一个研究方向。这样的研究将极大地方便我们对已有的数 据进行聚类分析，例如，如果我们将高维数据降维到或者映射到低维的向量 空间，尤其是在二维或三维的情况下，我们甚至可以直接通过对图象的空间 距离的观察来对数据进行聚类分析。通常，这样的映射问题就等价于这样一 个新问题： 即在二维或三维空间中找到n 个点，要使得这n 个点，点与点 之间的距离尽可能的和它们在m 维空间中距离相匹配。但是这样的问题的 解通常不是唯一的，因此问题的关键便成了如何在众多解中寻求一个最优化 的解，这往往取决于我们采用什么样的距离误差函数和使用了什么方法来使 该误差取得它的最小值。 这样的距离问题在众多场合有着广泛的应用，例如在分子结构学，生物 制药业等，但是这个问题是非常难求解的。通常我们只能得到一个局部近似 解，并且还需要一个很好的初始值用于迭代求解。 最近，这种距离几何化近似的方法被应用于分析和二维映射高维分子数 据库的研究中。一般都是s 锄o n 法来求解。在s 猢o n 法中采用了最速下降 法( s t e e p e s td e s c e n tm i n i i z a t i o na l g o r i t l l i ) 和一个随机的初始值。 这个算法有两大缺点：一是采用这种算法收敛的速度很慢，二是当问题的维 数增大是得到的结果往往和原始的距离匹配程度较差。因此一个好的初始值 浙江大学硕士学位论文 对整个问题的求解很重要，但是事实上又是比较难取得的。k o h o n e n 又提出 了一种采用神经网络来进行自我映射的方法( s e l f o r g n i z i n gm a pm e t h o d ) ， 这种方法也有它自身的缺陷。p r o f x i e 也提出了一种t n 法，他采用了截 断式牛顿法，并引入了采用有限差分法计算得到的二阶偏导阵。 现有的算法降维后的核心算法一般都是以下几种方法： 1 使用最速下降法( s t e e p e s td e s c e n tm i n i i z a t i o na l g o r i t h m ) 最速下降法( s t e e p e s td e s c e n tm i n i m i z a t i o na l g o r i t h m ) 是最简单 的梯度下降法。该方法对于下降方向的选择总是选择使函数值f 下降速度快 的方向，根据数学知识，我们知道这个方向是就函数f 的梯度方向的反方向。 搜索起始点则是随机选取任意一个点“o ，然后顺着我们找到的梯度反方向 下降，直到我们和所需要的解足够的接近，它的迭代公式为： + l = 雌一k v 厂( 吆) = t k g ( t ) 其中g 心j 为在任意一点以上的梯度。 由此可见，最速下降法是非常简单，便于实现的一种快速迭代法。而且 该方法也很稳定，只要最小值存在，该方法总能在有限的迭代次数中逼近它。 但是，即使它有着上述的一些明显的优点，它还是有一个非常严重的缺点： 它的收敛速度非常缓慢。经常开始于一个非常不错的收敛速度，但是在收敛 过程中，下降速度越来越慢，甚至几乎停滞，如下图所示： 浙江大学硕士学位论文 一一 叭一 x 一 ， 但是如果给定一个很好的迭代起始点，即使整个系统的不是良性的，它 还是可以取得一个较好的收敛速度。由此可见，该算法的优劣很大程度上依 赖于迭代起始点的选择。因此，如果对整个系统最小值的所处的区域有大致 的判断，则最速下降法是一种行之有效的不错的收敛方法。然而在实际应用 中，事实上，我们很难知道其全局最小值的位置，所以最速下降法通常必须 要与其他的优化方法合并使用才能达到较好的效果。不然很难保证其收敛速 度和是否能收敛到全局最小而不是一个局部最小。 2 共轭梯度算法( c o n j u g a t eg r a d i e n tm e t h o d ) 在上面最速下降方向( 负梯度方向) 上调节迭代步长。虽然这是性能函 数减小最快的方向，但并不能产生最快的收敛。究其收敛速度不快的原因， 我们可以从上图发现每次迭代它都要转一个直角的弯进行下去。针对最速下 降法该缺点的原因，对搜索方向寻找进行了改进，使它不在是简单地沿着梯 度的反方向进行，而是沿着共轭梯度方向进行，通常会产生比最速下降方向 快得多的收敛。 共轭( c o n j u g a t e ) 的意思是指两个不相等的向量： d 一和日，它们于任 意对称正定矩阵正交。例如：q 为一个对称正定矩阵，则有： 浙江大学硕士学位论文 q 嘭= 。 采用共轭这个概念是为了让我们的搜索方向吐不依赖于其它所有的搜 寻最小值的方向嘭。 把这个概念引入共轭梯度算法，我们就可以通过梯度向量来构造我们的 共轭集，以保证每次新的搜索方向是和以前的搜索方向正交的。它的步骤通 常如下： 共轭梯度算法在第一次迭代前选择初始搜索方向为最速下降方向( 负梯 度方向) p 0 2 一g q 下一个搜索方向要保证与以前的所有搜索方向共轭( 新的搜索向量与以 前的搜索向量集合正交) ，通常把新的最速下降方向和以前的搜索方向组合产 生： 见= 一舐+ 反n l| i = o ，1 ，2 ， 其中常量屏常有多种不同的计算方法，因此就有若干个不同版本的共 轭梯度算法，以下是两个常用的计算展的方法： 屈：粤趣 ( 1 ) 繇一i 一 屏= 粤 ( 2 ) 颤一1 暑 迭代步长k 通常也由下式给出： 九。：毒坠 ( q 矾) ，但当q 未知或者计算太复杂的时候，步长k 也可以通 过线搜索方法( l i n e s e a r c h ) 获得 浙江大学硕士学位论文 有了搜索方向和迭代步长后，我们就可以根据迭代公式： 进行迭代，寻求系统的最小值。 从以上可以看出，共轭梯度算法比最速下降法有了较大的改进，能处理 一些大型收敛问题，但是在迭代的初始常会产生一些无意义的，或者不高效 的搜索方向，而且随着迭代的进行，搜索方向将会逐渐失去它的共轭性 ( d i r e c t i o nc o n j u g a c y ) 。 3 生物遗传算法 遗传算法是模拟生物在自然环境中的遗传和进化过程而形成的一种自适 应全局优化概率搜索算法。它最早由美国密西根大学的h o l l a n d 教授提出， 起源干6 0 年代对自然和人工自适应系统的研究。7 0 年代d ej o n g 基于遗传 算法的思想在计算机上进行了大量的纯数值函数优化计算实验。 与传统的优化算法相比，遗传算法主要有下述几个特点：以决策变量的 编码作为运算对象；直接以目标函数值作为搜索信息；同时使用多个搜索点 搜索信息，这就是通常所说的遗传算法特有的隐并行性；使用概率搜索技术。 1 3 新的算法 在本文中，我们采用了一种新的算法来实现求解。该算法由两部分组成： 第一部分是采用奇异值分解( s i n g u l a rv a l u ed e c o m p o s i t i o n ，s v d ) 定 理来求得一个低维的映射值。奇异值分解定理。我们将会在后面部分做具体 介绍，这是一个在图象处理、编码和解码中，被广泛应用来进行数据压缩的 方法。 这个分解定理的优点是，只需要给定一个输入数组( 高维) ，并不需要映 射初始值，并且它的计算复杂度为o ( n 幸n 蜘) ，空间复杂度为o ( n 书m ) 。我们 从中发现s v d 映射的精度取决于奇异值的大小：若前两个奇异值的数值要远 浙扛大学硕士学位论文 大于其他的奇异值，那么进行二维映射就能取得较高的精度。这个规律可以 推广到三维甚至更高维，即：如果前n 个奇异值的数值要远大于其他的奇异 值，那么映射到n 维以下的维度就能取得较高的精度。 第二部分是对初始值的改进。在这里，我们把由s v d 产生的初始值作为截 断式牛顿迭代法( t r u n c a t e dn e w t o nm i n i m i z a t i o ni t e r a t i v em e t h o d ) 的 迭代初始值进行迭代一直到达到我们所设定的精度或者达到最大的迭代次 数，并采用稀疏矩阵精确求解得到二阶偏导阵( h e s s i a n 矩阵) 。截断式牛顿 迭代法通过求解牛顿方程来产生搜索方向，有较快的收敛和稳定性等优点。 相对于其他算法的采用梯度和用有限差分法所求得的二阶偏导阵来说，会在 性能上有较大提高，并且考虑在在算法中计算二阶偏导阵是个开销很大的过 程，因此我们采用了不完全压缩的二阶偏导阵，从而减小了开销，从结果看 出仍然可以取得较满意的结果。 第三部分是整个软件的用户图形接口，使用了o p e n g l 编程技术在l i n u x 环境下为用户的操作和使用提供了一个良好的交互平台，给研究带了极大的 便利。尤其是在有了蛋白质数据仓库( p r o t e i nd a t ab a n k ，p d b ) 以后， 可以很方便的直接读取p d b 数据文件重新分子图像，直观的展现在用户面前， 而不是以前那样只出来一堆数据，缺乏可操作性。 根据以上思想，我们编写了一个测试程序用以测试和使用，事实证明它 可以很高效地求得距离误差函数的最优解。 浙江大学硕士学位论文 第二章问题的数学建模及算法的数学描述 2 1 数学模型的建立 假设一个数据库s 中含有n 个化合物( 药成分) ，每个化合物又被m 个属 性( d e s c r i p t o r ) 所描述。因此，我们可以定义数据集s 如下： s = ，五，瓦) ， 其中向量z = ( t 。，五：，硝，) 7 表示数据集s 中的第i 个化合物 ( c o m p o u n d ) ，其中的实数 ) 则是表示与各个属性相关的数据。 我们也可把数据集s 写成矩阵的形式，用行表示化合物的个数( n ) ，用 列表示各项属性( m ) ，可以得到如下的矩阵： 州墨如，科：f o l 矗l 矗。j 对于大型的数据集，这个矩阵具有典型的n m 的结构，n 甚至可以上万。 为了衡量各个化合物之间的相似性，我们需要定义一个量化距离来测量 它们之间的差异，通常我们都采用如下几种方法来定义距离“： 1 ) 明氏( m i n k o w s k i ) 距离。 d 。( 口) ：f 杰f x 。一x ，。1 4 “9 口l 当9 2 l 时，为绝对距离； 当9 2 2 时，为欧氏距离； 当g 。3 时，为切比雪夫距离。 当各变量的测量值相差悬殊时，采用明氏距离并不合理，需要先对数据 浙江大学硕士学位论文 标准化，然后用标准化后的数据计算距离。 明氏距离特别是其中的欧氏距离是人们较为熟悉的，也是使用最多的距 离。但明氏距离存在不足之处，主要表现在两个方面：第一，它与各指标的 量纲有关；第二，它没有考虑指标之间的相关性，欧氏距离也不例外。 2 ) 马氏距离 设表示指标的协差阵即： 却 。，其中旷击萎( 驴圳矿i ) “= l ，卯 暑= 去喜z 。，弓= 去言x 。 如果。1 存在，则两个样品之间的马氏距离为 d ；( m ) = ( x ，一_ ) 。( 工，一) 这里o 为样品xr 的p 个指标组成的向量，即原始资料阵的第f 行向量。 样品z ，类似。 顺便给出样品z 到总体g 的马氏距离定义为 d 2 ( x ，g ) = ( z 一) - 1 ( x 一) 其中为总体的均值向量，为协方差阵。 马氏距离既排除了各指标之间相关性的干扰，而且还不受各指标量纲的 影响。除此之外，它还有一些优点，如可以证明，将原数据作一线性交换后， 马氏距离仍不变等等。 3 ) 兰氏距离 郴，= 古善制“乩，聆 此距离仅使用于一切z f 0 的情况，这个距离有助于克服各指标之间量 纲的影响，但没有考虑指标之间的相关性。 浙江大学硕士学位论文 计算任何两个样品x 。与x ，之间的距离d f ，其值越小表示两个样品接近 程度越大，吼值越大表示两个样品接近程度越小。如果把任何两个样品的距 离都算出来后，可排成距离阵d ： d = d ，吐2 ，矾。 d 2 1 ，d 2 2 ，- - d 2 n d 。l ，d 。2 ，- - d 肼 其中讲t2 d ：z 2 “2 0 。d 是一个实对称阵，所以只须计算上三角 形部分或下三角形部分即可。根据d 可对胛个点进行分类，距离近的点归为 一类，距离远的点归为不同的类。 在这里我们采用最基本的明氏( m i n k o w s k i ) 距离：取g 2 2 的情况 一 6 f = 慨一一岷( 靠一嘞) 2 ， ( 2 1 ) yi t l 其中il ff 表示e u c l i d e a n 范数。因此对于一个有n 个化台物的数据 集来说，如果只考虑i 0且盯，+ l = = a 。= o 。 若= 。，) 7 ，则根据( 2 7 ) ，我们可以把置表示成m 维空间 1 4 浙江大学硕士学位论文 r 1 的正交基向量 唯 ：。的线形组合： 置= 6 。“。唯= a 。咋， i = l ，2 ，n( 2 8 ) - 1- l 因为我们知道a 。一一a 。= o 。因此置可以重写为一个新的向量： 置= ( a i ，0 2 “f 2 ，d ，o ，0 ) 7 。 ( 2 9 ) 基于上式，我们就可以定义一个由r ”空间产生的低维子空间r “，记 其中的向量为r ，则： i = ( a l “m a 2 巩2 ，c 。“如。) 7 ， k = 0 ，1 ，2 ，( 2 1 0 ) 通过这样的映射就可以得到我们所需要的低维向量，但它未必能满足我 们需要的误差精度，这就需要我们采用截断式牛顿迭代法来提高精度，使我 们的误差函数e 取得最小值。 这个分解定理一个很明显的优点就是，只需要给定一个输入数组( 高 维) ，并不需要映射初始值，并且它的计算复杂度为o ( n 女n m ) ，空间复杂 度为0 ( n 女m ) 。并且我们发现s v d 映射的精度取决于奇异值的大小：若前两 个奇异值的数值要远大于其他的奇异值，那么进行二维映射就能取得较高的 精度。这个规律可以推广到三维甚至更高维，即：如果前n 个奇异值的数值 要远大于其他的奇异值，那么映射到n 维以下的维度就能取得较高的精度。 2 2 2 截断式牛顿迭代法 截断式牛顿迭代法是一种通过反复迭代来提高精度的方法，提到迭代法， 就必然会考虑到两个问题：一个是是否能达到全局最值，另一个是收敛速度 是否快。 由于截断式牛顿迭代法是一种非线性的全局优化法，因此可以避免落入 局部最值，而且它采用的是二阶偏导数矩阵，h e s s i a n 阵，其定义如下： 浙江大学硕士学位论文 日( ，( 而，而，矗) ) = a 2 ，a 2 ，a 2 ，a 2 厂 研瓤鸥铂驰阮瓯 a 2 ，a 2 厂a 2 ， a 2 厂 氓粥钙眠 ；ii； a 2 魄铆 a 2 ， 瓯 a 2 ， 瓯如 a 2 厂 锑 ( 2 1 1 ) 因此比起同类的迭代法，如采用一阶偏导阵的最速度下降法( s t e e p e s t d e s c e n ta l g o r i t h m ) 显然要收敛速度快很多。但是二阶导数阵往往会随着数 据集的增大而迅速增大o ( n 2 ) ，因此如何减少计算它所需要的时间开销和存 储开销也是一个很重要的问题。这部分内容将在我们算法的实现部分介绍。 截断式牛顿迭代法主要包括两个大循环：外循环和内循环。 外循环主要实现对点集( r ) 的修正，是主循环，其迭代表达式为： 】，“= y + 九t p ， k = o ，l ，2 ， ( 2 1 2 ) 其中：p 和都空间月”“中，因此其矩阵的大小都为n l 。w ，为 下降方向，k 为下降步长，y 0 为迭代初始值。 在本算法中由又内循环迭代寻得，这里采用的是带预处理的共轭梯度 法，通过求解牛顿方程肿p = 一g ，得到我们的下降方向，其中h 为误差函数e 的h e s s i a n 阵，g 为误差函数的梯度。 下降步长k 则是由直线搜索( 1 i n e s e a r c hm e t h o d ) 求的。 2 3 算法终止的条件 当迭代所得的p 满足下式的时候： 浙江大学硕士学位论文 0 9 ( y ) 8 e 。( 1 + f e ( ，) i ) ， ( 2 1 3 ) 其中：g ( 矿) 为误差函数e 在p 处的梯度向量，。是一个很小的正数( 我 们采用1 0 5 ) ，我们就认为所得的低维空间的点集p 能很好的近似并保持原 有的距离关系，并终止牛顿迭代。这是因为我们知道最小值的点处的梯度为 o ，所以梯度越接近o ，说明达到的效果越理想。 浙江大学硕士学位论文 第三章算法的实现 3 1 程序的数据结构框架 我们的整个程序数据结构是定义于t e m p l a t en u m e r i c a lt 0 0 1 k i t ( 简 称t n t ) 的结构之上。t e m p l a t en u m e r i c a lt 0 0 1 k i t 是由美国m a t h e m a t i c a l a n dc o m p u t a t i o n a ls c i e n c e sd i v i s i o n ，n a t i o n a li n s t i t u t eo fs t a n d a r d s a n dt e c h n o l o g y 用c + + 开发的一个用于科学计算的抽象类库，它为用户提供 了一些标准的在科学数值计算中常用的基本数据结构，如多维矩阵 ( a r r a y l d ，a r r a y 2 d ，a r r a y 3 d ) ，向量阵( v e c t o r ) 等，以及这些基本数据结构 之间的运算，经过自己的改写可以对矩阵运算进行简化，使用方便如 m a t i 。a b 。 浙江大学硕士学位论文 3 2 稀疏矩阵的存储结构和计算算法 在我们的测试程序中采用的是由3 个一维数组或者( 链表) 来实现一个 稀疏矩阵的存储，考虑的h e s s i a n 是对称的，我们首先可就可以节省一半的 存储空间，只存储它下三角的非零元素。我们如下定义3 个一维数组： 数组m a ：存储h e s si a n 所有非零元素的值。 数组m j ：存储h e s s i a n 所有非零元素的列下标，存储的顺序于m a 中的 值的顺序相对应。 数组m i ：存储h e s s i a n 矩阵中每一行中第一个非另元素在数组m a 中的 位置。 通过这3 个数组我们就可以找到任意一个非零元素在原h e s s i a n 中的位 置。 1 9 浙江大学硕士学位论文 例如：如下一个对称的稀疏矩阵： m = l2 2 4 05 3o o o o 0 o 3 5 o 6 o o8 7 0 0 9 o o o o 7 o 0 9 10 0 2 若采用新的存储方式只存储它的下三角中的非零元素，即： m = 1o 24 o 5 3 0 o o 0 0 0 o 0 0 6 o 0 8 7 0 o 9 0 o 0 0 o 0 o o 10 o 2 并以行的顺序来存储，我们可以得到如下三个一维数组： 删= 1 2 456387l9 2 】； 脚= 【1 l2 2314 3 5 4 6 】； m = 【l 2 4 6 8l o 1 2 】。 矩阵肘第一行的第一个非零元素是1 ，它存储在m 的第一个位置 m a ( 1 ) ，所以m i ( 1 ) = 1 ；而矩阵m 第2 行的第一个非零元素是2 ，它存储在 m a 的第二个位置m a ( 2 ) ，故m i ( 2 ) = 2 ；而第3 行的第一个非零元素5 存储在 在m a 的姒( 4 ) ，因此有m i ( 3 ) = 4 ，以此类推。 从上我们可以总结出这种存储方式的几个性质： 1 数组m a 和m j 的长度等于整个稀疏矩阵中非零元素的个数：而数组m i 的长度等于该稀疏矩阵的维数加一。 2 埘( f + 1 ) 一埘( o = 第i 行非零元素的个数，既第第i 行非零元素存储在 m a 的第m i ( i ) 到m i ( i + 1 ) 一1 的位置；数组m i 最后一个元素埘0 + 1 ) 一l = 浙江大学硕士学位论文 该稀疏矩阵中含有的非零元素的总数。 当稀疏矩阵的维数和含有的零元素增大时，我们可以看到这样的存储方 式可以有效的减少存储空间。 接下来来我们来看看这种稀疏矩阵的存储方式在计算仍然有其简便性。 我们以在该算法中常出现的正方形对称矩阵乘一个列向量的运算为例 子，即运算q = h d 其中h 为一个n 术n 的对阵稀疏矩阵，并以上述方式存储； d 是一个n l 的列向量。 由于对阵稀疏矩阵h 已经表示为3 个一维数组m a ，m j ， f i ，因此若我 们还按照h 原来的行坐标和列坐标( i ，j ) 来计算：吼= 嘭，显然我们 j _ 1 需要花费很多时间在寻找每一个元素巩在新的存储方式下处于数组m a 的 位置下标，最差的情况便是我们遍历了h 矩阵第i 行所有非零元素存储在 m j 中的列坐标后发现没有和我们所寻找的元素峨相对应的列坐标时，我们 才发现它是零元素，根本无须做运算。因此，我们考虑到h 矩阵的对称性和 这种存储方式自身的优点，我们不难发现对于每一个非对角线上的元素由于 = 的缘故，所以存储在m a 数组中的每一个元素峨( f ，) 其实对q 向 量中吼，g ，都有贡献， g = t ，= l q i = h 矿d i = h4 d j ij = l 在遍历所有非零元素数组m a 的时候，每个原非对角元素的值e ，分别 嘭，吐做乘法运算，并分别累加到吼，吼，如此只需遍历m a 的所有元素一 遍即可求得向量p ，具体算法如下： 浙江大学硕士学位论文 定义双精度型向量p ，并将向量q 置o v e c t o r q = o ： 遍历数组m a f o r ( i = o ；i n ：i + + ) ( f o r ( k = m i i ：k m i i + 1 ：k + 十) f 计算g f = 吼+ t ，的列下标j 存储数组m j 中 q i = q i + m a k ：# d m j k ： i f ( m j k ! = i ) 若不是对角线上的元素，则计算它对g ，的贡献 q i = qj + hu d i q m j k = q m j k + m a k 州 i ： 这只是计算简化中的一例，使用该种存储方式还有其它地方可以大大简化计 算量，这对大型数据集的计算是很有益处的。 3 3 误差函数e 梯度的算法实现 从前面的数学建模分析中，我们得到了我们的距离误差函数 月一jn1n l ” e ( 巧，e ，) = 言嘞( d ( 巧，) 2 6 ；) 2 = 去k ( ( 1i t l j = j + 1 r l _ l ，= j + 1t = l ( 3 1 ) 其中： 嚣；i ：，而q 是有一个很小的正数。 4 y 俺l ， 、 l i 浙江大学硕士学位论文 由于 i ，蔓，) 都是向量，所以我们把误差函数写成向量形式更容易 使我们求得我们的误差函数的梯度和二阶偏导数阵。 如w 令勺( j ，) = 饥一) 2 6 ；，v 。= ，m ：，】r ，则误差函数e i = l 可以简化成： 刚弘，= 丢琵啸 ( 3 2 ) 其中： 勺= ( v ，一v ，) 7 ( v 。一v ，) 一6 ；，同前。 由此，我们对e 做求导运算， g ：v e ：( 婺，罢，婺，- ，兰， 明1哕i 枷吼lb 。 = 去，f ，v ( 吩2 ) _ r f _ 1j = 件l 1h ln = 去嘞o v 吩 ( 3 3 ) 其中v 勺2 ，罢，兰) t 饥，饥伽7 孙他彬爱 浙江大学硕士学位论文 令g f ：丢v o ： o o 0 ( v ，一v ，) o o 一0 ( v 一v ) 0 0 ，从式子我们可以发现每个只有第i 个子矩阵和第j 个子矩阵非零，其它均为零，而且这两个子矩阵在数值上只 有一个符号差异，这样就给我们的计算带来了极大便利，每对( i ，j ) 只要 计算1 个子矩阵即可，那是一个l o w $ l 的矩阵。 因此，误差函数e 的梯度如下： h ln g = v e = 岛。 ( 3 4 ) ，= 1 = f + l 其对应的算法程序如下： v e c t o r g v a l u e ( a r r a y 2 d y ，v e c t o r d i s t ，i n tn ，i n t m d i ) v e c t o r g ( n 棚d i ) ： i n tu = 0 ，u u ，v v ： d o u b l et e m p o ，t e i i l p l ，t e m p 2 ，d e r j ，s u m ，r i j ： f o r ( i n ti = 0 ：i n l ：i + + ) f o r ( i n tj = i + 1 ：j n ：j + + ) s u m = 0 0 ： f o r ( i n tk = 0 ；k 卫上 鹾一i 毛。t 瓠t 瓠l ( 3 9 ) 通过满足上述两式，并对h e s s i a n 矩阵限制在一个带状内，其余均为零。 如下图所示： 箕尝旦。 到陲零螽岛? 囊。 l i 器i 。 。 扩， 【i ：0 袅誉 宽 度 占 一 宽度6 这样我们新的h e s s i a n 矩阵就成了稀疏矩阵，再进行运算就可以大大缩 小计算量和存储所需要的空间。并且带的宽度对整体精度并没有太大影响， 这在我们后面的实验结果中可以看出。 浙江大学硕士学位论文 3 5 核心算法流程图 3 5 1 外循环 浙江大学硕士学位论文 3 5 2 内循环 浙江大学硕士学位论文 第四章实验结果及分析 我们的实验是运行于r e d h a t1 i n u x9 o 下。 4 1 随机初始值和采用s v d 产生初始值对迭代收敛的影响 下图是在一个包含1 6 0 0 化合物，每个具有9 个属性的数据集上进行的 实验，我们分别采用了s v d 法生成初始值和随机产生初始值这两种方法来观 察s v d 法对整个迭带过程的影响。在这里，由于随机产生的初始值可能搜寻 不到正确的下降方向从而导致跌代过早终止，我们只是随机打乱了s v d 产生 的初始值的顺序，下图是实验结果： 图4 1 随机初始值和s v d 产生初始值对迭代性能的影响 从上图我们可以看出，s v d 法产生的初始值下降的比随机打乱顺序的情 浙江大学硕士学位论文 况下降的快，并且曲线更稳定，震荡少，而且到9 0 次跌代以后迅速下降趋 势，而采用随机初始值的情况下，梯度没有什么下降的趋势了。 在下表的比对中我们可以看出更显著的效果，下表来自p r o f x i e ( u n i v e r s i t yo fw i s c o n s i na tm i l w a u k e e ) ，他首先采用使用s v d 的方法 产生低维初始值： 表4 1 随机初始值和s v d 产生初始值对迭代性能的影响 从上表可以更明显的看出s v d 产生初始迭代值的优势，在总的迭代次数 和计算时间上都有明显的改进。 4 2 不同目标函数对迭代收敛的影响 在前面我们简单的介绍过我们采用新构造的误