LU分解MatLab算法分析

最近矩阵分析老师出了一道题目作为作业，是一道程序题，题目是对A矩阵做LU分解，要求能对A实现PA=LU的分解，并最终输出L，U，P矩阵。先来解读下题目，寥寥几句话，里面囊括的信息量却不少，然后这些都得自己去琢磨。首先对A矩阵能做LU分解，即能把A分解成这种形式A=L*U（U是上三角矩阵，是由A矩阵经过高斯消元后得到的，L是下三角矩阵，其对角线全为1，其他非零元素为在消去（i,j）位置元素过程中主元所乘的系数），条件有3，一是矩阵A必须为方阵，A如果不是方阵，就不要想着对它做LU分解啦，这是基本条件，牢记啊！二是矩阵A必须可逆，换种说法就是A必须为非奇异矩阵，这两种说法是等价的，而这又等价于A是满秩的，A是满秩又等价于A的行列式值非0，好绕，矩阵就是这样，很多定理其实是等价的，但是你得记住，不然在推导一些定理或公式的时候会犯一些基本的常识性错误。三是矩阵A在高斯消元过程中必须没有出现0主元，也就是说只有在对A进行高斯消元过程中没有出现进行交换这种情况下，A才能分解成L*U这种形式，如果对A进行高斯消元，中间某一步出现0主元，需要进行行交换了，这种情况下就不要想对A进行LU分解啦，因为不满足条件3啊！那么问题来了，假如出现了有0主元这种情况，我又想对A进行LU分解，那应该怎么办？这就引出了带行变换的LU分解，也就是本文的主题。根据书上的定理，对任意一个非奇异矩阵（等价于可逆矩阵）都存在一个置换矩阵P使得P*A可以分解成L*U这种形式，即PA=LU。想想其实这定理也是不言自明的，刚才A不是要进行行变换才能继续高斯消元吗？而LU分解前提又是高斯消元过程中不能出现行交换，那好，我事先对A矩阵在高斯消元过程中需要交换的行给交换掉，形成一个新的矩阵B，那我对B高斯消元那肯定就不会出现需要行交换的情况，这就满足了LU分解第三个条件了，这样B不就可以进行LU分解了吗？是的，PA=LU这种形式的LU分解采用的就是这种思想。那么现在的问题是，怎么在代码中实现对A矩阵的LU分解，并输出P，L，U矩阵呢？我在网上搜了一下，发现结果大都不尽如人意，大多数程序吧只能说做A=LU这种形式的分解，一旦说A不满足条件3，那就死翘翘了，这种程序先不论其能否运行成功，结果是否正确，其鲁棒性也太差了！用个时髦点的词来说就是太low了！通用性太差了！不光如此，代码也没什么注释，可读性很差，让人怀疑是不是写给别人看的，尤其像我这种编程渣渣，看个代码费老半天劲都看不懂说什么。于是，我决定按自己的想法来走。首先从最简单的情况考虑，这也是我们做研究、做学术、做工程必须要时刻牢记心中的一点，很多人喜欢一上来就把所有问题、把最复杂的情况、把方方面面都给考虑到，然后再开始实现他的想法，我自己也有这个习惯，但是，这并不是一个好习惯，一上来就好高骛远、就想着高大上这本质上是一种急功近利的表现，那样的话你会陷入到各种各样的技术细节当中，你会想半天却仍然写不出半点实质性的东西出来，所以最好的办法是，先考虑最简单、最核心的情况，这样不仅大大降低问题的复杂度，同时也为将来进一步扩展程序、解决更复杂的情况打下了一个坚实的基础。在这个例子中，最简单的情况就是矩阵A在高斯消元过程中不需要进行行交换，也就是说A可以分解成A=L*U这种形式。这种情况下，代码如下。functionLUDecomposition(A,n)%Aisthesquarematrix,nistheorderofAL=eye(n);%LettheLmatrixbeanidentitymatrixatfirstfori=1:n-1forj=i+1:nL(j,i)=A(j,i)/A(i,i);A(j,:)=A(j,:)-(A(j,i)/A(i,i))*A(i,:);endendU=A%AbecomesUmatrixafterGausseliminationL可以试着令A=[222;477;61822]，调用函数获得L矩阵为[100;210;341]，U矩阵为[222;033;044]，用笔验算下，这个结果是正确的。代码运行结果如图所示：这部分代码的主要思想是这样的，矩阵A的阶次为n的话，A在高斯消元后有n个非零主元。在消元过程中，A共需要消掉n-1个主元下面所有的元素，注意，第n个主元已经是矩阵的最后一个元素了，它的下面和右边都没有其他元素了，所以不存在说对第n个主元下面所有元素消去的情况。这就获得了我们代码的第一个for循环，从第1行主元开始消元，一直到第n-1行主元。而在获得每一行主元过程中，需要对该行主元下面所有元素都消去，假如现在要获得第i行主元的话，就是说要对该主元所在列的第i+1行到第n行元素都消掉，那么这就获得了我们代码的第二个for循环，从消去第i+1个元素开始一直到第n个元素。前文说过，消掉第（j，i）个位置元素过程中，主元所乘系数就是L矩阵第（j，i）位置的元素，所以有L(j,i)=A(j,i)/A(i,i);。然后的话，就是把A矩阵第j行减去第j行乘以L（j，i），这样就可以消掉第（j，i）个元素了，就是这行代码A(j,:)=A(j,:)-(A(j,i)/A(i,i))*A(i,:);。最后，执行完两层for循环后，A矩阵就成为了U矩阵，L矩阵也从最初的单位阵成了L矩阵。好了，我们已经实现了最简单的情况了，下面考虑复杂点的情况，就是说对A进行PA=LU这种形式的分解。假如A在消元过程中出现0主元了，那么怎么办？很简单，只需要从该0主元下面所有元素中找到一个非0元素，然后将其所在的行与该0主元所在的行进行交换就行了，不要忘了，对A矩阵两行进行了交换，对应到P矩阵中的操作是相应的两行也要进行交换，因为我们是通过P矩阵两行交换后然后左乘A矩阵使得A矩阵两行进行交换的。A矩阵交换第i行和第k行元素对应到L矩阵中相应两行的消元系数也应该交换位置，就是说L矩阵的第i行和第k行元素也要交换位置，当然，主对角线上的1是不需要交换的，因为他们并不是消元系数。交换完成后，继续执行消元操作，其步骤和上面考虑的最简单的情况就是一样的了。就这样，我们就实现了PA=LU这种形式的分解。令A=[12-34;4812-8;2321;-3-11-4]，代入函数运算得L矩阵为[1000;-3100;2-0.210;2-0.210;403.751],U矩阵为[12-34;05-88;006.4-5.4;000-3.75]，P矩阵为[1000;0001;0010;0100]，用笔验算下，结果与函数运行结果是一致的，当然了，这个函数我只是代了3，4个不同的A矩阵进去而已，可能样本数量不够多，但目前来说我觉得应该没什么问题了，如果有问题欢迎反馈给我。这部分代码如下：functionAdvanceLUDecomposition(A,n)%Aisthesquarematrix,nistheorderofA,AmustbeinvertibleD=A;%StorematrixAinD,forlateruseL=zeros(n);%LettheLmatrixbeanzeromatrixatfirstP=eye(n);%Letthepermutationmatrixbeaidentitymatrixatfirstfori=1:n-1forj=i+1:nifA(i,i)==0%Azeropivotappearson(i,i)position,weneedtofindanonzeroentrybelowittobethenewpivot,withrowexchangefork=n:-1:i+1%findanonzeroentrybelowthe(i,i)entryintheicolumn,startfromthelastrowifA(k,i)~=0%Wehavefoundanonzeroentry,tochooseitasthenewpivot,weneedrowexchangek--iL([ik],:)=L([ki],:);%PermuteiandkrowinLmatrixA([ik],:)=A([ki],:);%PermuteiandkrowinAmatrixP([ik],:)=P([ki],:);%PermuteiandkrowinPmatrixbreak;endendendL(j,i)=A(j,i)/A(i,i);A(j,:)=A(j,:)-(A(j,i)/A(i,i))*A(i,:);endendU=A%AbecomesUmatrixafterGausseliminationL=L+eye(n)%AllentriesonthediagonalofLmatrixmustbe1P%outputthepermutationmatrixB=L*U%verifyiftheproductofLandUequalstoP*AC=P*D%DistheoriginalAmatrix,checkitoutinrow2%IfBequalsC,thenitmeansthealgorithmworkscorrectly%somekeypointsandtheromsaboutLUfactorization%Theorem1AnonsigularmatrixAnxnpossessesanLUfactorizationifand%onlyifanonzeropivotdoesnotemergeduringrowreductiontoupper%triangularformwithtypeIIIoperations.%Theorem2ForeachnonsigularmatrixA,thereexistapermutationmatrixP%suchthatPApossessesanLUfactorizationPA=LU%Remember,theconceptofnonsingularmatrixisforsquarematrix,itmeans%thatthedeterminantisnonzero,andthisisequivalentthatthematrixhas%full-rank%Basedontheseconditions,thefirstthingaboutthematrixAonwhichwe%conductLUfactorizationisthatAmustbeasquarematrix.Thesecond%thingisAmustbeinvertible,whichisequaltothestatementthatAis%non-singular代码运行结果如图所示：最后补充一点，为什么要进行LU分解呢？这个问题很关键，很多人也许并不关注这个问题，我们学习很多时候都是只关注实现方法，却并不关心它存在的意义，这种学习是永远无法深入的，只能是停留在表面上，学习就应该多问为什么，多质疑这个东西存在的价值，存在的意义有多大，这样才能促使你去深入了解这个方法的优点和缺点，从而改进、完善它。简单点来说就是LU分解大大降低了算法复杂度，我们求解一个方程组Ax=b的时候，一般来说无非就两种方法，要么是高斯消元法，要么是先求A的逆矩阵，然后再乘以b获得x，而第二种方法比第一种方法要复杂并且限制更多，所以一般是用高斯消元法。高斯消元法解一个方程组算法复杂度是(n^3)/3，并且每获得一个新的b，要接x，都得执行复杂度为(n^3)/3的操作。而LU分解有什么好处呢？在第一次LU分解的时候，也就是说获得L和U的时候，其算法复杂度其实也是n^3,但是，一旦我们获得了L和U矩阵后，每次我们获得一个新的b要求对应的解x，算法复杂度就会大大降低，粗略来说就是n^2，把复杂度降低了一个级别，对于大型系统来说，这是非常了不起的一个改进，运算性能会大大提升。而实际应用中，这样的方法也是非常有意义的，实际系统中，A矩阵相当于系统里的各种滤波和变换操作，x相当于系统的输入，b相当与系统的输出，我们一般是获得了输出b，然后想求得输入x，只要系统不变，那么知道b，又知道了L和U矩阵，我们只需要对每