正交表的構造方法及Matlab實現*

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

龐善起 鹿姍姍

正交表的構造方法及Matlab實現*

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

龐善起 鹿姍姍

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

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

正交表構造方法與原理

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的全部標準向量,共有

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

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

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

2.混合水平正交表

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

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

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

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

3.高強度正交表

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

(1)每一列中每個元素出現的次數相同;

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

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

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

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

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

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

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

?

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

這樣就構造出了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的向量,下文同。

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

正交表構造的Matlab實現

為了更加直觀詳細地說明Matlab程序的可行性,下面我們將選取幾個較為簡單的正交表作為例子說明利用Matlab程序實現正交表的構造方法。

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列所構成的22=4個互異的有序組,把這2 列所構成的每個有序組均換成它對應的自然數,從而生成正交表L8(4×23),結果如表2,Matlab程序與解釋如下:

表1 正交設計表L8(27)

clc;clear

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

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

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);%有序組與自然數建立一一對應

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的正交表。結果如表3,和表4,Matlab程序與解釋如下:

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

clc;clear;

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

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

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

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

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的正交表

%第三部分：正交分劃的遞歸構造

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

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,列數為N,行數為q^(N-1)的正交表。比如修改初始條件,使得q=4,t=3,N=4便可得到水平為4,強度為3,列數為4,行數為64的正交表;修改初始條件,使得q=5,t=3,N=4，便可得到水平為5,強度為3,列數為6,行數為125的正交表。

表3 正交設計表L8(24)

表4 正交設計表L8(24)

結 論

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

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

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

[3]張應山.正交表的數據分析及其構造.上海:華東師范大學,2006.

[4]杜蛟,王守印,溫巧燕,等.強度m的對稱正交表的遞歸構造.應用數學學報,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進行正交表的構建.河南科技學院學報(自然科學版),2008,36(4):115-118.

[8]陳遠方,林曦晨,徐利華,等.常用正交表的構造原理及SAS實現.中國衛生統計,2012,29(4):470-474.

(責任編輯：劉 壯)

國家自然科學基金(11571094)