贺深泽 (Email：heshenze@gmail.com )
整理说明
本文原稿独立作于2005年，发表于《数学的实践与认识》（Vol.39，No. 3, 99-107）。 经查，零相关与弱相关试验设计的概念与国际上一些数学家的 orthogonal Latin hypercube 及 nearly orthogonal Latin hypercube 类似。和文献相比，我的表述不大相同。例如，正交性度量引入临界值。由于特殊原因，发表时未修未校。现予修订，整理成 HTML 版本，予以公布，也为方便本博客读者。
HTML 写数学表达式很不方便。有兴趣的读者请查找原文。
摘 要
本文把选择交换向量分量的方法引入蒙特卡诺法，搜索一个与已知向量零相关的向量的速度大大增加， 构造零相关阵列成为可能。 运行数 n 在范围 4≤n≤33 ( n≠4k+2)内的 n×(n-1) 维弱相关阵列已经被构造出来，其中包含一个零相关子阵。 它们的试验点均衡分散，回归系数估计简单且相关性小。 文中给出了一个包含 7×3 维零相关子阵的7×6维弱相关矩阵及其正交化矩阵的例子。
关键词
试验设计，正交设计，零相关设计，弱相关设计，选择向量分量交换法。
Abstract
After the method of choosing to exchange vector components is introduced into the Monte Carlo method,
the speed of searching a vector which is zero-correlation with a existing vector is greatly increased,
thus the construction of zero correlation arrays is become a possible.
The n×(n-1) weak correlation arrays have been constructed, for the runs n at the range 4≤n≤33( n≠4k+2),
a zero correlation sub-array be included in them.
Their experimental points are well-balanced and the estimation of the regression coefficients are simple
and the correlation between the estimates is the weak.
An example of 7×6 weak-correlation array which included a 7×3 zero-correlation sub-array is provided in the paper.
Keywords:
Experimental designs, Orthogonal designs, Zero-correlation designs, Method of exchanging Selectively components of a vector.
一. 引 言
设一个 n×p 维设计矩阵为 X=(x1,…,xm), 其中 xi=(x1i,…,xni)T， 统计学定义其相关矩阵为 R=(rij),其中，
暂时给设计X添加一个1 向量 x0=1 做第一列，仍记为 X。统计学[2,3,5,7] 有如下的结果，对于线性模型
当X是正交的，即, 列向量的内积 (xi,xj)=0 (i≠j；i,j=1,2,…,m)，且 x-i=0, x-j=0, 则，A 是对角阵。
对原始设计的变量xi做标准化线性变换
如果X是正交的，R将是单位矩阵，A1为0 向量。如果 X 的任何两个变量之间零相关，则 R 也会是单位矩阵，A1也是0 向量。这就显示了零相关设计的重要意义。
以 (1,2,…,n)作各个变量的基本水平设计，对 n>2，田口玄一以n2个试验构造出了n2×(n+1)维设计矩阵，称为正交表[8,6]。 按代数学观点，当应用于定量因子时，它们是零相关的，并不正交。应用这类矩阵作连续量的回归设计时，当以原点为中心关于原点对称地进行水平设计时，任何两个列向量的内积(xi,xj)=0, i≠j; i,j=1,2,…,n+1，满足代数学上的正交性定义。这类表的最大特点是具有均衡分散齐整可比的特点，每个列向量的每个水平的重复次数相等；任意两个列向量a,b的水平交互配合所构成的数对(ai,bj) 具有相同的次数。这是正交设计的第二种定义方式。用于非数字和不连续量的试验，不但能够估计出每个列的效应，还可以估计出每个列的每个水平的效应。
正交表也存在一些问题，例如用作连续量的回归设计时，二水平试验的试验点密度太低，单由一组试验所能得到的结果很有限。初学者不会处理，往往把过程一律当作线性过程处理，试验次数虽然少，分析推断得不到正确的结论，很容易陷入困局。想“撒大网，捉大鱼”，结果没有捉到鱼。多水平试验试验次数太多，难以接受。
方开泰的均匀设计试图解决正交设计的这一问题 [9]。现有两类均匀设计表，由方开泰在70年代末用同余乘法计算出来的设计表(以下简称同余表)和由马长兴用蒙特卡诺法计算出来于1999年发表的设计表[10](以下简称新均匀表)。同余表的偶数表与奇数表只有最后一个试验不一样，由奇数表删除最后一个试验即得到对应的偶数表。同余表属于高相关性的设计。偶数表的秩不超过列数的一半，它们最多只能安排列数的一半那么多的试验因子。多则必包含完全相关的列向量。奇数表的高相关列数也大于列数的一半。矩阵规模越大，相关性越高。新均匀表在相关性方面有了很大的改进，但其结构不尽合理，某些设计表中还存在高相关的情况。
本文首先讨论在田口玄一正交表之外零相关设计的构造方法及其基本性质，然后构造弱相关矩阵。
二．零相关矩阵存在的条件与构造方法
2.1 零相关矩阵的存在条件
用h=(1,2,…,n)T记基本列向量。h的分量的不同排列方式有n！种，都是n维实欧氏空间Εn中的向量，且都是Εn的格点。它们构成一个向量集合S。
设一个零相关矩阵由 p 个列向量构成 xi=(x1i,…,xni)T∈S,(i=1,2,…,m),x1=h。rij= 0. rij=0 等效于叉乘和等于0,
定理：零相关矩阵存在的必要条件是 n 为大于 3 的奇数或 4 的倍数。当 n=4k+2（k为非负整数）时不存在零相关矩阵。
证明：要想(2.1)在S中有解, cf 必须是整数。因为诸 x 的分量都是非负整数，它们的内积一定是非负整数，如果cf 不是整数，(2.1) 恒不能成立。因此，当cf 不是整数时，零相关矩阵对任何 p (列数) 不存在，此时(2.1)在S中无解。
由于cf= nx-ix-j= n(1+ n) 2/4，分母为4， 若 n 为 4 的倍数，4 是 cf 的分子的一个因数，cf 必为整数。 若 n 为奇数，(n+1) 必为偶数，(1+n) 2必包含因子4，cf 的分子中必包含因子4，分母是分子的约数。
n 不能是 3，当 n=3 时，S 中只有 6 个向量，可以直接验证，任何两个的相关系数不等于零。
综上所述，零相关矩阵存在的必要条件是 n 为大于 3 的奇数或 4 的倍数。
当 n 为偶数但不是 4 的倍数时必为 2 的奇数倍数，(n+1) 必为奇数，平方后仍为奇数。此时，cf 的分子必不是 4 的整数倍数，因而必不是整数。因此，此时零相关矩阵不存在。
这个条件不充分。例如，当 n=3 时，cf(=12) 是正整数，零相关矩阵不存在。对于3< n <34, n≠4k+2, 零相关矩阵都已经被构造出来。但尚不能证明对于更大的 n 一定存在。
2.2 零相关矩阵存在的不唯一性
在 S 中与h零相关的向量应当满足条件
2.3 零相关矩阵的序贯构造方法
即使零相关矩阵存在，矩阵有多大是未知的。因此，必须序贯地逐列构造。 假定已有 p 个列为ci=(cli,…,cni)T,(i=1,2,…,m)， 其中 c1= h。待求新列 x =(x1,…,xn)T应当满足的条件是
2.3.1 全排列筛选法
1．c1= h；(整理时注：新版本(2017)的第一列均采用 h 的一个随机排列。)
2．对S的成员逐一检验，如果x∈S满足(2.2)，以x为第2列；
3．假设已找到 i 个列，如果S中的一个成员满足(2.3)，则以该列为第 i+1列。依此类推，直到不存在一个向量满足(2.3) 时，结束选取。此时称它已经零相关饱和。
随着 n 增大，计算工作量爆炸性地增加。全排列筛选法所需时间太长。
2.3.2 随机法
随机法通常叫蒙特卡诺法。随机地从 S 中选取一个向量 x，如果它与各已知列向量皆为零相关便以它作为所求的一个新列向量。 否则，随机地从 S 中另取一个向量，重复上述过程。
即使是 m=2，任选一个向量满足(2.2)其概率也非常小。当 n=7 时,一次选中的概率为184/7！≈0.0365；n=8 时,一次选取的概率大约为936/8!≈0.0232。更不用说满足(2.3)了。
2.3.3 选择向量分量交换法
引理: 设 c=(c1,…,cn)T∈S, x=(x1，…,xn)T∈S , 记 g=(c, x)。交换 x 中的两个分量，内积的增量为 c 与 x 对应分量的差之积，两对分量的差同号时，符号为负，异号时符号为正。
设 x 中被交换的分量的下标分别为 i,j(不失一般性，设 i<j),x' 为交换得到的新向量，新老内积之间的关系为
这意味着当 c, x 中的两对分量的差同号时，交换后内积减小，异号时内积增加。
由于交换 x 中任两个分量之后，a=|(c, x)-cf|或者减小，或者增大。因而可以通过选择分量进行交换，使 a 减小。 在有限次交换之后，a 可以持续减小，甚至收敛于0。随机地选择一个向量恰好 a=0 的概率非常小，经过若干次选择分量交换使之等于 0 的概率就大得多。
假设已找到 k 个列向量，第 k+1 列向量与其前 k 个列向量有 k 个 ai=|(ci, x)-cf|, 交换 x 中两个分量后，某些ai可能减小，某些ai 可能增大。选择使 ai 一致地变小收敛于 0 的分量进行交换，速度更快。
三．零相关矩阵的基本性质
3.1可堆叠
维数相同的两个零相关矩阵上下堆叠之后仍然是一个零相关矩阵。
设 A=(a1,…,ap)和 B=(b1,…,bp)是两个n×m 维零相关矩阵， 堆叠构成一个2n×m 维矩阵C=(c1,…,cm)。容易得到
维数相同的不同零相关矩阵堆叠后，试验点密度增加，均匀性互补。
（整理时注：这个命题不够准确，详见 堆叠法 )
3.2 不可拼接
饱和的零相关矩阵拼接（或嵌入）任何一个属于S的列向量都不可能再是零相关矩阵。如果存在与它的已有列向量零相关的列向量，则它就不是零相关饱和的。因此，饱和的零相关矩阵不能拼接使用。
3.3 剪裁性质
一个零相关矩阵的部分列向量所构成的矩阵仍然是一个零相关矩阵。但部分行向量所构成的矩阵不是零相关矩阵。
3.4 多态性
行(列)可移动 零相关矩阵的行（或列）向量平行移动、交换之后仍然是零相关矩阵;
列的多选性 在零相关饱和之前，列有多选性，选择不同，后面的列向量结构不同。甚至矩阵的列向量数目不同。例如，n=8 时可能只有 3 个列向量，也可以有 4 个列向量，n=9 时可能找到 4 个列向量，也可能是 5 个列向量, 取决于搜索方法、搜索深度和选取路径。
多种组合方式 利用列向量的移动、矩阵的堆叠和剪裁等性质可以组合得到多种形态的零相关矩阵。
3.5 正交性
当一个零相关矩阵所有变量的水平设计等距离且均值x-i=0,(i=1,2,…,m)，此时设计是正交的。 如果水平设计为等距离的而均值不等于0，则设计不是正交的, 但线性变换zij=xij-x-i之后为正交。 如果变量的水平设计不满足上述条件，有可能破坏零相关特性，更不满足正交性条件。
3.6 均匀性
在概率论中[2,4]有结果：n个随机变量ξi(i=1,2,…,n)在区域 D={(ai,bi),i=1,2,…,n}上的分布如果是均匀的， 每个随机变量ξi 的数学期望Eξi=(ai + bi)/ 2，方差是一个 n 阶对角矩阵，对角线上的元素依次为各随机变量的方差。 对于离散随机变量，这个命题其逆不真。但如果是随机取样，在x-i=(ai + bi)/2,(i=1,2,…,m) 的条件下，零相关设计不均匀的概率非常小。在这一意义上，零相关设计是均匀的。（整理时注：考虑到均匀的统计学特性，使用术语均衡性更恰当。
零相关矩阵正交化[3]后，列向量的长度不变，这就意味着试验区域的几何形状和体积不变，试验点充满了整个试验空间。因此，在欧氏空间的几何分布的意义上，它们也是均衡的。
四． 弱相关设计矩阵的构造方法
有两种情况引起构造弱相关设计的要求：当 n=4k+2 时不存在零相关矩阵，可以用弱相关矩阵来近似；当零相关列数太少，可以补充一些弱相关列来满足设计的需要。
弱相关矩阵可以看作为含有微弱误差的零相关矩阵，微弱的相关性可以忽略，全部回归系数的最小二乘估计相关性微弱。
设一个弱相关矩阵的 p 个列向量为ci=(c1i,…,cni)T,(i=1,2,…,p),c1= h，若x=(x1,…,xn)T∈S，与已有 p 个列向量的 p 个叉乘和的绝对值分别为
本文方法构造的弱相关矩阵各级子阵的最大相关系数(mcc)逐列递增不减。要使设计获得小的相关性，可以空闲出右边的某些列，直至只用零相关列。
五．弱相关矩阵示例
以n=7为例，弱相关矩阵 W7 (其中包含一个7×3 的零相关子阵) 示于表 1，其相关矩阵示于表 2。
有时希望构造 n×n 弱相关方阵，正交化后一定是相关的。
六．弱相关设计的分析计算方法
有条件时，用回归分析程序处理数据，可以参考各种有关回归分析的专著 [3,5,7]。
6.1 直接比较
由于设计是均衡的，把从 n 个试验的结果的直接比较中得到的较好的条件认定作为总体的一个优化条件是有优势的。 它不是推理的结果，具有现实性，靠近总体的最优值的可能性很大。在这一意义上，弱相关设计用于含有非数字和不连续量的试验也是可行的。
6.2 线性回归分析
由于变量间弱相关，回归系数估计值之间弱相关，从叉乘和矩阵即可以得到回归系数的估计值，微弱的相关性可以忽略。 回归计算时，微弱的相关性会完全消除。微小的回归系数删除之后，不必重新建模，重新计算。因此不必做逐步回归。 如果事前的设计为正交的，逐步回归分析结果与从叉乘和矩阵中得到的结果相同，还提供了回归常数的准确估计。
为演示计算过程，下面模拟一个线性过程。假设
由于叉乘和矩阵中与试验变量有关的非对角元素大多为0或接近于0，可以按表 4 的方式估计回归系数，要计算的只有那几个粗体字。
如表 4 所示，由于设计方案相关性非常微弱，接近正交，均值也提供了回归常数的准确估计。
由于x1与x3有微弱的相关性, 这两个变量的回归系数有微弱的相关性, 估计值与原值有微小的误差。
可以相信，如果采用弱相关矩阵w8(它有四个零相关列)，结果会更好。因此，对于n=4k+2 的那些弱相关设计，宁可用比它大 1 或小 1 的那个设计。
6.3 非线性回归分析
如果过程为非线性的，为了建立非线性回归模型，需要导出新变量。 导出变量通常是原试验变量的一个函数，它与原变量必然会存在一定的相关性。 应该用相关矩阵检查回归模型，如果包含相关或甚高相关的回归变量应该调整， 删除相关性太高的回归变量让位给另一些回归变量，以利提高效率。 非线性回归一定要用向后逐步回归计算，否则，不能有效地消除回归系数之间的相关性。分析推断的效果不佳。
工艺优化最适宜使用梯度法，在小扰动条件下，线性模型回归系数提供了梯度方向的良好估计。
七． 结论与存在的问题
为了说明强相关设计存在的问题，可以用如下的模拟。表5 是一个强相关设计方案（简称Q），这个设计的模板是均匀表 u*6(64)。 试验数据发生器仍为(6.1)。表 6 为试验样本的相关矩阵，表 7 为其叉乘和矩阵。因为矩阵不满秩，回归扫描出错。 表 7 表明，由叉乘和矩阵得到的四个效应的估计值(最底下一行)与真值相差很大，完全相关的x0与x3的效应估计绝对值相等，符号相反。 但容易验证，以下两个方程
y = 71.8 + 5.12x0 – 0.73x1 + 9.62x2 ------(7.1)
y = 71.8 – 0.73x1 + 9.62x2 – 5.12xx3------(7.2)
目前已经在单台微机上构造出了直到 33×32 维弱相关设计（整理注：2017 版已经更新到 n=37,而且，正交列数普遍增加）。 试验次数少，效率高。模拟展现了其易使用，计算简单，估计误差小的特点。预示其有着良好的应用前景。
现有结果，当 n>9, 零相关子阵没有超过 5 个列向量。究竟应该含有多少个零相关列向量，需要继续研究。
参 考
1．George E.P.Box, J.Stuart Hunter & William G.Hunter Statistics for Experimenters (Second Edition ) Wiley-Interscience (2004)
2．J.S.Milton, Jesse C.Arnold Introduction To Probability and Statistics (Third Edition) Irwin McGraw-Hill （1995）
3．T.Hastie , R.Tibshirani & J.Friedman The Elements of Statistical learning （Springer 2001）
4．复旦大学数学系主编 《概率论与数理统计》（第二版） 上海科学技术出版社 （1961）
5．茆诗松 丁元 周纪芗 吕乃刚编著《回归分析及其试验设计》华东师范大学出版社（1981）
6．中国科学院数学研究所统计组《常用数理统计表》科学出版社（1974）
7． R.I.Jennrich，杨自强译 逐步回归《数字计算机上用的数学方法》卷Ⅲ(4)，科学出版社（1981）
8． 田口玄一 《实验计划法》第三版 丸善株式会社 东京 （1976）
9． 方开泰 均匀设计 应用数学学报 卷3（4）（1980）
10．方开泰，马长兴 著 《正交与均匀试验设计》 科学出版社 (2001)