矩阵分解的常用方法

2024-09-18

矩阵分解的常用方法（通用7篇）

矩阵分解的常用方法篇1

0 引言

随着国民经济的快速发展,电力系统迅速向多级化、大电网、交直流联合输电及大区联网运行发展,电力系统的动态分析即使在离线状态下也变得十分困难。实际上对一个大电力系统的动态研究最感兴趣的只是其中某一个区域,称为研究系统,而距此区域较远的区域,研究中只考虑其对研究区域的影响,可作降阶及简化,这种拟作简化的区域称为外部系统。保留研究系统不变,而对外部系统在保证其对研究系统的动态影响不畸变的条件下,进行简化的过程称为动态等值[1]。

同调等值法是动态等值的基本方法之一,其关键步骤为识别同调机群。同调机群的识别方法有很多,如相关系数法[2]、状态空间法[3]、模糊聚类法[4,5,6]、模糊相关自组织数据分析算法[7,8]等。对于大系统而言, 这些方法一般比较复杂, 计算量较大,难以满足工程实用化的要求。当前,从实时的曲线对发电机进行分群[9,10]的方法已经比较成熟,但是该类方法对于功角曲线的数据时间窗宽度有一定要求,且数据处理方法本身存在一定的局限性,给分群结果带来较大的影响。随着新理论的不断涌现,文献[11]提出轨迹特征根的概念,文献[12,13]将此思想进一步推广,并在小故障下验证了该思路的实用性。在此基础上,文献[14]用故障后系统的轨迹特征根进行机组分群趋势的预测,但该方法能否适应复杂系统和连锁故障场景及多摆失稳的情况,有待进一步研究。文献[15]分别提出了基于小波变换和人工神经网络的同调机群识别法,但这些方法仍克服不了计算复杂、耗时的缺点,且神经网络建模技术所强调的训练误差最小化的做法,易引起模型对样本数据的过拟合,从而导致模型的泛化能力较差。

文献[16]提出了一种基于主成分分析(principal component analysis,PCA)方法的同调机群识别方法。该方法通过构造量测样本数据的线性组合,重新组成一组新的相互无关的综合变量,最终获得量测数据的主成分分量。从本质上说,该方法是通过降维达到分群目的,原理简单、计算量小,可以满足大系统在线应用的要求。然而,PCA完全是从数学角度分析量测样本矩阵,因此只能抽取样本的代数特征,从而有可能导致在分解结果中产生负值,而负值往往无法从物理角度进行解释,这样就失去了对原始问题的最佳解释,也就是说,提取特征的过程丢失了部分信息,有可能会导致最终结果的不准确。

针对上述问题,非负矩阵分解(non-negative matrix factorization,NMF)提供了一种新的解决方案[17]。NMF是由Lee和Seung在1999年提出的一种新的矩阵分解算法[17],与PCA等其他特征提取算法相比,该算法对分解结果非负的限制使得NMF具有实现上的简便性、分解形式的可解释性等优点。为简单方便地对分解结果进行聚类分析,采用了K均值聚类算法[18]。K均值聚类算法是一种能用在许多聚类问题上的快速迭代算法,算法简单且聚类效果较好,是到目前为止应用比较成熟的一种聚类分析算法。

基于上述分析,本文提出一种基于NMF的电力系统同调机群识别方法。通过NMF对发电机转速数据降维,降维结果归一化后再利用K均值聚类对其进行聚类,以达到同调机组的分群目的。

1 NMF和K均值聚类算法

1.1 NMF问题描述

NMF问题可描述为[19]:对一个n×m阶的非负矩阵V,可以将其分解为一个n×r阶的非负矩阵W和一个r×m阶的非负矩阵H的乘积。

分解后可理解为,原矩阵V中的列向量是对矩阵W中的所有列向量的加权和,而权重系数是矩阵H对应列向量中的元素。故称W为基矩阵,H为系数矩阵(权矩阵)。这种基于基向量组合的表示形式具有很直观的语义解释,它反映了人类思维中“局部构成整体”的概念[20]。一般情况下,r的选择值要比n小,即r满足(n+m)r<mn,这时用系数矩阵代替原始数据矩阵,就可以对原始数据矩阵降维(且是非线性的维数约减),得到数据特征的降维矩阵,从而可以减少存储空间,节约计算资源[21]。

1.2 NMF实现算法

NMF的实现可以表述成最优化问题,常用的目标函数有2个,其中,i和j分别表示矩阵的行和列。

1)矩阵V与矩阵M(M=WH)的欧氏距离的平方,即

$∥ V - Μ ∥^{2} = \sum_{i j} (V_{i j} - Μ_{i j})^{2} (1)$

式中:Vij和Mij分别为矩阵V和M中的元素。

2)矩阵V与矩阵M推广的K-L离散度,即

$D (V ∥ Μ) = \sum_{i j} (V_{i j} \lg \frac{V_{i j}}{Μ_{i j}} - V_{i j} + Μ_{i j}) (2)$

式(1)与式(2)都是当且仅当V=M时达到最小值0。

因此,优化问题转变成:在约束W≥0,H≥0下,寻找使得式(1)或者式(2)表示的目标函数达到最小值时所对应的矩阵W和矩阵H。文献[17]提出了一种乘法迭代算法,即从任意非负初始值出发,交替更新矩阵W和矩阵H,直到它们的改变足够小。具体算法步骤如下。

步骤1:对非负矩阵W和H分别随机赋初值。

步骤2:更新W和H。

对于式(1),更新法则为:

$W_{i k} \leftarrow W_{i k} \frac{A_{i k}}{B_{i k}} (3) Η_{k j} \leftarrow Η_{k j} \frac{C_{k j}}{D_{k j}} (4)$

式中:Wik,Hkj,Aik,Bik,Ckj,Dkj分别为矩阵W,H,A,B,C,D的元素,A=VHT,B=WHHT,C=WTV,D=WTWH;k为降维后矩阵的维数。

对于式(2),更新法则为:

$\begin{array}{l} W_{i k} \leftarrow W_{i k} \frac{\sum_{j} \frac{Η_{k j} V_{i j}}{Μ_{i j}}}{\sum_{j} Η_{k j}} (5) \\ Η_{k j} \leftarrow Η_{k j} \frac{\sum_{i} \frac{W_{i k} V_{i j}}{Μ_{i j}}}{\sum_{i} W_{i k}} (6) \end{array}$

步骤3:重复步骤2直至收敛。

本文以式(1)作为NMF的目标函数。

1.3 K均值聚类算法

聚类是将数据集划分为若干个类(cluster)的过程,使同一个聚类中数据对象(样本点)具有较高的相似度,而不同聚类中的对象相似度较低。聚类方法中的K均值聚类算法具有复杂度较低、效率高、可扩展等优点。其基本思想是首先确定若干初始聚类中心,然后逐步改变或者调整这些中心,使聚类趋于合理[22]。

K均值聚类算法流程如下。设X=[x1,x2,…,xi,…,xN]为原始数据,k为聚类数。首先,从N个数据对象中随机选取k个对象作为初始聚类中心,其他对象则根据它们与这些聚类中心的相似度(距离)分别分配到最相似的类中。假设cj为第j个类的聚类中心,则xi与cj的相似度为:

$\begin{array}{l} d (x_{i}, c_{j}) = [(x_{i 1} - c_{j 1})^{2} + \dots + \\ (x_{i k^{'}} - c_{j k^{'}})^{2} + \dots + (x_{i n^{'}} - c_{j n^{'}})^{2}]^{\frac{1}{2}} (7) \end{array}$

式中:xik′和cjk′分别为xi与cj的第k′个属性。

然后计算每个更新的类的聚类中心,假设第j个类中的样本为[xj1,xj2,…,xjm′],即包含m′个样本,则该类的聚类中心为cj=[cj1,…,cjk′,…,cjn′],其中cjk′可根据式(8)求得:

$c_{j k^{'}} = \frac{x_{j_{1} k^{'}} + x_{j_{2} k^{'}} + \dots + x_{j_{m^{'}} k^{'}}}{m^{'}} (8)$

不断重复上述过程,直到标准测度函数收敛为止(从表现形式上看即更新后的聚类中心与更新前一致),一般采用均方差作为标准测度函数,其形式为:

$J = \sqrt{\frac{\sum_{j = 1}^{k} \sum_{l = 1}^{m^{'}} ∥ x_{j_{l}} - c_{j} ∥^{2}}{Ν - 1}} (9)$

最后得到的聚类具有如下的特点:聚类内部尽可能紧凑,不同类之间尽可能分开。

2 基于NMF的同调机群识别

在发电机的分群问题上,将具有相似动态行为的机组划分为一组机群,需要处理的信息是发电机的转子角等指标在研究时间内的相互接近程度[4]。为了便于与基于PCA的同调机群识别方法的分群结果相对比并清晰展示NMF的思想,分群指标选取文献[16]中采用的发电机实际转速。

发电机转速矩阵作为原矩阵,其自身是一个非负矩阵,而NMF的乘法迭代算法保证了基矩阵W和系数矩阵H中无负数元素,满足了一定物理意义的约束。从另一个角度来看,测取的发电机转速可以看做是表征发电机动态行为的物理量,参照NMF在文本领域应用中的解释,可将分解得到的基矩阵中的列向量理解为对发电机转速产生影响的各个因素的集合,系数矩阵H则描述的是每个因素集合的重要程度。将列向量通过系数矩阵中的权值进行线性组合得到原矩阵,这就是本文发电机转速进行NMF的直观的“局部构成整体”的感受。通过对系数矩阵H的聚类,也就是对因素重要程度的聚类,可以实现同调机群的识别。

文献[16]采用PCA对发电机转速降维,然后分别取前2个主成分和前3个主成分进行机组分群,结果表明用后者进行分群比前者更精确。本文对NMF和PCA进行比较时,将降维的维数定为3,这样有利于作三维图,使结果可视化,确定机组的同调特性,并进行分群结果对比。

设系统有m台发电机,利用仿真测取各发电机N个转速数据,记为

$V (n) = [v_{1} (n), v_{2} (n), \dots, v_{m} (n)]^{Τ} n = 1, 2, \dots, Ν (10)$

接着对N×r阶的非负矩阵W和r×m阶的非负矩阵H随机赋非负初值,其中r=3。依照式(3)和式(4)所示的更新法则更新W和H,重复上述步骤直至式(1)所示的目标函数达到最小值。降维所得到的r×m阶的矩阵H即为进行同调机群识别所需要的数据。可以看出,NMF解决了源数据维数高的问题,实现了对源数据的特征提取。

对于降维矩阵H,先进行归一化,再用K均值聚类算法对它进行聚类。通过对提取特征的聚类来确定各机组的同调特性,实现机组的分群。

基于NMF的同调机群识别的流程如图1所示。

3 算例分析

本文采用的算例以New England 10机39节点系统作为分析对象。通过算例中不同线路端点发生故障的2个方式来验证基于NMF的同调机群识别方法的有效性,并通过K均值聚类的标准测度函数J的值来比较基于NMF和基于PCA的方法的分群效果。在设置的2种不同故障方式下用PSD-BPA分别进行时域仿真,测取各发电机300个转速数据作为源数据进行算例分析。

3.1方式1

方式1在Bus9-Bus8的Bus9侧发生三相瞬时接地短路故障(定义为Bus9*-Bus8),0.10 s切除故障。下面分别采用基于NMF和基于PCA的同调机群识别方法,对10机39节点系统的机组进行分群。

将方式1中的源数据分别用NMF和PCA进行降维。NMF实现算法的目标函数值为7.996 4×10-9。NMF实现降维的过程实际上是一个优化问题,目标函数值越小,基矩阵W和系数矩阵H的乘积更接近于原矩阵V,W和H保留V的信息更多,也体现了“局部构成整体”的感受。降维数据归一化后利用K均值聚类的聚类结果如图2和图3所示。聚类结果均分为5类,即(G30),(G39),(G31,G32),(G34,G38),(G33,G35,G36,G37),由此可以看出,基于NMF和基于PCA的同调机群识别方法的结果是一致的。分群结果如表1所示。

此故障方式下各机3 s内的原始功角摇摆曲线如图4所示, G30为参考发电机。可以看出,(G31,G32),(G34,G38),(G33,G35,G36,G37)的同调特性较强,此方式下将10机系统分为5群符合系统发电机的功角轨迹趋势,同调机群结果为(G30),(G39),(G31,G32),(G34,G38),(G33,G35,G36,G37),验证了基于NMF和基于PCA的同调机群识别方法在该方式下的分群是合理的。

方式1分别用NMF和PCA对源数据降维并归一化后,K均值聚类的标准测度函数的值如表2所示。NMF的标准测度函数的值为0.24,PCA的标准测度函数的值为0.36。可见,方式1中NMF之后利用K均值聚类的标准测度函数的值小于PCA之后利用K均值聚类的标准测度函数的值。标准测度函数所描述的是每个数据样本与所属类的聚类中心之间的均方差,因此NMF使得同调机群的聚类内部更紧凑,表明此方式下相比于PCA,NMF可使同调机群聚类的效果更好。

3.2方式2

方式2在Bus28-Bus26的Bus28侧发生三相瞬时接地短路故障(定义为Bus28*-Bus26),0.10 s切除故障。将方式2中的源数据分别用NMF和PCA进行降维。NMF实现算法的目标函数值为4.539 4×10-10。降维数据归一化后利用K均值聚类的聚类结果如图5和图6所示,基于NMF和基于PCA的同调机群识别方法的结果是一致的。分群结果如表3所示。

此故障方式下各机3 s内的原始功角摇摆曲线如图7所示,G30为参考发电机。可以看出,(G31,G32,G33,G35,G36)的同调特性较强,此方式下将10机系统分为6群符合系统发电机的功角轨迹趋势,验证了基于NMF和基于PCA的同调机群识别方法在该方式下的分群是合理的。

方式2分别用NMF和PCA对源数据降维并归一化后,K均值聚类的标准测度函数的值如表4所示。2种方法的标准测度函数值的比较同样表明此方式下NMF比PCA的同调机群聚类效果更好。

综上所述,由方式1和方式2这2种不同故障的仿真算例可以看出,NMF较好地保留了原矩阵的主要信息,基于NMF进行同调机群识别也是可行有效的。标准测度函数值的比较表明,基于NMF的同调机群识别方法使同调机群的聚类内部更紧凑,即聚类效果更好。

4 结语

本文提出了基于NMF的同调机群识别方法,利用NMF对发电机转速源数据进行降维,得到具有非负性质的低维映射矩阵,归一化后进行K均值聚类,确定机组的同调特性,实现机组分群。基于不同方式下的算例测试结果验证了本文算法的有效性。本文的模型和方法具有以下特点。

1)低维矩阵较好地保留了原矩阵的信息,简化了分群工作。

2)在某些数据样本中,负值往往无法从物理角度进行解释,而NMF对分解结果非负的限制至少满足了一定的物理意义的约束,并使其具有分解形式的可解释性。

3)分群结果准确,聚类效果更好。通过与基于PCA方法的同调机群识别方法的分群结果的比较,清晰展示了本文方法的思想,并且通过标准测度函数值的比较表明了本文方法在聚类效果上的优势。

摘要：为了解决源数据维数较大的问题,提出了一种基于非负矩阵分解(NMF)的同调机群识别方法。采用发电机角速度作为源数据,使用NMF算法对其进行降维。由于此低维矩阵具有非负性质,因而该模型在消除冗余数据、降低维数的同时,保留了原始问题的实际意义。对低维矩阵归一化,再利用K均值聚类算法对其进行聚类,达到同调机群的分群目的。通过New England 10机39节点系统比较了基于NMF和主成分分析方法的分群效果,验证了基于NMF的同调机群识别方法的有效性。

关键词：同调机群,分群,非负矩阵分解,K均值聚类

矩阵分解的常用方法篇2

关键词：社交网络,社团结构,非负矩阵分解

1 概述

社交网络如微博、微信、博客等在人们的生活中发挥的越来越广泛的应用,这些社交网络也在网络舆情等方面发挥了重要作用。这是因为在社交网络中,用户之间相互连接、影响,从而促进了信息像洪水一样的快速而广泛的传播,与传统媒体、门户网站不同之处就在于,他没有一个明显的传播主渠道。在社交网络中,用户经常只是和少部分其他用户信息交互频繁,而与其他大部分用户联系很少,这样,在用户之间就形成了许多明显的圈子,即社团结构。社团内的用户之间相互联系密切,相互共享信息或者进行合作,如微博、facebook、Twitter等网络应用中,有共同兴趣的节点相互分享视频、评论等信息,形成一种社团结构[1]。在这些圈子内,总会存在诸如名人、超级用户(如微博中的大V)和活跃用户,是社团中的活跃分子。相比之下,用户与社团外部之间的联系就相对稀疏的许多,但是,许多舆情也通过这些稀疏的连接将信息从一个社团散布到其他社团,从而引起其他社团广泛的关注。因此,在社交网络上,既要挖掘这些社团结构,获得用户之间的关系结构,便于对网络进行管理,又要掌握网络中社团活动的重要用户,掌握信息广泛传播的关键节点,只有这样,才可以全面掌握社交网络的活动规律,正确预见潜在的网络事件,挖掘诸如犯罪关系网等重要信息,保障网络的健康有序。

社交网络是一种复杂网络,它规模庞大,结构复杂,表现在节点数目巨大,网络结构呈现多种不同特征。连接具有多样性,节点之间的连接权重存在诧异。节点集可能属于非线性动力学系统,例如节点状态随时间发生复杂变化。自2002年Gir⁃ven和Newman提出社区挖掘的概念至今,新的理论、方法层出不穷,新的应用领域不断涌现.它不仅吸引了来自计算机科学、物理、数学、生物和社会学等领域的研究者,成为最具挑战性的多学科交叉研究课题之一;而且已被应用于社会网络分析,如恐怖组织识别、组织结构管理、网络推荐等方方面面。采用非负矩阵分解(NMF)的方法进行社团的发现。

2 基于非负矩阵分解的社团发现的方法

非负矩阵分解是近年来出现的一种新方法。它对矩阵分解中的元素都添加了非负的限制,即大于或等于0,这与事物的整体是部分之和的直观概念相符合,在应用中具有很强的可解释性,但是,对于非对称矩阵的非负分解,目前的算法效果并不太好,而在实际中,很多网络连接关系是有向的,它们对应的矩阵就是非对称的。我们通过对目标矩阵进行简单变换,添加必要的约束项,来实现基于非负矩阵分解的有向加权网络的社团发现,同时进一步分析用户在社团中的重要性,判断用户在社团中所处的作用和特点,从而可以达到有效分析和解剖网络结构的目的。

2.1 社团发现模型

对于一个网络连接关系,一般使用图来描述。其中:表示第i个节点,表示第i个节点到第j个节点关系的大小。其中,X是社团指示矩阵,其元素表示第i个节点属于第j个社团的大小。S是社团关系矩阵,其元素表示第i个社团与第j个社团的密切程度,n表示节点的个数,k表示社团的个数,一般来讲,社团的个数要远远小于节点的个数。

可见,对公式进行转置运算,并不改变节点属于哪个社团的特性,只是将社团之间的关系结构进行了融合。于是,对于非对称的目标矩阵,可以转化为对称矩阵。如下所示:

令,则:

由此,我们可以构建基于矩阵非负分解的社团发现目标函数:

|‖F表示Frobenius范数,简称F范数,‖‖1是l1范数,用来控制模型的稀疏度,让节点尽可能被清晰的划分。是可调节的平衡参数,0表示规模为n×k,元素全为0的矩阵,同样I表示元素全为1的矩阵。E是单位矩阵。

2.2 模型的优化求解

定义矩阵Nk×k是一个元素全部为1的矩阵,公式(3)可以改写为:

根据文献,我们可以获得迭代规则如下

其中,N表示元素全为1的k×k维矩阵。

3 实验

3.1 Zachary karate数据集

Zachary karate数据集是一个真实网。它是Zachary空手道俱乐部成员关系网络,是社会学分析等领域中最常用的一个小型检测网络之一。从1970到1972年,Wayne Zachary用三年时间观察了美国一所大学空手道俱乐部成员间的社会关系,并构造出了社会关系网(Zachary’s karate club network)。网络中的每个节点分别表示某一个俱乐部成员,节点间的连接表示两个成员经常一起出现在俱乐部活动(如空手道训练、俱乐部聚会等)之外的其他场合,即在俱乐部之外他们可以被称为朋友。调查过程中,该俱乐部因为主管John A(.节点34)与教练Mr.Hi(节点1)之间的争执而分裂成2个各自为核心的小俱乐部。这个矩阵经常被用来检测算法的有效性。

3.2 其他真实数据集

我们还利用其他真实数据集来做实验比较,如:American College Football,Neural Network数据集,以模块度值[11]作为衡量指标,实验表明,我们的算法都得到了较精确的社团划分结果。除了Dolphin Social Network数据集外,我们的算法都要好于其他基于非负矩阵分解的算法。

4 结论

本文针对网络中节点复杂的连接关系,改进了以往社团划分算法只能解决无向网络问题的缺点,构建了基于有向加权图的非负矩阵社团划分模型,提出了模型求解的算法,并且提出了社团中节点的重要性等性质的定量分析,在真实的网络数据集上的实验表明,我们提出的算法正确的给出了社团划分结果,并对节点在社团中的性质进行了详细的分析,其结果符合人们的直观认识。我们的算法对于研究社交网络中复杂的人际关系提供了一种有用的方法,对于网络的结构分析,信息在用户之间的传播有参考意义。

参考文献

[1]曹贵强.纸币识别及红外鉴伪技术的研究[D].辽宁科技大学,2014.

矩阵的三角分解及其应用篇3

1. 定义

如果方阵A可分解成一个下三角矩阵L和一个上三角矩阵U的乘积, 则称A可作三角分解或LU分解。如果方阵A可分解成A=LDU (1.1) , 其中L是单位下三角矩阵, D是对角矩阵, U是单位上三角矩阵, 则称A可作LDU分解。

2. A的LDU分解和LU分解求法

(1) A的LDU分解求法:取A的一阶主子式, 作LDU分解A1= (1) (a11) (1) , L1= (1) , D1= (a11) , U1= (1) , 用 (1) 式- (4) 式确定v1, l1, 从而L2=βl1L110β, D2=β0D1 d20β, U2=β0V1 1v1β, A2=L2D2U2。然后重复使用 (1) - (4) 式得到A的顺序主子式的LDU分解。Ak=LkDkUk, k=1, 2, …, n, 当k=n时, 即完成A的LDU分解。

(2) A的LU分解求法:先按照 (1) 的步骤求出A的LDU分解, 再令DU的乘积为U’, 即求出A的LU分解。

3. LU分解的推广和改进

矩阵A的LU与LDU分解都需要假设A的前n-1阶顺序主子式非零。如果这个条件不满足, 可以给A左或右乘以置换矩阵, 就把A的行或列的次序重新排列, 使这个条件满足, 从而就有如下的带行交换的矩阵分解定理。

定理1设A是n阶非奇异矩阵, 则存在置换矩阵P使PA的n个顺序主子式非零。该定理的证明可在矩阵论教科书中查到。

4. LU分解在解方程组中的应用

设A是n阶非奇异矩阵, 则存在置换矩阵P使PA=LDU=LU’ (1.3) 。

其中L是单位下三角矩阵, D是对角矩阵, U是单位上三角矩阵, U’是上三角矩阵。

如果方程组AX=b (1.3) , 系数矩阵A是n阶非奇异矩阵且△k≠0, (k=1, 2, …, n-1) , 则存在三角分解A=LU.于是得到与原方程组同解的以三角矩阵为系数矩阵的方程组βUx=yLy=b (2) (1) , 方程组的 (1) 式解出代入 (2) 式解出, 这就是线性方程组的三角分解法。如果方程组的系数矩阵A的某个顺序主子式△k=0, (k=1, 2, …, n-1) , 可用定理1的推论考虑与其同解的方程组PAX=Pb (1.4) , 即也可用三角分解法求解。

AX=b的解。

解:A的LU分解为

则Ly=b的形式为

于是可求AX=b解向量为

二、矩阵的QR分解

定义如果非奇异矩阵A能够化成正交矩阵Q与非奇异上三角矩阵R的乘积, 即A=QR, 则称A可QR分解。

定理2如果A是n阶非奇异矩阵, 则存在使A=QR成立 (2.1) , 且除去相差一个对角元素的绝对值 (模) 全等于1的对角矩阵因子外, 分解是唯一的。

证记A的列向量α1, α2, …, αn方阵A是非奇异矩阵, 所以列向量α1, α2, …, αn线性无关, 对A的列向量α1, α2, …, αn, 正交化得正交向量β1, β2, …, βn, (α1, α2, …, αn) = (β1, β2, …, βn) K,

于是R=diag (β1, β2, …, βn) K, 则Q是正交矩阵, R是非奇异上三角矩阵, 且有A=QR。

关于矩阵分解问题的讨论篇4

关键词：矩阵分解,正交阵,对角化,逆矩阵,对称阵

一、矩阵分解之和

1. 对称阵和反对称阵之和

定理1:任何一个n阶矩阵A, 均可表为一个对称阵与一个反对称阵之和, 即A=B+C, 其中B′=B, C′=C。

定理2: (1) 定理1表示法是唯一的;

(2) B与C的乘积是不可以交换的, 即BC≠CB, 并且不是对称阵也不是反对称阵;

(3) BC=CB的充要条件是A′A=AA′, BC=-CB的充要条件是A′A=AA′。

2. 纯量阵与迹为零的矩阵之和

定理3:任一n阶方阵均可表为一个纯量阵与一个迹为零的方阵之和。

分析:假若A=X+Y, 而A= (aij) , 令X为纯量阵, 即X=xIn, Y= (yij) 是迹为零的方阵。即。

现在的问题是求x, yij, i, j=1, …, n。其中i≠j时, yij=aij显而易见, 所以求yii (i=1, …, n。) 与x即可。由A=X+Y, 得Tr (A) =Tr (X) +Tr (Y) , 然而Tr (X) =nx, Tr (Y) =0, 故。又由A=X+Y, 得aij=x+yii, 于是。

3. 幂零阵与对角阵之和

定理4:复数域上任一个n阶矩阵A可分解成A=B+C, 其中B为幂零阵, C相似于对角阵, 且BC=CB。

4. 其它矩阵之和

定理5:任一秩为r的矩阵均可表为r个秩为1的矩阵的和。

二、矩阵分解之积

1. 满秩分解

定理6: (满秩分解定理) :任一秩为r>0的m×n矩阵A均可表为m×r列满秩阵B与r×n行满秩阵C之乘积即A=BC。

2. 幂等阵的分解

定理7:设A为秩r的n阶幂等阵, 则:

(1) 存在n阶可逆矩阵P, 使P-1AP=。

(2) 存在两个n阶实对称矩阵B、C, 使A=BC, 即任一幂等阵必可分解为两个实对称阵的乘积。

3. 任一方阵分为两个对称阵之积

定理8:设A为复数域上任一n阶方阵, 试证:

(1) A=T-1AT其中T为n阶对称的可逆矩阵。

(2) A=ST其中S与T皆为n阶对称阵, 并且它们之中至少有一个是可逆的。

4. 三角分解

定理9:设A= (aij) 是n阶方阵, 它的顺序主子式全不为零, 则存在:非奇异下三角形矩阵B与非奇异上三角形矩阵C, 使A=BC。

5. QR分解

定理10:设A为任意n阶实阵, 则必有分解式A=QR, 其中Q为n阶正交阵, R为主对角元皆为非负数的n阶上 (下) 三角形阵。

推论1:设A为n阶实非奇异矩阵, 则A可以分解成A=QR。其中Q为正交阵, R是一个对角线上全为正实数的上三角矩阵, 并且这种分解是唯一的。

推论2:正交阵A的QR分解可分解为有限个镜象乘积的形式。

6. 正定阵的分解

定理11:设A为n阶正定矩阵, 则:

(1) A=T′T, 其中T为n阶实可逆阵;

(2) A=B2, 其中B为n阶正定矩阵;

(3) A=R′R, 其中R为主对角元为正的上三角形n阶实阵:A=T′T= (QR) ′ (QR) =R′R。

推论3:n阶正定阵A, 对于任意给定的正整数m, 则存在唯一的正定实对称矩阵B, Bm=A。

7. 可逆阵的分解

定理12:设A的n阶实可逆阵, 则A必可表为一个正交阵Q与一个正定阵S之积, 即A=QS, 且此种表示是唯一的。

定理13:设A为n阶可逆阵, 且A的特征根互异, 则在复数域上有n阶B, 使A=B2。

定理14:设A为n阶非奇异实矩阵, A可分解为:A=Pdiag (λ1, …, λn) T。

其中P, T皆为n阶正交阵, λi>0, i=1, 2, …, n, 且λ12, …, λn2是矩阵A′A的全部特征根。

8. 奇异值分解

定理14:设A为秩是r的m×n实矩阵, 则A有分解式A=其中P, T分别是m阶与n阶正交阵, ai>0, i=1, 2, …, r。

9. 其它几种形式的矩阵分解

定理15:任何一个n阶方阵均可表为In+aijEij这种形状的矩阵乘积, 其中皆为1×n矩阵) , 通常称为矩阵单位, aij为数, i, j=1, 2, …, n。

关于矩阵分解问题分为两种情形, 矩阵分解之和与矩阵分解之积。在理论与应用上都是十分重要的。在矩阵分解之和的讨论中, 我们应掌握几种特殊矩阵, 如对称阵、反对称阵、幂等阵、幂零阵与纯量阵, 以及它们的性质, 并用构造法解决问题。在证明中要用到一个重要的分解式, 其中P, Q为可逆阵, r为A的秩, 这种分解式的应用十分广泛, 我们可以结合初等变换, 求矩阵的秩, 逆矩阵等应用方面起重要作用, 也是矩阵构造的基础。

矩阵分解之积中主要包括三角分解、QR分解、奇异值分解、满秩分解及可逆阵的分解等。若一矩阵A可分解为两个矩阵的乘积, 或可以相似对角化, 则可以利用分解式计算A的幂次, 或A的多项式。我们也可以利用分解式计算矩阵的秩。若矩阵A可分解为两个或几个矩阵的和 (如分解为单位阵的和) , 我们则可以用分解式计算与A可交换的矩阵, A的幂次, 或A的特征值。求解方程组的难易程度与系数矩阵A的结构特征有关。如果A是分块三角矩阵或通过合同变换化成分块三角矩阵, 那么求解方程组是较方便的。对一般的非奇异矩阵, 我们通常可化为下三角矩阵和上三角矩阵乘积, 然后求逆, 可利用满秩分解求广义逆。矩阵的奇异值分解问题在最优化问题, 特征值问题, 最小二乘方问题, 广义逆矩阵问题等方面都涉及广泛。

参考文献

[1]杨尚骏, 蔡茜.关于矩阵分解为稳定矩阵之和的讨论.大学数学, 2003, 19, (3) :60-62.

[2]时贞军.二次规划的矩阵分解算法.曲师范大学学报, 1991, 17, (1) :54-57.

[3]陈焕艮.正则环上矩阵分解.数学杂志, 1999, 19 (4) :405-407.

[4]钱雅利.矩阵分解与矩阵行列式估计的证明.黑龙江教育学院学报, 1992, 15, (1) :91-92.

[5]黄刚石.一个受限非负矩阵分解方法.东南大学学报, 2004, 4, (2) :189-193.

矩阵分解理论及应用篇5

关键词：等价分解,三角分解,奇异值分解,Fitting分解

在代数学中,矩阵分解就是把矩阵分解成某种规范型[1,2,3]。在实际应用中,利用矩阵的某些分解来解决一些实际的工程数学问题有明显的效果。如计算某些特大、特殊矩阵时,矩阵的三角分解非常有作用,可以大大简化很多计算过程[4]。矩阵的QR分解在状态估计具有很好的计算效率[5]。Schur分解作为一种强有力的工具,在处理矩阵等式和矩阵不等式[3]的过程中有着非常重要的作用。奇异值分解是研究数学的一种重要方法,并在最优化问题、特征值问题、广义逆矩阵计算、谱估计、控制理论等领域[3,6,7,8],有极其重要的作用。矩阵的Fitting分解可看作是复矩阵的Jordan分解在一般域上的推广,它在运筹与控制论方面有至关重要的作用。

由于许多文献都是研究矩阵分解的某一特性,因此本文将系统地归纳和总结矩阵分解理论及其一些理论应用。本文把矩阵分解大致分为等价分解、三角分解、谱分解、奇异值分解和Fitting分解等五类。特别地,给出了矩阵等价分解的一种新的证明。在理论应用方面,利用等价分解求解矩阵方程AXA=A;利用矩阵三角分解求解线性方程组;展示了谱分解在极值原理中的应用;利用矩阵的奇异值分解求解Moore-Penrose逆、解决了最小二乘问题并且利用该分解给出了最小二乘问题的通解;利用Fitting分解证明了Drazin逆的存在和唯一性。

1 矩阵的等价分解

1.1 矩阵的等价分解

定理1.1.1(等价分解)若,则存在m阶的可逆阵P及n阶可逆阵Q使得

其中r=rank(A)

证明:设是的基,将其扩充成

定理1.1.4(满秩分解)若,则存在列满秩阵和行满秩阵使得,其中r=rank(A)。

1.2 矩阵等价分解的应用

矩阵的等价分解在矩阵秩、矩阵方程等许多问题的研究中有广泛的应用。下面举例说明:

例1.2.1设

解:由矩阵等价分解知存在m阶可逆阵p及n阶可逆阵Q满足

因此矩阵方程的全部解为

其中XÁ,XÂ,XÃ为任意的适当矩阵。

2 矩阵的三角分解

2.1 矩阵的三角分解

2.1.1 LU分解

三角分解法是将方阵分解成一个上三角阵和一个下三角阵,这样的分解法又称为LU分解法。

定理2.1.1(LU分解)对一任意方阵A,均可以分解为两个三角阵的乘积,即:PA-LU,其中P为置换阵,L为下三角阵,U为上三角阵。

定义2.1.2(Doolittle分解)若方阵A有分解:A=LU,其中L为单位下三角阵,U为上三角阵,则称该分解为Doolittle分解。

定理2.1.3设A缀Rnxn,若A的前n-1阶顺序主子式不为0,则Doolittle分解可以实现,即A=LU,其中L为单位下三角阵,U为上三角阵,并且分解式唯一。

证明:对阶数n用数学归纳法来证。

则L为单位下三角阵,U为上三角阵,故A=LU。

下面证明唯一性:因为A的顺序主子式,若A可逆,则,其中为单位下三角阵且可逆,为上三角阵且可逆,故有,上式左边为上三角阵,右边为单位下三角阵,从而

若A不可逆,则其中L,L1为单位下三角阵且可逆,U,U1为上三角阵不可逆且秩为n-1,于是有为单位下三角阵,为使上式成立,则必须有,否则上式矛盾。从而,进而有

定义2.1.4(Crout分解)若方阵A有分解:A=LU,其中L为下三角阵,U为单位上三角阵,则称该分解为Crout分解。

2.1.2 对称阵的Cholesky分解

定理2.1.5设A为对称阵,则存在唯一分解:A=LLT,其中L为单位下三角阵。

证明:由Doolittle分解,A有唯一分解:A=LU,则

定理2.1.6设A为对称正定阵,则存在唯一分解:A=LDLT,其中L为单位下三角阵,D为对角线元素大于零的对角阵。

2.1.3 三角分解的应用

例2.1.7试用三角分解求解方程组

解:由Doolittle分解有:

2.2矩阵的正交三角分解

2.2.1 QR分解

定理2.2.1若,则存在酉矩阵Q及上三角阵R使得。

定理2.2.2设,则A可唯一分解为,其中,U的各列标准正交,R是r阶正线上三角阵。

推论2.2.3若A为可逆阵时,则A可唯一分解为

其中U,U1是酉矩阵,R是正线上三角阵,R1是正线下三角阵。

推论2.2.4若,则A可唯一分解成,其中L是r阶正线下三角阵,,U的各行标准正交。

推论2.2.5若则A可以分解为,其中,,且U1的各列标准正交,U2的各行标准正交,R1是r阶正线上三角阵,L2是r阶正线下三角阵。

证明:根据定理1.1.4和定理2.2.2和推论2.2.4可证。

2.2.2 Schur分解

定理2.2.6若,则存在酉矩阵U使得,其中T为对角线上的元素都是A的特征值的上三角阵。

证明:用数学归纳法。当A的阶数为1时,定理显然成立。现在假设A的阶数为时,定理成立,下面考虑A的阶数为k时的情况。

下面我们利用Fitting分解给出矩阵Drazin逆的存在性与唯一性的证明。

3矩阵的谱分解

矩阵的谱分解在极值问题的研究和有关特征值的不等式方面有着重要的作用。

3.1矩阵的谱分解

定理3.1.3为正定阵的充分必要条件是存在n个线性无关的向量,使得

定理3.1.4(谱分解)设A为n阶单纯阵

证明:由单纯阵及特征值特征向量的定义可证。

3.2矩阵的极分解

定理3.2.1 A是正定阵的充分必要条件是存在可逆阵B,使得证明:由正定阵的Cholesky分解易证。ÁA=B B

定理3.2.2 A是正定阵的充分必要条件是存在正定阵B,使得,k为任意正整数。

证明:由推论3.1.2可证。

定理3.2.3(极分解)若n阶可逆阵,则存在一个正定Hermite阵和一个酉矩阵U,使得

证明:由定理3.2.1和定理3.2.2易证。

3.3矩阵谱分解的应用

矩阵谱分解在极值原理中的应用:

4矩阵的奇异值分解

奇异值分解是另一种正交矩阵分解法;也是最可靠的分解法。

4.1矩阵的奇异值分解(SVD分解)Á

定义4.1.1对任意A,设的非零特征值为,则称

显然

定理4.1.2(SVD分解)设,有奇异值,则存在m阶酉矩阵U和n阶酉矩阵V使得,

4.2矩阵奇异值分解的应用

4.2.1用奇异值分解求解Moore-Penrose逆ÁÂÁÁÂÁ

定义4.2.1设满足下面的方程,则称X为A的Moore-Penrose逆,简称M-P逆,常记为。

定理4.2.2对任意存在且唯一。

证明:若

若,则A有奇异值分解,其中为酉矩阵。令,容易验证X满足(i)-(iv)条,这就证明了A+的存在性。下证A+的唯一性。

4.2.2用奇异值分解来讨论最小二乘问题

4.2.3

另外

因而。于是推出是Ax=b的最小二乘解。从而线性方程组Ax=b的所有最小二乘解可表示为

及解中的最小长度解。

5矩阵的Fitting分解

矩阵的Fitting分解将矩阵分解为可逆部分和幂零部分的直和。

5.1矩阵的Fitting分解

定义5.1.1对n阶矩阵A,若存在正整数满足:AK=0,则称A为幂零阵。

定理5.1.3(Fitting分解)设则存在可逆阵T使得

其中D为可逆阵,N为幂零阵。

证明:对n用数学归纳法。当n-1时,显然。假设A的阶数小于n时结论成立。下面证明当A的阶数是n时结论仍然成立。

若A=0或A可逆,结论显然成立。

若。令

由归纳假设,存在可逆阵P1使得其中D为可逆阵,N1为幂零阵,于是

则由N1幂零及引理5.1.2知N幂零。

5.2矩阵Fitting分解的应用

定义5.2.1设

称X是A的一个Drazin逆,简称D-逆。

定理5.2.2设,则A的Drazin逆存在且唯一。ÁÁA CÁÎ

证明:若A可逆,易见A-1是A的唯一的Drazin逆。若A不可逆,,设A的Fitting分解为

则容易验证是A的Drazin逆,并且在A的Fitting分解给定的情况下,A的Drazin逆是确定的。因此,为了展示A的Drazin逆是唯一的,只需证:对于A的任意两个Fitting分解

6结论

本文归纳和总结了矩阵分解的类型及其相关的理论应用,把矩阵分解大致分为等价分解、三角分解、谱分解、奇异值分解和Fitting分解等五类,在证明过程中注重理论的前后衔接。在很多实际工程问题中,矩阵分解都具有很大的优势,但因作者水平有限,本文仅仅只是略提了一下有关矩阵分解在许多实际工程方面的应用,以后可以做更深一步的研究。

参考文献

[1]罗家洪.矩阵分析引论[M].广州:华南理工大学出版社,1998.

[2]苏育才,姜翠波,张跃辉.矩阵理论[M].北京:科学出版社,2006.

[3]张显,仲光苹,高翔宇.系统与控制中的矩阵理论[M].哈尔滨:黑龙江大学出版社,2011.

[4]王岩,王爱青.矩阵分解的应用[J].青岛建筑工程学院学报,2005,26(2):90-93.

[5]Joanne A.Foster,John G.,Martin R.Davies and Jonathon A.Cham-bers.An Algorithm for Calculating the QR and Singular Value Decom-positions of Polynomial Matrices[J].IEEE Transactions On Signal Processing,March,2010,58(3):1263-1274.

[6]Soo-Chang Pei,Ja-Han Chang,Jian-Jiun Ding,and Ming-Yang Chen.Eigenvalues and Singular Value Decompositions of Reduced Bi-quaternion Matrices[J].IEEE Transactions On Circuits and Systems-I:Regular Papers,2008,55(9):2673-2685

[7]VIRGINIA C.KLEMA AND ALAN J.LAUB.The Singular Value Decomposition:Its Computation and Some Applications[J].IEEE Transactions On Automatic Control,1980,25(2):164-176.

基于矩阵分解的物资管理系统优化篇6

物资计划是物资管理的重要环节,科学的物资计划可以降低物资管理的成本。物资管理系统不仅需要管理大量不同种类的物资,还需要同时满足所有需求单位的需求,因而运营难度很大。每个需求单位对每种物资的需求量往往难以确定,很多时候,物资管理系统并不清楚每个需求单位对每个种类物资的真实需求量,它只知道部分需求单位对部分物资的需求量。鉴于此,需要根据已知的需求单位对物资的需求量推断出未知的需求单位对物资的需求量。我们通常会利用推荐系统中常用的协同过滤模型来预测未知的需求单位对物资的需求量,从而准确估算所有需求单位对所有物资的需求量,进而科学地对物资进行计划,达到优化物资管理系统的目的。

2 协同过滤模型

协同过滤模型(Collaborative Filtering)是经典的个性化推荐模型(Recommend System),它可以根据用户的历史行为向用户推荐合适的商品。推荐系统和物资管理系统有一定的相似之处,因此,我们将协同过滤模型应用在物资管理系统的优化中。协同过滤模型能根据需求单位的历史行为数据对需求单位对某种物资的需求状况进行预测。协同过滤策略通过分析已有的需求单位和物资之间的关联模式,来预测新的需求单位和物资之间的关联。如果需求单位A和需求单位B在历史上对某种类型的物资的需求量是相似的,那么协同过滤模型就认为需求单位A和需求单位B的需求相似,从而认为需求单位A和需求单位B对其他种类的物资的需求量也是相似的。基于矩阵分解的协同过滤模型使用潜在因素/特征来描述需求单位和物资的特征,这些特征往往是难以明确地进行定义的。我们将在下一节详细介绍基于矩阵分解的协同过滤模型。

3 基于矩阵分解的协同过滤模型及其学习算法

3.1 基于矩阵分解的协同过滤模型

基于矩阵分解(Matrix Factorization)的协同过滤模型是使用最广泛的协同过滤模型之一。对于基于矩阵分解的协同过滤模型,每个需求单位和物资都有一个向量(Vector),用于表示它的潜在特征,这些潜在特征(Latent Features)一般是不能具象化的。基于矩阵的协同过滤模型,通过分析历史上需求单位对物资的需求量的模式来学习这些潜在特征。如果某个需求单位的潜在特征和某种物资的潜在特征相关性非常高,那么协同过滤模型就会认为这个需求单位对这种物资的需求量很高。基于矩阵分解的协同过滤模型对于现实中的许多场景都是十分有效的。

基于矩阵分解的协同过滤模型根据部分需求单位对部分物资的需求量来构建一个需求单位-物资的矩阵,矩阵的其中一个维度是需求单位,另一个维度是物资。举个例子,假设有M个需求单位、N种物资,那么模型所构建的需求单位-物资的矩阵X的维度就是M×N,即X∈RM×N。矩阵X中的元素表示第i个需求单位对第j种物资的需求量,我们称这个值为观测值(Observed Value)。如果第i个需求单位对第j种物资的需求量是未知的,那么xij就是缺失的,我们称这个值为缺失值(Missing Value)。

假设潜在特征的数量是D,基于矩阵分解的协同过滤模型将需求单位-物资矩阵X∈RM×N分解成两个矩阵U∈RM×D和V∈RN×D:

其中,矩阵U和矩阵V分别表示需求单位和物资的潜在特征,矩阵U中的一行ui表示第i个需求单位的潜在特征,矩阵V中的一行vj表示第j种物资的潜在特征。ui中的每一维表示第i个需求单位在每一个潜在特征中的权重,同理vj中的每一个维表示第j种物资在每一个潜在特征中的权重。向量ui和向量vj的内积(Inner Product)表示第i个需求单位和第j种物资的相关度,或者说是第i个需求单位对第j种物资感兴趣的程度。潜在特征的数量D往往远小于需求单位的数量M和物资的数量N。

基于矩阵分解的协同过滤模型关键在于设法求出矩阵U和矩阵V中每个元素的值。一旦求解出矩阵U和矩阵V中每个元素的值,我们就可以轻易地求出每个需求单位对每种物资的关联度。换言之,只要求解出矩阵U和矩阵V,我们就可以填充需求单位-物资矩阵X中的缺失值。

基于矩阵分解的协同过滤模型的解法和奇异值分解算法(Singular Value Decomposition,SVD)非常相似。奇异值分解算法在信息检索领域被广泛应用于挖掘文章的潜在语义,虽然SVD算法可对矩阵进行分解,但直接将其用于分解需求单位-物资矩阵X并不是十分合适,因为矩阵中大量的缺失元素将无法被处理,还可能会导致模型过拟合(Over-fitting)。若为了解决稀疏问题,将矩阵中所有缺失的元素全部设置为0显然也是不合理的。

为了解决稀疏的问题,许多工作只对观测值进行拟合,并通过添加正则化项(Regularized Term)来避免模型的过拟合。协同过滤模型通过最小化向量ui和向量vj的内积与观测值xij的平方误差(Square Error)来学习向量ui和向量vj中元素的值。算法的目标函数(Objective Function)为:

式中,集合D为所有观测值的集合;λ为正则化项的权重,是一个超参数(Hyper-parameter),用于控制正则化的程度。通过最小化目标函数F(D),我们可以学习出矩阵U和矩阵V中每个元素的值。

3.2 基于矩阵分解的协同过滤模型的学习算法

基于矩阵分解的协同过滤模型的学习算法主要有两种:随机梯度下降法(Stochastic Gradient Descent,SGD)和交替最小二乘法(Alternative Least Square,ALS)。

随机梯度下降法每次随机地从观测值的集合D里挑出一个需求单位和物资的二元组(i,j)。给定第i个需求单位和第j种物资以及这个需求单位对该物资的需求量xij,对(xij-uivjT)2+λ(‖ui‖2+‖vj‖2)关于ui和vj中的每一个元素求偏导,并沿着求出来的偏导数的反方向优化ui和vj中的元素。随机梯度下降法收敛速度快,但最后收敛到的解并不能保证是全局最优解,有可能是局部最优解。

基于矩阵分解的协同过滤模型的另一种常用的学习算法是交替最小二乘法。在每次迭代中,交替最小二乘法固定矩阵U中的元素的值,求解矩阵V,再固定矩阵V中的元素的值,求解矩阵U。在进行多轮迭代以后,矩阵U和矩阵V会收敛到一个稳定的解。

4 结语

科学地做好物资计划是物资管理系统运行的关键,因此,我们利用了推荐系统中的协同过滤算法来优化物资管理系统。基于矩阵分解的协同过滤系统通过分析历史上需求单位对物资的需求量的模式,用矩阵分解的方法学习需求单位和物资的潜在特征。这种方法所需要的仅仅是部分需求单位对部分物资的需求量,就可以预测全部需求单位对全部物资的需求量。

参考文献

[1]夏培勇.个性化推荐技术中的协同过滤算法研究[D].青岛:中国海洋大学,2011.

[2]刘青文.基于协同过滤的推荐算法研究[D].合肥:中国科学技术大学,2013.

[3]荣辉桂,火生旭,胡春华,等.基于用户相似度的协同过滤推荐算法[J].通信学报,2014,35(2):16-24.

[4]冷亚军,陆青,梁昌勇.协同过滤推荐技术综述[J].模式识别与人工智能,2014,27(8):720-734.

矩阵分解的常用方法篇7

以求解线性方程组的矩阵三角分解 (LU分解) 为代表的数值计算在现代科学研究和工程技术中得到广泛应用。当今的计算机在高速的CPU和低速的存储器之间设置了二级或更多级的高速缓冲存储器 (简称高缓, cache) , 高缓的存取速度很快, 但容量较小, 将计算机当前正在运算的数据或经常需要访问的数据置于高缓中, 让数据存取的速度适应CPU的处理速度, 从而提高计算机整体的运算效率。对于大型方程组, 其系数矩阵很大, 标准LU分解直接对大矩阵进行计算, 一次参与运算的数据量很大, 不能很好地利用当今计算机的高缓, 运算效率不高。本文对LU分解运用分块算法, 将大矩阵的运算分解成分块子矩阵的运算, 使各部分小的分块矩阵的数据运算能够更多地在高缓中进行, 同时采取矩阵转置、循环展开等优化技术, 对分块算法中各个模块的计算设计了较高效的算法, 提高了算法的总体运算速度。

1 LU分解分块算法

LU分解用于求解线形方程组Ax=b。从高斯消去法出发, 其系数矩阵A可分解为单位下三角阵L和上三角阵U之积:

或:

A=LU (2)

利用分解后的L、U, 通过回代可方便地求出方程组的解。由矩阵乘法可推导出:

式 (3) 即为标准LU分解的计算公式, 整个数据计算在大矩阵A的范围内进行。

下面介绍LU分解分块算法。将矩阵A分解成2×2块, 方程A=LU可写成:

由矩阵乘法有:

A11=L11U11 (5)

A12=L11U12 (6)

A21=L21U11 (7)

A22=L21U12+L22U22 (8)

LU分解分块算法的步骤如下:

① 对式 (5) 中的分块矩阵A11进行LU分解, 可得L11、U11。

② 将①中得到的L11求逆, 并代入式 (6) , 有U12 = L11 -1A12 。

③ 将①中得到的U11求逆, 并代入式 (7) , 有L21 = A21 U11 -1。

④ 改写式 (8) , 设A′22=A22-L21U12, 则式 (8) 成为A′22=L22U22, 即为A′22的LU分解。

对A′22的LU分解可继续采用上述分块过程进行计算。这样, 只要每次选取适当的分块大小, 分块后的各子矩阵间的运算就有可能在高缓中进行, 从而提高运算速度。

2 LU分解分块算法的设计与实现

综上所述, LU分解分块算法的实现涉及六个模块的计算, 分别是分块矩阵的LU分解、单位下三角阵求逆、上三角阵求逆、单位下三角阵乘矩形阵、矩形阵乘上三角阵和矩形阵乘矩形阵。

2.1 分块矩阵LU分解

分块矩阵 (A11) 的LU分解计算采用式 (3) 所示的标准LU分解公式。与标准LU分解算法不同之处在于: (1) 形参的个数不同。对于分块分解, 需提供分块矩阵入口地址p, 阶数d以及在原矩阵内的偏移量displace和原矩阵的阶数n。 (2) 选主元的范围不同。标准LU分解是在大矩阵A的范围内选主元, 而块内LU分解是在分块矩阵的范围内选主元。

2.2 单位下三角阵求逆

可证明单位下三角阵的逆矩阵仍为单位下三角阵, 上三角阵的逆矩阵仍为上三角阵且对角线上的元素分别变为原来的倒数。

因为L11为单位下三角矩阵, 求逆后对角线上元素仍为1, 所以在计算时, 只需对对角线以下部分求逆。设单位下三角阵为L= (lij) d×d, 其逆矩阵为B= (bij) d×d。由矩阵乘法 $\sum_{k = j}^{i} (l_{i k} \times b_{i j}) = 0 (j < i)$ , 有 $b_{i j} = - \sum_{k = j}^{i - 1} (l_{i k} \times b_{k j}) (i = 2, 3, \dots, d$ ; j=2, 3, …, i-1) , 对角线上元素为bii=lii=1。在该算法中可直接将算出的bij保存在原lij的位置上。

2.3 上三角阵求逆

设上三角矩阵为U= (uij) d×d, 其逆矩阵为C= (cij) d×d。同样有 $c_{i i} = \frac{1}{u_{i i}} ‚ \sum_{k = i}^{j} (u_{i k} \times c_{k j}) = 0$ , 得到 $c_{i j} = - \frac{1}{u_{i i}} \sum_{k = i + 1}^{j} (u_{i k} \times c_{k j}) (i = d, d - 1, \dots, 1$ ; j=d, d-1, …, i+2) 。计算出的cij直接保存在原uij的位置上。

上面分别设计了针对单位下三角阵和上三角阵的求逆算法, 而没有采用常用的原地全选主元高斯约旦消元法, 这样做可以充分利用三角矩阵的特点, 避免无效的计算。且高斯约旦消元法对矩阵求逆时, 要多次全选主元, 进行很多次比较, 新设计的两个算法中不需要进行比较。

2.4 矩阵相乘

矩阵相乘是LU分解分块算法中的核心运算, 提高矩阵相乘的运算速度是提高整个分块算法效率的关键。本文采用了矩阵分块、矩阵转置和循环展开等方法来提高矩阵相乘的运算速度。

矩阵分块将大矩阵逐次一分为二, 直到划分到分块矩阵小于指定的大小。把这些分块矩阵当作元素来处理, 分块矩阵内部用普通三重循环来计算。设两相乘矩阵为A、B, 则过程如下所示:

$(\begin{matrix} A_{11} & A_{12} \\ A_{21} & A_{22} \end{matrix}) (\begin{matrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{matrix}) = (\begin{matrix} A_{11} B_{11} + A_{12} B_{21} & A_{11} B_{12} + A_{12} B_{22} \\ A_{21} B_{11} + A_{22} B_{21} & A_{21} B_{12} + A_{22} B_{22} \end{matrix})$

测试表明, 采用矩阵分块相乘技术可使运算速度得到较大提高。但对矩阵相乘仅仅采用分块还不够, 有必要引入矩阵转置和循环展开技术。

矩阵转置对于矩阵Am×pBp×n=Cm×n, 无论分不分块, 最终都要用三重循环实现。对C语言来说, 是采取按行存放元素, 最内层循环内A、C是顺序存取, 而B是不连续的。当矩阵的阶数较大时, 访问B中上下相邻两个元素跨度很大, cache 内可能放不下, 访问时不能命中, 仍要访问内存。这样的访问会进行m次, 当m很大时, 会降低效率。所以将B转置, 在最内层循环中用A的一行乘B的一行, 这样A、B、C都是顺序存取, 提高了数据访问速度。

循环展开循环展开可减少测试循环终止条件的次数, 它不仅可以减少循环的开销, 还可以在两个方面提高程序的性能。一是使程序利用多个执行部件并行计算, 充分发挥现代计算机多执行部件的特性。二是循环展开有助于重用装入寄存器中的数据, 降低了寄存器与cache以及cache 与内存之间对通信带宽的需求。展开多层循环的外层循环可以使一些变量得到重用。

下面介绍三处矩阵相乘的具体实现。

2.4.1 单位下三角阵乘矩形阵

计算U12 = L $_{11}^{- 1}$ A12 , 先将A12转置, 然后相乘, 结果保存于一临时矩阵T中。计算时需使用一大小为m的额外向量空间, 来保存A12转置后矩阵当前循环到的一行。乘运算从右边矩阵中的第一行开始, 逐行向下进行。采用循环展开技术, 第一层和第二层循环各展开1次。此时需要考虑两矩阵的行数为奇数的情况。当左边矩阵行数为奇数时, 计算结果存入右边矩阵的最后一列中;当右边矩阵行数为奇数时, 计算结果存入右边矩阵的最后一行中。

由于L $_{11}^{- 1}$ 的对角线上元素全为1, 所以它们不需和T中的对应元素相乘, 可直接将这些元素作为最内层循环乘积项累加的变量初始值, 这样可减少m×n次乘操作。又L $_{11}^{- 1}$ 为下三角阵, 只需将其对角线以下元素与T中对应元素相乘即可。

2.4.2 矩形阵乘上三角阵

在计算L21 = A21 U $_{11}^{- 1}$ 时, 将结果保存在左边矩阵 (A21) 中, 使用一大小为n的额外向量空间来保存左边矩阵当前循环到的一行, 从左边矩阵的第一行开始, 逐行向下运算。计算时, 第一层和第二层各循环展开1次。同样要考虑两矩阵的行数或列数为奇数的情况。当左边矩阵行数为奇数时, 计算结果存入左边矩阵最后一行;右边矩阵列数为奇数时, 结果存入左边矩阵最后一列。

在该模块中, 矩形阵A21的列数大于等于上三角阵U1l的列数, 且U11对角线以下元素全为 0, 所以没有对其转置, 只需将三角阵对角线以上元素与左边矩形阵中对应元素相乘即可。

2.4.3 矩形阵乘矩形阵

式A′22=A22-L21U12中, L21U12为两矩形阵相乘。计算时, 对L21进行分块, 将U12转置并分块。由2.4.1节中得到的结果T可直接作为L21所乘的对象, 并在做乘法时将所得的积取反, 相当于完成减运算, 这些操作有利于提高整个算法的效率。

对L21、T分块时, 取m=n=d, 即分块矩阵为d阶方阵。这样一直分块下去, 直至L21最后所剩一块矩阵的行数为N%d+d, 列数为d;T最后所剩一块矩阵的行数为d, 列数为N%d+d, 其中N为原矩阵A的阶数。

分块子矩阵相乘时采用循环展开技术, 第一层和第二层循环各展开1次。此时由于参与计算的分块矩阵均为稠密型矩阵, 可直接展开不必象有三角阵参与那样考虑边界处的处理, 计算结果最终保存在L21分块对应的位置上。

3 测试与分析

3.1 测试结果与比较

两种算法在PC机上测试, CPU为Pentium 4, 主频2.8GHz, 内存512MB, 编译系统采用VC++6.0。在验证算法计算正确的基础上, 依次取矩阵的阶数为300～6000, 对两种算法进行运行和比较。对于分块算法, 取每次分块大小d=150。测得的算法运行时间如表1所示, 表中第四行为分块算法比标准算法效率提高的百分比。

表1 (续)

由表1可见, LU分解分块算法明显优于标准算法, 在大矩阵情况下, 运算速度提高50%以上。

3.2 测试结果分析

标准LU分解算法和LU分解分块算法中的运算时间主要用在乘操作和除操作上, 可分别计算出两种算法中乘操作和除操作的运算次数, 限于篇幅在此只给出结果, 作简要分析。

分块算法 (取每次分块大小为d) 的乘运算量约为n3/3+ (n%d-0.5) n2, 除运算量约为 (d+1) n/2。标准算法的乘运算量为n (2n-1) (n-1) /6, 除运算量为n (n-1) /2。当矩阵阶数n很大时, 两算法的乘运算量相差不大, 而分块算法除运算量比标准算法降低了一阶。计算机做除法最耗时, 且标准算法没能充分利用高缓。分块算法在较好地利用高缓的同时, 通过矩阵分块、矩阵转置和循环展开等技术加速了矩阵乘法的运算, 所以运算效率显著提高。

LU分解分块算法比标准LU分解算法多用约d×n的存储空间。

4 结束语

本文将分块算法引入LU分解, 并加以有效实现, 使算法的运算速度得到较大提高。由算法分析和测试可知, 在分块算法中, 两矩形阵相乘的用时最多, 约占总时间的70%。因此, 要想进一步提高运算速度, 必须提高两矩形阵相乘的效率。可以对矩阵相乘作并行处理, 事实上在LU分解分块算法中, 有些模块是可以并行计算的。再者就是研究出达到最佳运算效率时总阶数n与分块大小d之间、cache大小与d之间的关系。以上这些问题将是以后进一步研究的目标。

摘要：对稠密型线性方程组的系数矩阵进行分块LU分解, 更充分地利用高速缓存, 提高运算效率。对LU分解分块算法进行了研究, 用VC++6.0对分块算法进行实现, 并与标准的LU分解算法进行比较。在大矩阵情况下, 分块算法比标准算法运算速度提高50%以上。

关键词：LU分解,矩阵分块,矩阵快速相乘,VC++6.0

参考文献

[1]陈建平.LU分解递归算法的研究[J].计算机科学, 2004, 31 (6) :141-142.

[2]李玉成, 朱鹏.BLAS的加速方法与实现技术[J].数值计算与计算机应用, 1998, 19 (3) :227-240.

[3]李玉成.LAPACK中的分块算法及其效果[J].数值计算与计算机应用, 2001, 22 (3) :172-180.

[4]李忠泽, 陈瑾, 龙翔, 等.基于PentiumPro的高性能BLAS的设计与实现[J].北京航空航天大学学报, 1998, 24 (4) :454-457.

[5]王小牛, 冯百明.基于存储的矩阵乘积优化算法[J].西北师范大学学报, 2005, 41 (1) :22-24.

[6]刘萍.数值计算方法[M].2版.北京:人民邮电出版社, 2007:85-88.

[7]Elmroth E, Gustavson F, Jonsson I, et al.Recursive Blocked Algorithms and Hybrid Data Structures.SIAMReview, 2004, 46 (1) :3-45.

[8]Toledo S.Locality of Reference in LUDecomposition with partial pivo-ting[J].SIAM Journal of Matrix Analysis&Application, 1997, 18 (4) :1065-1081.

【矩阵分解的常用方法】推荐阅读：