沈家銘

(四川大學(xué) 數(shù)學(xué)學(xué)院， 四川 成都 610065)

一種判斷LR分解是否存在的優(yōu)化算法

沈家銘

(四川大學(xué) 數(shù)學(xué)學(xué)院， 四川 成都 610065)

使用定理直接判斷方陣A是否存在LR分解的計(jì)算難度系數(shù)為O(n3)，文中在此基礎(chǔ)上提出了一個(gè)改進(jìn)算法,將計(jì)算難度降為O(n2)。給出了設(shè)計(jì)思路及推導(dǎo)方法，并用Matlab實(shí)現(xiàn)了兩種算法，通過(guò)例題驗(yàn)證了計(jì)算效率的提升。

LR分解; 方陣; 算法; Matlab

0 引 言

在有應(yīng)用背景的數(shù)學(xué)問(wèn)題中，很多求解問(wèn)題最終都可歸結(jié)為線性方程組的求解。因此，解線性方程組在科學(xué)與工程計(jì)算中有著十分重要的作用。對(duì)其系數(shù)矩陣進(jìn)行三角分解，是一種解線性方程組的有效方法，LR分解是求三角分解的一個(gè)強(qiáng)有力的算法,受到了人們的極大關(guān)注。在不進(jìn)行行列置換的前提下，為了判斷方陣A的LR分解是否存在，是進(jìn)行LR分解的前提[1-4]。文中給出了一種快速判斷方陣LR分解是否存在，同時(shí)計(jì)算LR分解的算法,大大提升了計(jì)算速度。

1 問(wèn)題描述

定義：LR分解：

即把矩陣A分解為一個(gè)下三角矩陣L和一個(gè)上三角矩陣R

而如何求證方陣A是否存在唯一的LR分解，有以下定理:

定理1 設(shè)A為n階方陣，A的k階主子式記為

則A的LR分解存在唯一的充分必要條件是

記

2 計(jì)算方法

為此，需要計(jì)算n-1個(gè)矩陣行列式的值。而容易證明，矩陣行列式的值就等于矩陣LR分解得到的R矩陣的對(duì)角元的乘積。

同時(shí)， d1≠0時(shí)可以得到A2的LR分解存在，由此算得d2，而若d2≠0，又已知d1≠0，所以A3的LR分解存在，以此類推，可以得到dk(k=1,2,…,n-1)的值。

由這個(gè)思路，可以得到算法1的計(jì)算步驟如下：

1)對(duì)Ak用Doolittle算法做LR分解；

2)由U計(jì)算dk；

3)若dk≠0，則k=k+1，否則返回錯(cuò)誤，不能進(jìn)行LR分解。按此思路即可得到算法1的Matlab程序代碼如下：

function[L,R,pp] =LR_P(A)

tic;

[n,m]=size(A);

R=zeros(n);

L=eye(n);

suM=0;

R(1,1)=A(1,1);

pan=R(1,1);

pp=0;

ifpan==0;

pp=-1;

return

end

fori=2:n;

forj=1:i-1;

fork=1:j-1

suM=suM+L(j,k)*R(k,i);

end

R(j,i)=A(j,i)-suM;

suM=0;

fork=1:j-1

suM=suM+R(k,j)*L(i,k);

end

L(i,j)=(A(i,j)-suM)/R(j,j);

end

R(i,i)=A(i,i)-L(i,1:i-1)*R(1:i-1,i);

pan=pan*R(i,i);

ifpan==0&&i～=n;

pp=-1;

return;

end

end

toc;

end

Doolittle分解算法代碼:

function[L,R]=LR_L(A)

[n,m]=size(A);

R=zeros(n);

L=eye(n);suM=0;

R(1,1)=A(1,1);

fori=2:n;

forj=1:i-1;

fork=1:j-1

suM=suM+L(j,k)*R(k,i);

end

R(j,i)=A(j,i)-suM;

suM=0;

fork=1:j-1

suM=suM+R(k,j)*L(i,k);

end

L(i,j)=(A(i,j)-suM)/R(j,j);

end

R(i,i)=A(i,i)-L(i,1:i-1)*R(1:i-1,i);

end

end

3 計(jì)算方法的改進(jìn)

n次LR分解的計(jì)算量非常大，是一個(gè)十分耗費(fèi)時(shí)間的過(guò)程，其時(shí)間難度為O(n3)，為此，尋找是否有優(yōu)化算法使得計(jì)算量減少[5-7]。

通過(guò)上述算法可以知道，每次LR分解的矩陣之間是有聯(lián)系的。它們之間的關(guān)系是：前一個(gè)矩陣增加一行一列就是后一個(gè)矩陣[8]，故進(jìn)行嘗試，已知矩陣A的LR分解式，是否能通過(guò)LR計(jì)算在末尾增加一行一列的矩陣A*的LR分解式L*，R*。

3.1改進(jìn)算法的推導(dǎo)

設(shè)矩陣A的LR分解式已知，不妨設(shè)增加一行一列以后的矩陣為：

其中

使用待定系數(shù)法，設(shè)

其中

則有

所以

即

可解得

然后可以得到

則通過(guò)A的LR分解式得到了在末尾添加一行一列A′的LR分解式A′=L′R′。

通過(guò)這種算法，使得在判斷A是否存在LR分解式，以及同時(shí)計(jì)算A分解式的LR的時(shí)間難度從O(n3)變?yōu)榱薕(n2)。

觀察解出來(lái)的β,γ的形式可以發(fā)現(xiàn)，每次對(duì)A加一行一列以后做LR分解，并不影響之前的分解結(jié)果，而且在算β,γ的時(shí)候推出來(lái)的遞推公式與Doolittle算法中對(duì)LR分解的算法遞推公式形式幾乎一致。通過(guò)比較得出了以下結(jié)論:改進(jìn)的算法(以上)與原算法區(qū)別只是計(jì)算的次序不同。

改進(jìn)的算法是從主對(duì)角元開(kāi)始，一圈一圈向外計(jì)算，也就是說(shuō)a11開(kāi)始，計(jì)算r11,然后計(jì)算r11周邊的元素r12,r22,l21,…。

每算完一圈，利用∏rii是否為0判斷能否進(jìn)行LR分解，若進(jìn)行到某一步時(shí)，該式為0，則該矩陣不能進(jìn)行LR分解，反之，則可以。

3.2改進(jìn)算法的實(shí)現(xiàn)

我們得到如下改進(jìn)的LR分解算法,稱之為算法2：

算法2的Matlab代碼如下：

function [L,R,pp] = LR_P2(A)

tic;

p=1;

[n,～]=size(A);

for i=1:n

if p==0

pp=-1;

break;

end

p=1;

A_=A(1:i,1:i);

[L,R]=LR_L(A_);//Dolittle算法，for j=1:i

p=p*R(i,i);

end

pp=0;

toc;

end

end

4 計(jì)算效率驗(yàn)證

使用算法1和算法2的代碼做了一個(gè)簡(jiǎn)單的例題來(lái)驗(yàn)證計(jì)算，檢驗(yàn)計(jì)算速度是否有提高。

對(duì)一個(gè)1 024×1 024的矩陣進(jìn)行LR分解，矩陣來(lái)源是一張1 024×768的點(diǎn)陣圖片。

在平臺(tái)為intel 2640qm,內(nèi)存8 GB，Matlab環(huán)境下，統(tǒng)一計(jì)算平臺(tái)下對(duì)1 024×1 024的矩陣用兩個(gè)算法進(jìn)行對(duì)比。

算法1：Elapsed time is 2 331 seconds。

算法2：Elapsed time is 12 seconds。

可見(jiàn)效率提升了194.25倍。

5 結(jié) 語(yǔ)

在定理的基礎(chǔ)上設(shè)計(jì)實(shí)現(xiàn)了一種高效的判斷LR分解是否存在以及計(jì)算LR分解式的算法。算例測(cè)試證明了算法可行性和效率的提高。

[1] 高惠璇.統(tǒng)計(jì)計(jì)算[M].北京：北京大學(xué)出版社,1997.

[2] 徐曉飛,曹祥玉,姚旭,等.一種基于Doolittle LU分解的線性方程組并行求解方法[J].電子與信息學(xué)報(bào),2010(8):2019-2022.

[3] G H Golub. Matrix computation 4th edition[M].北京：人民郵電出版社,2014.

[4] 王群英.矩陣分解方法的探究[J].長(zhǎng)春工業(yè)大學(xué)學(xué)報(bào):自然科學(xué)版，2011，32(1):95-101.

[5] 戴華.矩陣論[M].北京：科學(xué)出版社，2001：110-145.

[6] 徐泰燕，郝玉龍.非負(fù)矩陣分解及其應(yīng)用現(xiàn)狀分析[J].武漢工業(yè)學(xué)院學(xué)報(bào)，2010(1):110-115.

[7] 黃麗嫦.矩陣的LU并行遞歸分解算法的設(shè)計(jì)研究[J].科學(xué)技術(shù)與工程，2012(5)：3626-3629.

[8] 周康，陳金.基于部分基變量的LP問(wèn)題矩陣算法[J].運(yùn)籌學(xué)學(xué)報(bào)，2012(6)：121-126.

AnoptimizationalgorithmtodeterminetheexistenceofLRdecomposition

SHENJia-ming

(SchoolofMathematics,SichuanUniversity,Chengdu610065,China)

SincethecalculatingdifficultycoefficientofLRdecompositionisO(n3),byusingthetheoremtojudgewhetherphalanxAexists,hereweputforwardanimprovedalgorithmtoreducethecalculatingdifficultycoefficienttoO(n2).Designandderivationmethodaregiven,andtwoalgorithmsareimplementedwithmatlab.Thecalculationefficiencyimprovementisverifiedwithexamples.

LRdecomposition;phalanx;algorithm;Matlab.

2014-07-18

沈家銘(1993-)，男，漢族，云南昆明人，主要從事統(tǒng)計(jì)計(jì)算方向研究,E-mail:493879714@qq.com.

O 151.21

A

1674-1374(2014)05-0593-04