河南師范大學數(shù)學與信息科學學院(453007)

龐善起 鹿姍姍

正交表的構(gòu)造方法及Matlab實現(xiàn)*

河南師范大學數(shù)學與信息科學學院(453007)

龐善起 鹿姍姍

正交試驗設計是研究多因素水平的試驗設計方法,它是一種高效率、快速、經(jīng)濟的試驗設計方法。正交表在衛(wèi)生、統(tǒng)計、醫(yī)藥、農(nóng)業(yè)、制造業(yè)等方面都有非常重要的應用,尤其是高強度正交表可以研究較多因素之間的相互作用,但是正交表的強度越高越難構(gòu)造;另外,混合正交表的魅力在于它有較大的靈活性,允許試驗因素具有不同水平數(shù),得到了越來越多的關注。因此如何快速、高效地得到需要的正交表是有必要和有意義的。

至今為止,已經(jīng)有許多組合數(shù)學家和統(tǒng)計學家致力于正交表的研究,例如楊子胥、龐善起、張應山、杜蛟等分別在文獻[1-8]中給出了不同類型的正交表的構(gòu)造方法,這些正交表不僅在試驗設計中有廣泛的應用,隨著研究的不斷深入,它們還可以被應用于編碼學、密碼學、計算機科學中。

正交表構(gòu)造方法與原理

1.Ltu(tm)型正交表[1]

(1,0,0,…,0,0,0)

1個

(bu-1,1,1,0,…,0,0,0)

t個

(bu-2,1,bu-2,2,1,…,0,0,0)

t2個

?

?

(b2,1,b2,2,b2,3,…,b2,u-2,1,0)

tu-2個

(b1,1,b1,2,b1,3,…,b1,u-2,b1,u-1,1)

tu-1個

因為每個bij為域GF(t)的所有元素,即有t種取法,所以以上的向量,就是Vu的全部標準向量,共有

構(gòu)造Ltu(tm)型正交表的方法和步驟如下:

第1步,先給出t階有限域GF(t);

第3步,以全體向量為行號,以全部的標準向量為列號,做出所有內(nèi)積即可得Ltu(tm)型正交表。

2.混合水平正交表

混合水平正交表[1]可由Ltu(tm)型正交表通過并列而產(chǎn)生。構(gòu)造混合水平正交表的方法和步驟如下:

第1步,可借助上文Ltu(tm)型正交表的方法得到Ltu(tm)型正交表;

第2步,從Ltu(tm)中任意取定k個獨立的列(k

第4步,再將這k列同行元素所構(gòu)成的tk個互異的有序組,按字典排列法與tk個自然數(shù)1,2,…,tk建立一一對應,而把這k列所構(gòu)成的每個有序組均換成它對應的自然數(shù),便得到一個新的tk水平列;

3.高強度正交表

一個第j 列的元素是0,1,…,qj-1的ω×n矩陣A,稱為一個強度為2的正交表,如果滿足以下條件[4]:

(1)每一列中每個元素出現(xiàn)的次數(shù)相同;

(2)在任兩列ai,aj(1≤i,j≤n)中的每一個數(shù)對(0,0),…,(0,qj-1),(1,0),…,(1,qj-1),…,(qj-1,0),…,(qj-1,qj-1)出現(xiàn)的次數(shù)相同。

只滿足上述條件(1)的矩陣叫做強度1的正交表。如果A的任何ω×d的子矩陣包含所有可能的1×d行向量,且這些行向量出現(xiàn)的次數(shù)相同,稱為強度為d的正交表。強度d≥3的正交表稱為高強度正交表[4-6]。

引理1 設M是一個選自Qn的任意t(1≤t≤qn-1)個行向量構(gòu)成的矩陣,則

是一個強度為1的正交表。

引理2 對于任意的整數(shù)1≤k≤q總存在Qn上的1-正交分劃A1,A2,…,Ak。

證明:由于Qn本身就是一個正交表,所以當k=1時定理顯然成立； 當k

由于Qn中有qn個q元n維向量,考慮將Qn分為行數(shù)相等的分劃容量為q的A1,A2,…,Ak,其中的每個Ai(1≤i≤k)中都有qn-1個向量,為了得到Ai:

?

第2q-1步,將剩下的qn-1個向量,組成矩陣Aq,它當然是一個強度1的正交表。

這樣就構(gòu)造出了Qn的一個1-正交分劃A1,A2,…,Aq。

引理3 若有Qn上的正交分劃A1,A2,…,Aq,如果Ai(i=1,2,…,q)都是qn-1行n列的強度為t的正交表,那么下列矩陣都是強度為t+1的正交表,并且是Qn+1的(t+1)-正交分劃。

其中的0,1,2,…,q-1均指與其對應的Ai有相同的行的元素全為0,1,2,…,q-1的向量,下文同。

高強度正交表的構(gòu)造可由引理1,2得到強度為1的正交表,和Qn的1-正交分劃,在這個前提下,再由引理3得到強度為2的正交表,和Qn+1的2-正交分劃,接著由引理3得到強度為3的正交表,和Qn+2的3-正交分劃,依次下去,逐步構(gòu)造出高強度的正交表。

正交表構(gòu)造的Matlab實現(xiàn)

為了更加直觀詳細地說明Matlab程序的可行性,下面我們將選取幾個較為簡單的正交表作為例子說明利用Matlab程序?qū)崿F(xiàn)正交表的構(gòu)造方法。

1.Ltu(tm)型正交表

clc;clear

t=2;u=3;%初始條件

p=t;V=[ ];K=[ ];L=[ ];D=[ ];L1=[ ];L2=[ ];

G=transpose(0:(t-1));%t階有限域

i=1;

while(i<=u)

R=kron(kron(ones(t^(i-1),1),G),ones(t^(u-i),1));

V=[V,R];

i=i+1;

end

V%所有的u維向量

U=rref(V′);%標準化

[m,n]=find(U==1);

[s,t]=find(U((m+1):u,n)<=1);

B=U(:,unique(n(t)))′;%標準向量

L=V*B′;

L=mod(L,p)%得到正交表L8(27)

在Matlab程序中只需修改初始條件t,u的值,我們就可以生成形如Ltu(tm)的正交表。

2.混合水平正交表

以L8(4×23)為例,則23=8,此時t=2,u=3,k=2,由上文的程序得到L8(27),取定2個獨立的列(2<3),去掉能用這2列線性表示的所有列,再將這2列所構(gòu)成的22=4個互異的有序組,把這2 列所構(gòu)成的每個有序組均換成它對應的自然數(shù),從而生成正交表L8(4×23),結(jié)果如表2,Matlab程序與解釋如下:

表1 正交設計表L8(27)

clc;clear

t=2;u=3;k=3;%初始條件

p=t;[z,z1]=size(L);%接上節(jié)程序

K=L(:,1:k);%默認1:k列,也可以換成其他列

i=1;

while(i<=z1)

d=rank([K,L(:,i)]′)%m個n維向量線性相關?秩

D=[D,d];

i=i+1;

end

D

a=find(D(1,:)>=(k+1))

L1=L(:,a)%得到正交表Ltu(tr)

%以下生成混合水平列

Y=unique(K,′rows′)%刪除K中重復的行,得到有序組

[y,y1]=size(Y);

j=1;

while(j<=y)

b=find(K(:,1)==Y(j,1)&K(:,2)==Y(j,2));

[z,z1]=size(b);

K(b,:)=(y-(j-1))*ones(z,k);%有序組與自然數(shù)建立一一對應

j=j+1;

end

L2=[flipud(K(:,1)),L1]%得到混合水平正交設計表

3.高強度正交表

以L8(24)為例,則Q={0,1},q=2,N=4,t=3,首先生成強度為1的正交表,和Q2的1-正交分劃,再生成強度為2的正交表,和Q3的2-正交分劃,接著生成強度為3的正交表,和Q4的3-正交分劃,這樣就可以生成2水平強度為3的正交表。結(jié)果如表3,和表4,Matlab程序與解釋如下:

表2 正交設計表L8(4×23)

clc;clear;

q=2;%水平(初始條件)

t=3;%強度(初始條件)

N=4;%列數(shù)(初始條件)

n=N-(t-1);%q水平強度為1的正交表的列數(shù)

Q=[];a0=[];a1=[];a=[];A1=[];A=[];L=[];L0=[];L1=[];

%第一部分：生成Q^n即所有的q元n維向量

i=1;

while(i<=n)

R=kron(kron(ones(q^(i-1),1),[0:(q-1)]′),ones(q^(n-i),1));

Q=[Q,R];%Q^n即所有的q元n維向量

i=i+1;

end

%第二部分：生成1-正交分劃

j=1;

while(j<=(q-1))

for i=0:(q-1);

a0=Q((1+q^(n-2)*(j-1)):q^(n-2)*j,:);

a1=mod(i+a0,q);%正交分劃

a=[a;a1];

end

Q=setdiff(Q,repmat(a,2,1),′rows′);%Q^n-a

j=j+1;

end

A=[a;Q];%生成強度為1的正交表

%第三部分：正交分劃的遞歸構(gòu)造

while(n<=N-1)%Q正交分劃q^(n-2)行向量的個數(shù)

L0=kron(ones(q,1),kron([0:q-1]′,ones(q^(n-1),1)));

j=0;

while(j<=q-1)

A1=circshift(A,[q^(n-1)*j,0]);%每次向下移q^(n-1)

L1=[L1;A1];

j=j+1;

end

A=[L0,L1];%所有正交分劃的并列

L1=[];

n=n+1;

end

%第四部分:生成強度為t的正交表

for i=1:q;

L=A((1+q^(n-1)*(i-1)):q^(n-1)*i,:)%強度為t的正交表

end

在Matlab程序中,保證n=N-(t-1)≥2的前提下,只需修改初始條件q,N,t的值,我們就生成q水平,強度為t,列數(shù)為N,行數(shù)為q^(N-1)的正交表。比如修改初始條件,使得q=4,t=3,N=4便可得到水平為4,強度為3,列數(shù)為4,行數(shù)為64的正交表;修改初始條件,使得q=5,t=3,N=4，便可得到水平為5,強度為3,列數(shù)為6,行數(shù)為125的正交表。

表3 正交設計表L8(24)

表4 正交設計表L8(24)

結(jié) 論

正交表研究正在引起眾多領域研究者的廣泛關注。本文在原有的理論基礎上使用Matlab 生成三種類型的正交表,但是仍有一些問題值得改進和研究，本文只給出了生成單一水平的高水平高強度正交表構(gòu)造方法,所以接下來還需要繼續(xù)尋找混合水平的高水平高強度正交表的構(gòu)造方法。

[1]楊子胥.正交表的構(gòu)造.濟南:山東人民出版社,1978.

[2]龐善起.正交表的構(gòu)造及其應用.成都:電子科技大學,2004.

[3]張應山.正交表的數(shù)據(jù)分析及其構(gòu)造.上海:華東師范大學,2006.

[4]杜蛟,王守印,溫巧燕,等.強度m的對稱正交表的遞歸構(gòu)造.應用數(shù)學學報,2012,35(2):232-244.

[5]Jiang L,Yin JX.An approach of constructing mixed-level orthogonal arrays of strength≧3.Science China(Mathematics),2013,6(56):1109-1115.

[6]Du J,Wen QY,Zhang J.New Construction of Symmetric Orthogonal Arrays of Strength t.IEICE Transactions on Fundamentals of Electronics Communications and Computer Sciences,2013,9(96):1901-1904.

[7]李學斌,余小領,李斌.運用Excel進行正交表的構(gòu)建.河南科技學院學報(自然科學版),2008,36(4):115-118.

[8]陳遠方,林曦晨,徐利華,等.常用正交表的構(gòu)造原理及SAS實現(xiàn).中國衛(wèi)生統(tǒng)計,2012,29(4):470-474.

(責任編輯：劉 壯)

國家自然科學基金(11571094)