王 芳，馬艷麗

(1.安徽外國(guó)語(yǔ)學(xué)院 信息與數(shù)學(xué)學(xué)院，安徽 合肥 231200；2.安徽新華學(xué)院 通識(shí)教育部，安徽合肥 230088)

1 相關(guān)定義及引理

定義1[1]若矩陣A∈Rm×n的秩為r，則A必有一個(gè)非奇異的子矩陣A11∈Rr×r，稱A11為矩陣A的最大線性無(wú)關(guān)子塊。

引理2[3]r維向量組的每一個(gè)向量添加n-r個(gè)分量成為n維向量。如果r維向量組線性無(wú)關(guān)，那么，n維向量組也線性無(wú)關(guān)。反言之，如果n維向量組線性相關(guān)，那么，r維向量組也線性相關(guān)。

2 尋找矩陣中最大線性無(wú)關(guān)子塊的兩種算法

由引理2可知對(duì)于n維向量組，判斷線性相關(guān)性，只需判斷n維向量組中的每一個(gè)向量去掉n-r個(gè)分量后所形成的r維向量組的線性相關(guān)性即可。因?yàn)閞+1個(gè)r維向量組一定線性相關(guān)，故在實(shí)際應(yīng)用中需考慮向量組中向量的個(gè)數(shù)。

(1)

其中

(2)

求出齊次線性方程組Λsx=0的一個(gè)基礎(chǔ)解系。其中

Λs=(h1,h2,h3,…h(huán)s)T

(3)

③矩陣B是一個(gè)行滿秩矩陣，行向量的先后順序與矩陣A中的有區(qū)別。而若已知矩陣A∈Rm×n，PA=B，其中矩陣B∈Rr×n的秩為r，P∈Rr×m是某置換矩陣的前r行，在MATLAB中輸入以下命令

B=P′*B;B(all(B==0,2),∶)=[]

可使B的行向量的先后順序與在矩陣A中的先后順序相等。

具體算法如下：

算法1

1)令size(A)=[m,n]，置初始量：m階單位矩陣Im×m，u=2，v=1，k=2，r=2,l=1,P、B為空矩陣，分別將矩陣Im×m和A的第一行移動(dòng)到矩陣P、B的第一行。

2)當(dāng)u小于等于n，令l1=1，取u階單位矩陣Iu×u，p(∶,u-1)=-a(u)×Iu×u(∶,1)+a(1)×Iu×u(∶,u)。若l1小于等于m-1，c=max(|A(l1,1∶u)×p|),如果c大于給定的ε，v增加1，并分別將矩陣Im×m和A的第l1行移動(dòng)到矩陣P、B的第v行。停止,l1增加1，如果l1大于m-1，u增加1，p(u,u-1)=0，c大于給定的ε，停止。

3)對(duì)t=u+1,……,n，取t階單位矩陣It×t，令p(t,t-1)=0，p(∶,t-1)=-a(t)×It×t(∶,1)+a(1)×It×t(∶,t),b=p,r=2;當(dāng)r小于等于v，取空矩陣Q，計(jì)算B(r,1:t)×b，并將其值賦予矩陣Q，令[c,i]=max{|Q|}，用g表示矩陣b的第i列。

對(duì)于j=1,2,……，t-r+1，計(jì)算b(∶,j)=b(∶,j)-Q(j)×g/Q(i)，去掉矩陣b的第i列，r增加1。

令size(A)=[m1,n1]，l=1。

若k小于等于t，l大于m1，取空矩陣Q，計(jì)算A(l,1∶t)×b，并將其值賦予矩陣Q，令[c,i]=max{|Q|}，用g表示矩陣b的第i列。

如果c大于給定的ε，v增加1，并分別將矩陣Im×m和A的第l行移動(dòng)到矩陣P、B的第v行，l減小1，m1減小1，令size(b)=[m2,n2]，n2-1小于等于0，停止運(yùn)行。

對(duì)j=1,……n2，計(jì)算b(∶,j)=b(∶,j)-Q(j)×g/Q(i)，去掉矩陣b的第i列，k增加1，終止，l增加1。

如果v-m大于等于零(B的行數(shù)大于等于A的行數(shù))，程序終止。

如果v-t大于等于零(B的行數(shù)大于等于A的列數(shù))，程序終止。

t增加1，令B=PTB，去掉B中的全零行；輸出B。

算法1所得到的矩陣B為行滿秩矩陣，它的行向量組是矩陣A的行向量組的一個(gè)極大線性無(wú)關(guān)組，假設(shè)size(B)=(r,n)，r即為矩陣A的秩，然后對(duì)BT采用算法1，再取轉(zhuǎn)置，即可得到最大線性無(wú)關(guān)子塊A11。顯然A11所對(duì)應(yīng)的行列式即為矩陣的最高階非零子式，若已知道rank(A)，則找到rank(A)維線性無(wú)關(guān)子塊即可中斷。

以上算法是先提取矩陣A的行向量組的一個(gè)極大線性無(wú)關(guān)組，然后形成最大線性無(wú)關(guān)子塊，接下來(lái)將考慮直接尋找最大線性無(wú)關(guān)子塊的算法，其主要思想是從二維最大線性無(wú)關(guān)子塊開始，已知r維線性無(wú)關(guān)塊(不一定是矩陣A的子式)，然后在其基礎(chǔ)上尋找r+1階線性無(wú)關(guān)塊。主要方法是在r的每一個(gè)行向量添加一個(gè)分量，形成含有r個(gè)r+1維向量的線性無(wú)關(guān)組，然后利用引理1尋找第r+1個(gè)r+1維向量，使其與已尋找到的r個(gè)r+1維向量線性無(wú)關(guān)，如果不存在，則互換矩陣A的兩列，并重新添加一個(gè)分量，重新尋找，直到找遍所有行與列即可停止，若已知rank(A)，則找到rank(A)維線性無(wú)關(guān)子塊即可中斷。

若已知矩陣A∈Rm×n且MAN=B，M、N是置換矩陣，B的秩為r,B的前r行、前r列構(gòu)成B的最大線性無(wú)關(guān)子塊，則A的最大線性無(wú)關(guān)子塊可通過(guò)在MATLAB中輸入以下命令得到：M(r+1∶m,∶)=0,M=MT;N(∶,r+1∶n)=0;N=NT;B=M*A*N;B(all(B==0,2),∶)=[];B(∶,all(B==0,1))=[]。

直接提取最大線性無(wú)關(guān)子塊的具體算法如下：

算法2

1)取初始量：m階單位矩陣Im×n,n階單位矩陣In×n，u=2，k1=0，hh=0，令v=min(m,n)。

2)若u小于等于v，取矩陣A第一行的前u列，u階單位矩陣In×n，令

p(∶,u-1)=-a(u)×Iu×u(∶,1)+a(1)×Iu×u(∶,u)，b=p,r=2。

若r小于等于u-1，取空矩陣G，計(jì)算A(r,1∶u)×b，并將其值賦予矩陣G，令[c,i]=max{|G|}，用g表示矩陣b的第i列。

對(duì)于j=1,2,……,u-r+1，計(jì)算b(∶,j)=b(∶,j)-G(j)×g/G(i)，去掉矩陣b的第i列，r增加,程序終止。

l=u，當(dāng)l小于等于m，令c=max(|A(l,1∶u)×b|)。

算法2輸出的矩陣B就是矩陣A的最大線性無(wú)關(guān)子塊A11。

3 數(shù)值實(shí)驗(yàn)

以下數(shù)值實(shí)驗(yàn)均在Intel(R)Core(TM)i5-8265U CPU @ 1.60 GHz 1.80 GHz內(nèi)存為8.00 GB的個(gè)人計(jì)算機(jī)上完成，所用軟件為MATLAB R2018a[4-6]。

方法1：

在MATLAB命令窗口輸入矩陣A=[1-1 2 1 0;2-2 4 2 0;3 0 6-1 1;0 3 0 0 1]，使用算法1可得行滿秩矩陣B：

令B=BT,對(duì)B使用算法1并取轉(zhuǎn)置可得：

方法2：

在MATLAB命令窗口輸入矩陣A=[1-1 2 1 0;2-2 4 2 0;3 0 6-1 1;0 3 0 0 1]，直接使用算法2即可得最大線性無(wú)關(guān)子塊：

數(shù)值實(shí)驗(yàn)2使用算法1、算法2提取矩陣(表1)，矩陣來(lái)自“Matrix Market”，皆為實(shí)數(shù)矩陣的最大線性無(wú)關(guān)子塊，其中ε取1.1×10-10，ti(i=1,2)依次表示兩種算法的CPU運(yùn)行時(shí)間(單位為s)。

表1 兩種算法的運(yùn)行時(shí)間 s

4 結(jié)論與討論

兩種算法雖可用來(lái)提取大型稀疏矩陣的最大線性無(wú)關(guān)子塊，但時(shí)間上并不占優(yōu)勢(shì)。所需CPU運(yùn)行時(shí)間第一種算法遠(yuǎn)遠(yuǎn)比第二種算法短，兩種算法CPU運(yùn)行時(shí)間與矩陣內(nèi)部的線性關(guān)系有很大聯(lián)系。兩種算法可清楚地看到矩陣內(nèi)部的線性無(wú)關(guān)性，即矩陣的行與列的局部線性無(wú)關(guān)性。另外如果需要求取矩陣的秩及提取最大線性無(wú)關(guān)子塊時(shí)的置換矩陣，只需輸出時(shí)稍作改動(dòng)即可。

我們?cè)谑褂没谇‘?dāng)分裂的預(yù)條件QMR算法和預(yù)條件GMRES算法[7-10]求解奇異線性方程組時(shí)，可使用這兩種算法提取最大線性無(wú)關(guān)子塊和置換矩陣來(lái)構(gòu)造預(yù)條件子。