Lattices and Dual Lattices in Optimal Experimental Design for Fourier Models

R.A. Bates, University of Warwick E. Riccomagno, University of Warwick R. Schwabe, Free University, Berlin H.P. Wynn, University of Warwick

Abstract Number theoretic lattices, used in integration theory, are studied from the viewpoint of the design and analysis of experiments. For certain Fourier regression models lattices are optimal as experimental designs because they produce orthogonal information matrices. This generalises early results by Hoel (1965). When the Fourier model is restricted, that is a special subset of the full factorial (cross–spectral) model is used, there is a difficult inversion problem to find generators for an optimal design for the given model. Asymptotic results are derived for certain models as the dimension of the space goes to infinity. These can be thought of as a complexity theory connecting designs and models or as special type of Nyquist sampling theory.

Keywords: experimental design, integer lattices, Nyquist sampling, complexity

1 Introduction

Integer lattices have been used extensively in numerical integration. The books by Niederreiter (1992), Conway and Sloane (1992), Sloan and Joe (1994) and Fang and Wang (1994) are basic references. Periodic, and indeed Fourier models are the natural models to be associated with the methods because of the periodic nature of the lattices themselves. In this paper we investigate lattices as experimental designs for fitting Fourier models. It will be seen that as with classical polynomial models orthogonality plays an important role and yields D–optimum designs in the sense of Kiefer and Wolfowitz. Orthogonality will be related to the so–called “dual lattice” and there is a strong connection with classical designs based on defining subgroups and aliasing.

The challenge of the theory, as we see it, is to solve “inversion” problems. That is, given a particular model find the lattice which as an experimental design allows orthogonal estimation of all terms in the model, without aliasing. The paper follows a very brief mentioning of the main idea in Bates, Buck, Riccomagno and Wynn (1996) and emphasis on the experimnetal design background in Riccomagno, Schwabe and Wynn (1995).

1.1 Fourier model

Let A^{+} be any set of d–dimensional integer vectors not containing 0 and such that h ∈ A^{+} implies −h ∈/ A^{+}. Then let A = A^{+} ∪ {−A^{+}} ∪ {0}. Let #(A^{+}) = m so that #(A) = 2m + 1. Although the results have general applicability we shall refer to the vectors h as frequencies to help motivation. The set A defines the following trigonometric regression model

E(Y (x)) = μ(x) = θ_{0} + \sqrt{2} ∑_{h∈A^{+}} [θ_{h} sin(2πh^{t}x) + ϕ_{h} cos(2πh^{t}x)]
= α_{0} + ∑_{h∈A^{+}} [β_{h}e^{2πih^{t}x} + δ_{h}e^{−2πih^{t}x}]

where θ_{h}, ϕ_{h} and β_{h}, δ_{h} are related by

α_{0} = θ_{0}
β_{h} = \frac{ϕ_{h} − iθ_{h}}{\sqrt{2}}
δ_{h} = \frac{ϕ_{h} + iθ_{h}}{\sqrt{2}}

The two forms of the model, real and complex, can be written in matrix form

μ(x) = X(x) θ = Z(x) α

where

θ = {θ_{0} : {θ_{h}}_{h∈A^{+}} : {ϕ_{h}}_{h∈A^{+}}}^{t}
X(x) = \{1 : \{\sqrt{2} sin(2πh^{t}x)\}_{h∈A^{+}} : \{\sqrt{2} cos(2πh^{t}x)\}_{h∈A^{+}}\}
α = {α_{0} : {β_{h}}_{h∈A^{+}} : {δ_{h}}_{h∈A^{+}}}^{t}
Z(x) = {1 : {e^{2πih^{t}x}}_{h∈A^{+}} : {e^{−2πih^{t}x}}_{h∈A^{+}}} .

The equation relating α and θ is

α = Hθ

where

H = ( \begin{matrix}1 & 0 & 0 \\ 0 & −\frac{i}{\sqrt{2}}I & \frac{1}{\sqrt{2}}I \\ 0 & \frac{i}{\sqrt{2}}I & \frac{1}{\sqrt{2}}I \end{matrix} )

where I is the m × m identity matrix and H is an unitary matrix: H^{∗}H = I equivalently H^{−1} = H^{∗}, where H^{∗} is the transpose and conjugate matrix of H.

1.2 Orthogonality conditions

The set of x_{i} values D = {x_{1}, …, x_{N}} is called here the experimental design and we are interested in the design–model pair (D, A). Reparametrizing the vector of the values of A in D we have for the vector of observations Y

E(Y) = Xθ = XH^{−1}Hθ = XH^{−1}α = Zα

where the ith row of X is X(x_{i}) and

Z = XH^{−1}

thus Z is the design matrix from the model in complex form. It is easily verified that the ith row of Z is Z(x_{i}) . A statistical version of the model which until now has only been expressed in expectation form is

Y (x) = μ(x) + ϵ

where ϵ are uncorrelated errors with zero mean and equal variance, σ^{2} < ∞ . Then the best linear unbiased estimator for θ, the usual least squares estimator, assuming X has full rank, is

θ̂ = (X^{t}X)^{−1} X^{t}Y

Since Z^{∗} = (XH^{−1})^{∗} = HX^{t}, the corresponding estimator for α̂ is given by

α̂ = (Z^{∗}Z)^{−1} Z^{∗}Y
= (HX^{t}XH^{−1})^{−1} HX^{t}Y
= H (X^{t}X)^{−1} X^{t}Y
= Hθ̂,

as expected from standard theory.

This paper is critically concerned with the choice of the experimental design D for a given model A. The theory of optimal experimental design can be used to show that D–optimality namely achieving max |X^{t}X| over all designs with a given sample size N is achieved for the models described here when X^{t}X = NI, namely for orthogonal designs. This is an extension of early work by Hoel (1965). See also Kiefer and Wolfowitz (1959). The D–optimality property of lattices was pointed out by Bates, Buck, Riccomagno and Wynn (1996). Riccomagno, Schwabe and Wynn (1995) present the full D–optimality theory. For the models of type A, above, namely which contain h and −h, orthogonality is equivalent to D-optimality and several other types of optimality. The full details appear in Riccomagno, Schwabe and Wynn (1995). Here we concentrate on the orthogonality.

Definition 1 A design with N points is called completely orthogonal if X^{t}X = NI where X is the design matrix.

Lemma 1 For a design model pair (D, A)

X^{t}X = NI if and only if Z^{∗}Z = NI

Proof. From the definitions Z^{∗}Z = HX^{t}XH^{−1} and X^{t}X = H^{−1}Z^{∗}ZH. □ Now Z^{∗}Z has a block structure

Z^{∗}Z = ( \begin{matrix}N & a_{1} & b_{1} \\ a_{2} & A & C \\ b_{2} & D & B \end{matrix} )

where for all h_{1}.h_{2} ∈ A^{+}

a_{1_{h_{1}}} = b_{2_{h_{1}}} = ∑_{x∈D} e^{2πjh_{1}^{t}x}
a_{2_{h_{1}}} = b_{1_{h_{1}}} = ∑_{x∈D} e^{−2πjh_{1}^{t}x}
A_{h_{1},h_{2}} = ∑_{x∈D} e^{2πj(−h_{1}+h_{2})^{t}x}
C_{h_{1},h_{2}} = ∑_{x∈D} e^{2πj(−h_{1}−h_{2})^{t}x}
D_{h_{1},h_{2}} = ∑_{x∈D} e^{2πj(+h_{1}+h_{2})^{t}x}
B_{h_{1},h_{2}} = ∑_{x∈D} e^{2πj(+h_{1}−h_{2})^{t}x}  (1)

and a_{1} = a_{2}^{∗} = b_{2}^{t} = b̄_{1}, A = B^{∗} and C = D^{∗}.

The conditions for complete orthogonality are, then, a_{1}^{t} = a_{2}^{t} = b_{1} = b_{2} = 0 and C = D = 0 and A = B = NI. On the diagonal, that is when h_{1} = h_{2}, the conditions are met since A_{h_{1},h_{1}} = B_{h_{1},h_{1}} = 1. Outside the diagonal they have to be derived from the conditions on the design. These conditions are discussed in the next section.

2 Experimental design

The designs we investigate are those given by so–called integration lattices. These are used extensively in algebraic number theory and numerical integration. One can consider this paper as an investigation of the role of lattices in functional estimatio n or approximation rather than integration.

Definition 2 An infinite lattice L_{∞} in IR^{d} is a discrete subset of IR^{d} closed under summation and subtraction. An infinite lattice is called an integration infinite lattice if it includes ZZ^{d}. We define a (finite) lattice L := L_{∞} ∩ [0, 1[^{d}.

In group theory an infinite lattice is an infinite Abelian group. A finite Abelian group, L, is also formed under the operation x_{1} ∘ x_{2} = {x_{1} + x_{2}}, where {x} is the point in [0, 1[^{d} whose coordinates are the fractional part of the coordinates of x. The order (number of elements) of a lattice will be denoted by N. We work with a particular class of integration infinite lattices formed by the so–called good lattice points or one–generator lattices, defined as follows: let g = (g_{1}, …, g_{d}) ∈ ZZ^{d} and N ∈ ZZ_{+}, with the constraint that each component of g is prime with respect to N, then

L_{∞} = \{\{\frac{kg}{N} \} ⊕ ZZ^{d} : k ∈ \{0, …, N − 1\}\}

where ⊕ denotes the direct sum. The integer vector g is called the generator of the lattice and N is the order of L_{∞} ∩ [0, 1[^{d}. We write the corresponding lattice in a more useful way as

L = \{\{\frac{kg}{N} \} : k = 0, …, N − 1\}
= \{(\frac{kg_{1} mod N}{N} , ..., \frac{kg_{d} mod N}{N} ); k = 0, ..., N − 1\}

which can be visualised as points on a line wrapped around [0, 1[^{d} interpreted as the d–dimensional torus. We shall write L_{g,N} = L for a one-generator lattice with generator g and order N.

A main observation of this paper is that the orthogonality conditions can be expressed in terms of the so–called dual lattice of the lattice L.

Definition 3 The dual lattice of a finite or infinite lattice L is defined as

L^{⊥} = {h ∈ ZZ^{d} such that h^{t}x ∈ ZZ for all x ∈ L}

Lemma 2 If L = L_{g,N} is a one–generator lattice, then

L_{g,N}^{⊥} = {h ∈ ZZ^{d} such that h^{t}g ≡ 0 (mod N)}

Proof. For a one–generator lattice h ∈ L_{g,N}^{⊥} if (and only if) h^{t}kg ≡ 0 (mod N) for all k ∈ ZZ. Obviously, the latter is equivalent to h^{t}g ≡ 0 (mod N). □

Notice that L^{⊥} is itself a lattice and in particular if h ∈ L^{⊥} then sh ∈ L^{⊥} for any integer s. In engineering terms L^{⊥} contains all higher order frequencies (harmonics) of any member.

The key to the theory is the next lemma whose proof can be found in Sloan and Joe (1994, Lemma 2.1) and in Niederreiter (1992, Lemma 5.2.1) based on elementary properties of group characters.

Lemma 3 Let L be an integration lattice with N points in the unit hypercube [0, 1)^{d}, then

\frac{1}{N} ∑_{x∈L} e^{2πih^{t}·x} = \{ \begin{matrix}1 & if h ∈ L^{⊥} \\ 0 & otherwise \end{matrix}

This lemma immediately gives conditions for complete orthogonality which we present in the following theorem. The proof is immediate from (1). Note that the complete orthogonality is more easily expressed via Z following Lemma 1.

Theorem 1 For D ≡ L, a one–generator lattice in [0, 1[^{d} with generator g, the Z–matrix satisfies Z^{∗}Z = NI, complete orthogonality, if and only if

(i) ±h_{1} ∈/ L^{⊥}

(ii) ±(h_{1} − h_{2}) ∈/ L^{⊥} with h_{1} ≠ h_{2}

(iii) ±(h_{1} + h_{2}) ∈/ L^{⊥} for all h_{1}, h_{2} frequencies in the Fourier regression model A^{+}; if and only if (i^{∗}) (h_{1} − h_{2}) ∈/ L^{⊥} for all h_{1}, h_{2} ∈ A with h_{1} ≠ h_{2}.

Note that (iii) requires in particular ±2h_{1} ∈/ L^{⊥}, which is in turn implied by h_{1} ∈/ L^{⊥} and h_{1} − h_{2} ∈/ L^{⊥} if N is odd. To avoid additional technicalities we will assume N odd. Then (iii) simplifies to h_{1} + h_{2} ∈/ L^{⊥} for h_{1} ≠ h_{2}. Or equivalently, Z^{∗}Z = NI if and only if for all h, h′ ∈ A the relation h + h′ ∈ L^{⊥} or h − h′ ∈ L^{⊥} implies h = h′ or h = −h′.

Since if L contains an element then it also contains its additive inverse. In practice, for all h_{1}, h_{2} ∈ A^{+}, one need only check

h_{1} ∈/ L^{⊥}
(h_{1} − h_{2}) ∈/ L^{⊥} (h_{1} ≠ h_{2})
(h_{1} + h_{2}) ∈/ L^{⊥}

and the last relation only for h_{1} ≠ h_{2} in case N is odd.

It is convenient to express the orthogonality as a pairwise property between different frequencies h_{1}, h_{2}.

Definition 4 A pair of frequencies (h_{1}, h_{2}) with h_{1}, h_{2} ∈/ L^{⊥} is called mutually orthogonal if h_{1}, h_{2} ∈/ L^{⊥} and ±h_{1} ± h_{2} ∈/ L^{⊥}.

Note that the above definition implies h_{1} ≠ h_{2} for orthogonal h_{1} and h_{2}. The following lemma explains the importance of the above definition. Its proof is straightforward by Lemma 3.

Lemma 4 If every pair h_{1}, h_{2} ∈ A^{+} is mutually orthogonal with respect to L^{⊥} then L is completely orthogonal for the model A = A^{+} ∪ (−A^{+}) ∪ {0}.

2.1 Cosets

There is a necessary and sufficient condition for a pair (h_{1}, h_{2}) to be mutually orthogonal which is similar to the condition of not being aliased in factorial experimental design.

It is of historic interest that the term alias occurs both in the time series/signal processing literature and in the experimental design literature. That this is true in the one–dimensional case is well–known. We believe that the present extension to lattice designs is new. For readers acquainted with the classical development of the alias table in experimental design the following development in terms of cosets will be familiar. Note however how the Fourier model expressed as A^{+} ∪ {−A^{+}} gives a double coset interpretation of aliasing. We need to recall some elementary group theory.

Let (G, +) be an Abelian group, U a subgroup of G and a be in G. The set K_{a} = K − a := {−a + u : u ∈ U}K_{a} = K − a := {−a + u : u ∈ U} is called the ath-coset of G with respect to U. More technically the equivalence relation on G defined by a ∼ b if −a + b ∈ U generates the quotient G/U whose elements, the residue classes, are the cosets of G with respect to U. The cosets form a partition of G. It turns out that the orthogonality conditions in Theorem 1 are related to pairs of cosets K_{a} ∪ K_{−a}. We denote by L_{h}^{⊥} the h–coset of L^{⊥}.

Theorem 2 The pair of frequencies h_{1} and h_{2} is mutually orthogonal if and only if there exist h ∈/ L^{⊥} and h′ ∈/ L^{⊥} such that h_{1} ∈ L_{h}^{⊥} ∪ L_{−h}^{⊥} and h_{2} ∈ L_{h′}^{⊥} ∪ L_{−h′}^{⊥} and the two double cosets are disjoint.

Proof. Let h_{1} and h_{2} be mutually orthogonal. Assume that h_{1} and h_{2} are in the same double coset L_{h}^{⊥} ∪ L_{−h}^{⊥}. Without loss of generality h_{1} ∈ L_{h}^{⊥} and, hence, h_{1}−h ∈ L^{⊥}. Then either h_{2}−h ∈ L^{⊥} or h_{2}+h ∈ L^{⊥} and, consequently, h_{1}−h_{2} ∈ L^{⊥} or h_{1} + h_{2} ∈ L^{⊥} in contradiction to the asumption which proves the direct part.

Conversely, if h_{1} and h_{2} are not mutually orthogonal, then either h_{1} − h_{2} ∈ L^{⊥} or h_{1} + h_{2} ∈ L^{⊥}. Hence, h_{1} ∈ L_{h_{2}}^{⊥} ∪ L_{−h_{2}}^{⊥} which shows that h_{1} and h_{2} are in the same double coset. □

Note that since our models always contain the constant term θ_{0} (α_{0}) we are only considering h_{1}, h_{2} ∈/ L^{⊥} since any h ∈ L^{⊥} will be aliased with the constant term. The theory can clearly be adjusted to the case with a constant term missing.

We exhibit the partition on ZZ^{d} defined by the cosets of the dual lattice of a one–generator lattice design L with N points, N odd. Since we choose N mutually prime with all the g_{i} for all i = 1, …, d, then without loss of generality we may assume g_{1} = 1. Thus the lattice is

L = \{\frac{k(1, g_{2}, …, g_{d}) (mod N)}{N} : k = 0, …, N − 1\}

and its dual is

L^{⊥} = {(a_{1}, …, a_{d}) ∈ ZZ^{d} : (a_{1}, …, a_{d})^{t}(1, …, g_{d}) ≡ 0 (mod N)}

The (h_{1}, …, h_{d})–coset is

L_{(h_{1},...,h_{d})}^{⊥} = {(a_{1} − h_{1}, …, a_{d} − h_{d}) ∈ ZZ^{d} : (a_{1}, …, a_{d}) ∈ L^{⊥}}

It turns out that there are N cosets forming a partition of ZZ^{d}. They are enumerated as

L_{(0,...,0)}^{⊥}, L_{(1,...,0)}^{⊥}, …, L_{(N−1,...,0)}^{⊥}

For notational convenience we write L_{k}^{⊥} = L_{(k,...,0)}^{⊥} for k = 0, …, N − 1 and L_{0}^{⊥} = L^{⊥}.

Theorem 3 For the lattice L, L_{(h_{1},...,h_{d})}^{⊥} = L_{k}^{⊥} if and only if h_{1}+g_{2}h_{2}+…+g_{d}h_{d} ≡ k (mod N)

Proof. Suppose h ∈ L_{k}^{⊥}. We have to prove that if h ∈ L_{(h_{1},...,h_{d})}^{⊥} then h_{1} + g_{2}h_{2} + … + g_{d}h_{d} ≡ k (mod N). Then h may be written as

\begin{aligned}h & = (a_{1} − k, a_{2}, …, a_{d}) & (a_{1}, …, a_{d}) ∈ L^{⊥} \\ & = (b_{1} − h_{1}, …, b_{d} − h_{d}) & (b_{1}, …, b_{d}) ∈ L^{⊥} \end{aligned}

This implies

(a_{1}, a_{2}, …, a_{d}) = (b_{1} − h_{1} + k, b_{2} − h_{2}, …, b_{d} − h_{d}) ∈ L^{⊥}

This condition expands to

(b_{1} − h_{1} + k) + g_{2}(b_{2} − h_{2}) + … + g_{d}(b_{d} − h_{d}) ≡ 0 (mod N)
−(h_{1} + h_{2}g_{2} + … + g_{d}h_{d}) + k ≡ 0 (mod N)

that is (h_{1} + h_{2}g_{2} + … + g_{d}h_{d}) ≡ k (mod N), which proves the assertion. Similar arguments prove that h ∈ L_{k}^{⊥}. □

Since N is odd there are \frac{N−1}{2} double cosets plus L^{⊥} itself, specifically

L_{0}^{⊥}, L_{1}^{⊥} ∪ L_{N−1}^{⊥}, …, L_{k}^{⊥} ∪ L_{N−k}^{⊥}, …, L⊥_{\frac{N−1}{2}} ∪ L⊥_{\frac{N+1}{2}} (k = 2, ..., \frac{N−3}{2})

This fact is of course clear from the underlying group theory.

Theorem 2 can be reformulated as h and h′ are mutually orthogonal if and only if

h_{1} + g_{2}h_{2} + … + g_{d}h_{d} ≠ ± (h′_{1} + g_{2}h′_{2} + … + g_{d}h′_{d}) (mod N)

Because of Lemma 1 if there are any two pairs of frequencies which are not mutually orthogonal Z (and X) becomes singular that is the model as a whole is not identifiable (a single off diagonal 1 makes Z^{t}Z singular). Thus models can have at most \frac{N−1}{2} terms and for every chosen k = 0, …, \frac{N−1}{2} at most one model–frequency satisfies the equation

h_{1k} + g_{2}h_{2k} + … + g_{d}h_{dk} ≡ k (mod N)

It is not surprising that there are at most \frac{N−1}{2} terms in A^{+} because there are at most N parameters.

As an example consider a two dimensional lattice, L, generated by (1, 5) and with N = 13 points. The design D = L is completely orthogonal for the model defined by A = (1, 0), (0, 1), (1, 1), (1, −1) because its frequencies generate respectively

L_{(1,0)}^{⊥} ∪ L_{(12,0)}^{⊥} ( 1 + 5 · 0 ≡ 1 (mod 13))
L_{(5,0)}^{⊥} ∪ L_{(8,0)}^{⊥} ( 0 + 5 ≡ 5 (mod 13))
L_{(6,0)}^{⊥} ∪ L_{(7,0)}^{⊥} ( 1 + 5 ≡ 6 (mod 13))
L_{(4,0)}^{⊥} ∪ L_{(9,0)}^{⊥} ( 1 − 5 ≡ −4 = 9 (mod 13))

It is a useful exercise to find other identifiable models for this design.

The lattice design generated by (1, 5) and N = 23 is suitable to identify the main effects for the model {(1, 0), (2, 0), (0, 1), (0, 2), (1, ±1), (1, ±2), (2, ±1), (2, ±2)}, but not the interactions since the frequencies (1, 2) and (2, 2) are aliased, being both in L_{11}^{⊥} ∪ L_{−11}^{⊥}.

3 Finding good generators: inversion

The previous section can be considered as a foundation for a more difficult problem: given a model, specified by a set of frequencies A, find the generator g = (g_{1}, …, g_{d}) and N for which the orthogonality conditions of Theorem 1 hold. In a moment we shall classify particular models of interest. It will then work out that certain models, for example without high–order cross frequencies can lead to designs (generators) for which N only grows polynomially or even linearly in the dimension d. This points to a complexity theory for the design model pair (A, D). Let us put this more precisely: what is the smallest sample size N^{∗} such that there exists a generator g^{∗} with D = L_{g^{∗},N^{∗}} satisfying the orthogonality conditions of Theorem 1? In many cases it appears that N^{∗} is hard to find explicitly and we shall find L_{g,N} with a suitable upper bound N ≥ N^{∗}.

Let A be in a particular class of models and let K(A) be some measure of the size of A such as the total number of frequencies in A, or the order of the model. Considering N^{∗} as a function of A and namely: N^{∗}(A), so that N(A) ≥ N^{∗}(A). A plot of N^{∗}(A) against K(A) gives a frontier below which no orthogonal design exists for A. Typically as K(A) increases N^{∗}(A) will increase. Thus the broad objective in terms of complexity is to find this frontier. For example in the case of no cross spectra we shall find a linear boundary. Finally as the dimension tends to infinity for a fixed class of models we can express K(A) in terms of d or present (N^{∗}, d) plots or (N, d) plots for an upper bound N ≥ N^{∗}.

Since N^{∗} represents a broad measure of the sampling density and K(A) will typically measure quantities like the largest fitted frequency or the number of distinct frequencies, the result can be thought of as a kind of Nyquist–Shannon sampling theory.

3.1 Linear models

In this section we consider models with no interactions that is cross spectral terms but including all the marginal frequencies up to a certain order. The set A for the so–referred additive or linear models is formed by points on the coordinate axis, and

A^{+} = \{ \begin{matrix}(h_{1}, 0, 0, …, 0) : & h_{1} = 1, …, m_{1} \\ (0, h_{2}, 0, …, 0) : & h_{2} = 1, …, m_{2} \\ … \\ (0, …, 0, 0, h_{d}) : & h_{d} = 1, …, m_{d} \end{matrix} \}

We use the notation F (d; m_{1}, …, m_{d}; 1) to refer to this model.

According to Theorem 3 we look for g = (g_{1}, …, g_{d}) and N such that

0, g_{k}h_{k} : h_{k} = ±1, …, ±m_{k} k = 1, …, d

are all different (mod N).

The recursive scheme

g_{1} = 1
g_{d} = m_{d−1}g_{d−1} + 1

and N = 2m_{d}g_{d}+1 gives an orthogonal design but the size N increases exponentially with the dimension d. In fact for m_{1} = … = m_{d} = m we have

g_{d} = \frac{m^{d} − 1}{m − 1} and N = 2m\frac{m^{d} − 1}{m − 1} + 1

To improve upon this arrange the model factors in such a way that m_{1} ≥ m_{2} ≥ … ≥ m_{d}. Define the recursive procedure

g_{1} = 1
g_{d} = (d − 1)m_{1} + 1

and N = 2m_{1}^{2}(d − 1) + 1 ≥ 2 max_{k=1,...,d}{g_{k}m_{k}} + 1 gives a design whose size N increases linearly with the dimension when

the orders m_{k} remain bounded. Note that a new rearrangement of the factors has to be made each time a new factor is added.

Let us comment again on the choice of N and its link to the Nyquist frequency as in Fourier Analysis. The Nyquist frequency in discrete Fourier Analysis is the highest frequency f^{∗} used in Fourier transform. If Δt is the sampling interval in time, then 2πf^{∗} = \frac{π}{Δt}. For all the frequencies f < f^{∗} the components with frequencies f^{∗} +k (k ∈ ZZ) are all aliased with each other and the sum of their contributions will appear collected, at frequency f. This is sometimes refered to as “Nyquist folding”. To avoid this the frequencies under study must be less than the Nyquist frequency, equivalently the sampling interval will be chosen to make this happen. In other terms we require 2πf < \frac{π}{Δt}. In the case of N equally spaced points Δt is \frac{1}{N} we have 2πf < πN i.e. 2f + 1 ≤ N.

In the additive case presented here the Nyquist frequency of the one dimensional case is replaced by max_{k=1,...,d}{g_{k}m_{k}}. This is also a suitable K(A) in the discussion in Section 3. We shall return to the general question of how certain aliased frequency contributions collect at other model terms in Section 4.

3.2 Two–factor interaction models

In this section we study the model F (d; 2, …, 2; 2) i.e. a two–factor interaction model with marginal frequencies equal to 0, 1 or 2. The frequency set is

A^{+} = \{ \begin{matrix}(k, 0, …, 0) & (0, k, …, 0) & … & (0, …, k, 0) & (0, …, k) \\ (k, h, …, 0) & (0, k, h, …, 0) & … & (0, …, k, h) \\ … \\ (k, 0, …, h) & (0, k, …, h) & … & & & : k, h = 1, 2 \end{matrix} \}

The lattice design generated by the sequence

g_{1} = 1
g_{d} = g_{d−1} + 4α_{d}
α_{1} = 1
α_{2^{j}} = 3α_{2^{j−1}} − 1
α_{k} = α_{2^{j}} for k = 2^{j}q, q odd

with N = 2g_{d} + 2g_{d−1} + 1 allows a model with all main effects for each factor (dimension) orthogonally to all second order interaction: a resolution IV design in experimental design. Note however that the interaction can be aliased. For example h_{1} = (1, 2, 0, …, 0) and h_{2} = (−2, 0, 1, 0, …, 0) are in the same double coset L_{11}^{⊥} since g_{1} + 2g_{2} = g_{3} − 2g_{1} = 11.

The given sequence of generators starts as

(1, 5, 13, 17, 37, 41, 53, 109, 113, 121, 125, 145, 157, 161, 325, 329)

and the α–sequence has a self–similar structure as shown by

1, 2, 1, 5, 1, 2, 1, 14, 1, 2, 1, 5, 1, 2, 1, 41,
1, 2, 1, 5, 1, 2, 1, 14, 1, 2, 1, 5, 1, 2, 1, 122
1, 2, 1, 5, 1, 2, 1, 14, 1, 2, 1, 5, 1, 2, 1, 41,
1, 2, 1, 5, 1, 2, 1, 14, 1, 2, 1, 5, 1, 2, 1, 365...

Standard methods for solving recurrence relations yield the solutions α_{2^{j}} = \frac{1}{2}(3^{j} +1) and g_{2^{j}+1} = 4 3^{j} +1. This implies that g_{d}, and N, increases like d^{γ} where γ = \frac{log 3}{log 2} ≈ 1.5858.

To study the whole parameter set of F (d; 2, …, 2; 2) we suggest the recursive sequence of generators

g_{1} = 1
g_{2} = 5
g_{d} = 4g_{d−1} + 2g_{d−2} + 1

i. e. g_{d} = \frac{1}{10}(\sqrt{6}(2+\sqrt{6})^{d}−\sqrt{6}(2−\sqrt{6})^{d}−2) ≈ \frac{\sqrt{6}}{10}(2+\sqrt{6})^{d} with N = 4(g_{d}+g_{d−1})+1 ≈ \frac{6}{5}(2 + \sqrt{6})^{d}.

Another computer–generated sequence up to dimension 14 is

(1, 5, 23, 60, 77, 173, 222, 409, 535, 634, 935, 1182, 1361, 1497)

and the corresponding sample sizes for d ≥ 2

25, 113, 333, 549, 1001, 1581, 2525, 4677, 6277, 8469, 10173, 11533.

These sample sizes are obtained by the formula 4(g_{d} + g_{d−1}) + 1 where 2(g_{d} + g_{d−1}) is the highest frequency in the model. Again we see the important Nyquist–like role played by the highest frequencies in driving up the sample size.

4 The fitted model

In the usual regression setup, assuming the original model is correct, we can estimate μ(x) at a particular x by the unbiased estimator

μ̂(x) = X(x)θ̂ = Z(x)α̂

Under the orthogonality conditions, the variance of estimation under the model A is easily seen to be

Var (μ̂(x)) = E (μ̂(x) − μ(x))^{2} = \frac{σ^{2}}{n} #(A)

where σ^{2} is the error variance. Under the Kiefer–Wolfowitz theory max Var(μ̂(x)) over [0, 1]^{d} is also minimised: G-optimality. The integrated variance is

I_{1} = ∫_{[0,1)^{d}} #(A)\frac{σ^{2}}{n} dx = #(A)\frac{σ^{2}}{n}

It is instructive to consider the consequences of taking the observations as arising as if the underlying model was given by μ(x) = f(x) for a wide choice of f(x). Thus suppose that f(x) is expressible in a standard Fourier representation

f(x) = ∑_{h∈Z Z^{d}} α(h)e^{2πih^{t}x}  (2)

where α(h) = ∫_{[0,1]^{d}} f(x)e^{−2πih^{t}x} dx. Sufficient conditions for uniform convergence are well known and omitted.

Under suitable conditions the set {e^{2πih^{t}·x} : h ∈ ZZ^{d}} forms an orthogonal basis of L^{2}(IR^{d}) and through the Fourier transform there is a one–to–one mapping from l^{2}(IR^{d}) to L^{2}(IR^{d}). We may thus fix an order on ZZ^{d} (the same as in Section 1) and write the Fourier decomposition of f(x) as the scalar product of the infinite complex vectors

Z_{∞}(x) = (1, …, e^{2πih^{t}x}, …, e^{−2πih^{t}x}, …)^{t}
α_{∞} = (α(0), …, α(h), …, α(−h), …)^{t}

and

f(x) = Z_{∞}(x)^{t}α_{∞}

Suppose we have a Fourier model A and a one–generator lattice design D = L with N points is chosen such that (D, A) is completely orthogonal. If we observe without error the observation vector is

Y = (f(x_{1}), …, f(x_{N}))^{t} = Z_{∞}α_{∞}

where the kth row of Z_{∞} is Z_{∞}(x_{k})^{t}. Then, using the orthogonality, the least squared estimator is

α̂ = (Z^{∗}Z)^{−1}Z^{∗}Y = \frac{1}{N} Z^{∗}Y = \frac{1}{N} Z^{∗}Z_{∞}α_{∞}

Then

\frac{1}{N} Z^{∗}Z_{∞} = \frac{1}{N} \{∑_{x∈L} e^{2πi(−h+k)^{t}x}\}_{h∈A,k∈Z Z^{d}}
= {δ_{h,k}}_{h∈A∪{0},k∈Z Z^{d}}

where

δ_{h,k} = \{ \begin{matrix}1 & if −h + k ∈ L^{⊥} \\ 0 & otherwise \end{matrix}

or equivalently

δ_{h,k} = \{ \begin{matrix}1 & if k ∈ L_{h}^{⊥} \\ 0 & otherwise \end{matrix}

This analysis generalises the Nyquist folding to lattices. Thus the fitted model, when there is no error, is

μ̂(x) = ∑_{h∈A} ( ∑_{k∈L_{h}^{⊥}} α(k)) e^{2πix^{t}h}

It is instructive also to give the integrated squared error. Employing the L^{2} orthogonality of the Fourier terms over [0, 1]^{d} we have

I_{2} = ∫_{[0,1]^{d}} (Ŷ (x) − f(x))^{2} dx = ∑_{h∈A} ( ∑_{k∈L_{h}^{⊥}; k≠h} α(k))^{2} + ∑_{k∈Z^{d}−A} α(k)^{2}

If there is observation error then I_{2} represents the bias term and the full integrated mean square error is

IMSE = I_{1} + I_{2}

It is clear from (4) that I_{2} is unbounded unless further conditions are added to the Fourier coefficients α(k) over and above L^{2}–conditions: ∑_{k∈Z Z} α(k)^{2} < ∞. A natural condition arising is

∑_{k∈Z Z} |α(k)| < ∞

which is equivalent to absolute convergence of f(x). Stronger conditions still are the E_{α}(c) conditions used in the integration theory of lattices (see Niederreiter, 1992, Definition 5.1).

Optimum design theory can be used to show that in addition to D–optimality for the Fourier model the orthogonal designs described here also minimize the integrated variance under the regression model. Further research is needed to find designs which guard against growth of the bias term I_{2}. In such work the relationship between the design D = L and the source of biases as given by (4) should play a major part and there is a trade off between I_{1} and I_{2}.

5 Discussion

This paper can be considered as the theoretical backgrounds for building Fourier models on integer lattices with special reference to the design model pair (D,A). Models built on space filling designs are becoming popular in building emulators, or surrogates, of large and expensive–to–run simulation codes in com puter aided design and similar areas. This body of research started with the work by Sacks et al. (1989). A recent paper on this is Bates, Riccomagno, Schwabe and Wynn (1995), having first been suggested in Bates, Buck, Riccomagno and Wynn (1996). A restriction of the Fourier models in the case where the vector x is some multidimensional set of inputs is the assumption of periodicity. For example μ(0, …, 0) = μ(1, …, 1) so that both x–values are clearly not included. One can try to avoid this restriction by enlarging, shifting, or randomising the start of the design. Also there are methods for periodising non–periodic function, see Sloan and Joe (1994). The authors’ preliminary conclusions from the practical use of these methods is that they are valuable, at least, for screening, that is detection of significant factors very much in the style of screening designs based on orthogonal arrays which have become popular in “Taguchi” methodology in Robust Engineering Design.

6 Acknowledgements

The first, second and fourth authors acknowledge the support of the UK Engineering and Physical Science Research Council. The research of the third author was supported by grants Ku719/2-1 and Ku719/2-2 of the Deutsche Forschungsgemeinschaft. The authors thank the British Council and the Deutscher Akademischer Austauschdienst for funding their joint ARC project which supports this collaboration.

References

[1] Bates, R.A., Buck, R.J., Riccomagno, E. and Wynn, H.P. Experimental design and observation for large systems. J. Royal Statistical Society, Series B, 58 (1996) 76–94.

[2] Bates, R.A., Riccomagno, E., Schwabe, R. and Wynn, H.P. Fast Fourier emulators based on space–filling lattices, Improceedings of ICIAM95.

[3] Conway, J.H. and Sloane, N.J.A. Sphere Packings, Lattices and Groups. Springer–Verlag, New York 1992.

[4] Fang, W.T. and Wang, Y. Number–theoretic Methods in Statistics. Chapman & Hall (London), 1994.

[5] Kiefer, J.C. and Wolfowitz, J. Optimum design in regression problem. Ann. Math. Statist., 30 (1959), 271–294.

[6] Hoel, P.G. Minimax distance designs in two dimensional regression. Ann. Math. Stat., 36 (1965), 1097–1106.

[7] Niederreiter, H. Random Number Generation and Quasi-Monte Carlo Methods. CBMS-NFS, SIAM (Philadelphia), 1992.

[8] Riccomagno, E., Schwabe, R. and Wynn, H.P. Lattice–based optimum design for Fourier regression. Annals of Statistics, 1995 (submitted).

[9] Sacks, J., Welch, W.J., Mitchell, T.J. and Wynn, H.P. Design and analysis of computer experiments. Statistical Science, 4 (1989), 409–435.

[10] Sloan, I.H. and Joe, S. Lattice Methods for Multiple Integration. Oxford Science Publications, Oxford, 1994.