Fourier Series Estimation on Periodic Domains

Fourier Regression on a Periodic Interval

Fourier series are a natural basis for regression functions whose domain is periodic: for example, a response that varies over time of day, season, or phase. A Fourier regression estimator represents the unknown function by a finite number of sine and cosine terms. Increasing the number of terms can reduce approximation error, but the additional estimated coefficients also increase variance.

We work on the unit circle, represented by $[0,1)$ with the endpoints identified. Assume the observations $(X_i,Y_i)$, $i=1,\ldots,n$, are independent and identically distributed, $X_i$ is uniform on $[0,1)$, and $Y_i=m(X_i)+\varepsilon_i$, where $E[\varepsilon_i\mid X_i]=0$ and $E[Y_i^2]<\infty$. The independence assumption here concerns the observation pairs, not just the errors. The target $m$ is square-integrable on the interval.

Definition: The real Fourier basis on $[0,1)$ is $\phi_0(x)=1$, $\phi_{k,c}(x)=\sqrt{2}\cos(2\pi kx)$, and $\phi_{k,s}(x)=\sqrt{2}\sin(2\pi kx)$, for integers $k\geq1$. These functions form an orthonormal basis of $L^2([0,1))$, with inner product $\langle g,h\rangle=\int_0^1g(x)h(x)\,dx$. Write $\beta_0=\langle m,\phi_0\rangle$, $\beta_{k,c}=\langle m,\phi_{k,c}\rangle$, and $\beta_{k,s}=\langle m,\phi_{k,s}\rangle$.

For an integer $K\geq0$, retain the constant and the sine and cosine terms through frequency $K$. The resulting dimension is $d_K=2K+1$. Uniform design makes each population coefficient an expectation of an observed response times its basis function, which motivates estimating each coefficient by a sample average.

Definition: For $j$ indexing $\phi_0$ and the Fourier functions with frequencies $1,\ldots,K$, define $\hat\beta_j=n^{-1}\sum_{i=1}^nY_i\phi_j(X_i)$. The truncated Fourier regression estimator is $\hat m_K(x)=\sum_{j}\hat\beta_j\phi_j(x)$, where the sum is over those $2K+1$ indices.

Coefficient Estimation and Integrated Risk

The coefficient estimates are unbiased: conditional mean zero gives $E[Y\mid X]=m(X)$, and uniformity of $X$ then gives $E[Y\phi_j(X)]=\int_0^1m(x)\phi_j(x)\,dx=\beta_j$. Since the observation pairs are identically distributed, each summand has the same distribution; since the pairs are independent, the summands for a fixed coefficient are independent. Thus $\operatorname{Var}(\hat\beta_j)=n^{-1}\operatorname{Var}(Y\phi_j(X))$.

The integrated squared error has an exact decomposition. Parseval's identity for this orthonormal basis gives the squared norm of a function as the sum of squared coefficients. The omitted coefficients contribute a deterministic approximation error; the retained coefficients contribute estimation error.

Theorem: Under the assumptions above, for every integer $K\geq0$, $$ E\!\left[\int_0^1\{\hat m_K(x)-m(x)\}^2\,dx\right] = \sum_{k>K}(\beta_{k,c}^2+\beta_{k,s}^2) +\frac{1}{n}\sum_{j:\,j\text{ retained}}\operatorname{Var}\{Y\phi_j(X)\}. $$ In particular, the variance term is at most $2(2K+1)E[Y^2]/n$.

Proof. Let $m_K=\sum_{j:\,j\text{ retained}}\beta_j\phi_j$ be the orthogonal projection of $m$ onto the retained span. The difference $\hat m_K-m_K$ belongs to that span, while $m-m_K$ is orthogonal to it. Therefore, $\|\hat m_K-m\|_2^2=\|\hat m_K-m_K\|_2^2+\|m-m_K\|_2^2$. Orthonormality expresses the first squared norm as $\sum_j(\hat\beta_j-\beta_j)^2$ over retained indices. Taking expectations gives $\sum_j\operatorname{Var}(\hat\beta_j)$, because each estimate is unbiased. Parseval's identity expresses the second squared norm as $\sum_{k>K}(\beta_{k,c}^2+\beta_{k,s}^2)$. Substituting the variance of each sample average proves the equality. Finally, $|\phi_j(x)|^2\leq2$ for every basis function, so $\operatorname{Var}(Y\phi_j(X))\leq E[Y^2\phi_j(X)^2]\leq2E[Y^2]$. Summing over $2K+1$ terms proves the bound. $\square$

The upper bound need not be the exact variance: the coefficient estimates can have different variances. Nor is independence between estimates of different coefficients needed for this integrated-risk identity. Orthonormality turns the squared norm into a sum of coefficientwise squares, so no cross-products between distinct estimated coefficients appear.

Worked Example: Exact Variance for a Constant Regression Function

Suppose $X_i$ are independent uniform observations, $Y_i=2+\varepsilon_i$, and the errors are independent of the designs and of one another, with $E[\varepsilon_i]=0$ and $\operatorname{Var}(\varepsilon_i)=1$. Here $m(x)=2$, so the only nonzero Fourier coefficient is $\beta_0=2$. Take $K=3$, giving seven retained coefficients.

For the constant basis function, $Y\phi_0(X)=Y$, so its variance is $1$. For each nonconstant basis function $\phi_j$, independence of $Y$ and $X$, together with $E[\phi_j(X)^2]=1$ and $E[\phi_j(X)]=0$, gives $\operatorname{Var}(Y\phi_j(X))=E[Y^2]E[\phi_j(X)^2]=E[Y^2]=2^2+1=5$. There are six such terms: sine and cosine at each of three frequencies.

The truncation error is zero because a constant is already in the retained span. The integrated risk is therefore $\{1+6(5)\}/n=31/n$. For $n=100$, this is $0.31$. The nonzero variance of the nonconstant estimated coefficients reflects random fluctuations in the empirical products $Y_i\phi_j(X_i)$, even though the true nonconstant coefficients are zero.

Coefficient Decay and Approximation Error

Smooth periodic functions tend to have small high-frequency Fourier coefficients. A useful way to state a quantitative smoothness condition is to bound their weighted squared sum. For $s>0$, suppose $\sum_{k\geq1}k^{2s}(\beta_{k,c}^2+\beta_{k,s}^2)\leq R^2$. This is a periodic Sobolev-type condition; the constant coefficient is not part of this bound because it does not affect the truncation tail.

Result: Under this coefficient condition, the squared Fourier approximation error after retaining frequencies through $K\geq0$ satisfies $$ \int_0^1\{m(x)-m_K(x)\}^2\,dx =\sum_{k>K}(\beta_{k,c}^2+\beta_{k,s}^2) \leq R^2(K+1)^{-2s}. $$

Derivation. For every omitted frequency $k>K$, $k^{2s}\geq(K+1)^{2s}$, and hence $\beta_{k,c}^2+\beta_{k,s}^2\leq(K+1)^{-2s}k^{2s}(\beta_{k,c}^2+\beta_{k,s}^2)$. Sum this inequality over $k>K$, then bound the resulting weighted tail by the full nonnegative weighted sum, which is at most $R^2$. The equality is Parseval's identity. This proves the stated bound.

Worked Example: A Decaying Fourier Series

Consider the periodic function with Fourier expansion $m(x)=\sum_{k=1}^{\infty}k^{-3/2}\cos(2\pi kx)$. Its cosine coefficient in the normalized real basis is $\beta_{k,c}=k^{-3/2}/\sqrt{2}$, and every sine coefficient and the constant coefficient are zero. The cosine series converges absolutely because $\sum_{k\geq1}k^{-3/2}<\infty$, so it defines a continuous periodic function.

For $s=1$, the weighted coefficient sum is $\sum_{k\geq1}k^2\beta_{k,c}^2=\frac12\sum_{k\geq1}k^{-1}$, which diverges. Thus this particular series does not meet the $s=1$ Sobolev condition; applying that condition to it would be an error. Its actual squared truncation error is $\frac12\sum_{k>K}k^{-3}$. Since $x^{-3}$ is decreasing, $\sum_{k=K+1}^{\infty}k^{-3}\leq\int_K^\infty x^{-3}\,dx=1/(2K^2)$ for $K\geq1$. Consequently its squared truncation error is at most $1/(4K^2)$. At $K=4$, the bound is $1/64$.

This calculation illustrates why the precise coefficient decay matters. The tail can be bounded directly from the specified coefficients, while a Sobolev bound is available only when its weighted sum is finite.

Balancing Approximation and Variance

Combining the risk identity with the Sobolev tail bound and the general variance bound gives $$ E\|\hat m_K-m\|_2^2 \leq R^2(K+1)^{-2s}+\frac{2(2K+1)E[Y^2]}{n}. $$ The first term decreases with $K$; the second increases with $K$. For large $K$, the two terms have orders $K^{-2s}$ and $K/n$. Balancing them gives $K$ of order $n^{1/(2s+1)}$, and an upper risk bound of order $n^{-2s/(2s+1)}$, with constants depending on the smoothness radius and the second moment. This is a rate calculation from the displayed bound, not a claim that the rough variance bound is always sharp.

Worked Example: Choosing a Truncation Level from a Risk Bound

Suppose a bound for a twice-smoothness calculation ($s=1$) has been simplified to the proxy $Q(K)=K^{-2}+K/n$, for integer $K\geq1$, and take $n=1000$. Treating $K$ temporarily as continuous, differentiate: $Q'(K)=-2K^{-3}+1/1000$. Setting this derivative to zero gives $K^3=2000$, or $K\approx12.60$.

Check nearby integers rather than rounding without comparison. At $K=12$, $Q(12)=1/144+12/1000\approx0.018944$. At $K=13$, $Q(13)=1/169+13/1000\approx0.018917$. At $K=14$, $Q(14)=1/196+14/1000\approx0.019102$. Thus $K=13$ has the smallest proxy among these nearby choices. Its order, $n^{1/3}$, matches the balance $K^{-2}$ against $K/n$. The proxy is a bound-based guide, not a guarantee that the actual unknown risk is minimized at exactly 13.

Assumptions and a Common Pitfall

Uniform design and an orthonormal basis make the population coefficient a simple expectation and avoid solving an empirical Gram-matrix system in this construction. If the design is not uniform, then $E[Y\phi_j(X)]$ generally does not equal the Fourier coefficient $\int_0^1m(x)\phi_j(x)\,dx$; the unbiasedness argument above no longer applies as written.

Pitfall: assuming independent errors alone justifies the variance calculation. It does not. If an error at one observation can depend on another observation's design point, the products $Y_i\phi_j(X_i)$ need not be independent, even if the errors are mutually independent and each error has conditional mean zero given its own design.
Correction: The result here assumes the observation pairs $(X_i,\varepsilon_i)$ are independent and identically distributed, in addition to the stated conditional-mean and moment conditions. Alternatively, a variance calculation would need explicit assumptions directly ensuring the relevant summands are independent and have the stated common variance. Without such conditions, covariance terms between summands may change the coefficient variance.
Pitfall: treating every coefficient tail bound as a Sobolev bound. A rate such as $R^2(K+1)^{-2s}$ requires the corresponding weighted coefficient sum to be finite. Check that condition, or derive a tail bound directly from the coefficients, as in the decaying-series example.

Key takeaway: Fourier truncation separates periodic regression risk into a frequency tail and the variance of estimated coefficients. Coefficient decay controls the tail; sample size and the number of retained terms control variance. The variance calculation depends on independence across observation pairs, not merely on marginal conditional error assumptions.

Check Your Understanding

Use the definitions and risk decomposition above to answer the following.

How many real Fourier coefficients are retained when the truncation frequency is $K=5$?
Why does $E[Y\phi_j(X)]=\beta_j$ hold under uniform design and conditional mean zero?
What coefficient condition gives a squared truncation error bounded by a constant times $K^{-2s}$?
For a constant regression function and homoskedastic independent noise, why can estimated nonconstant Fourier coefficients still have positive variance?
Why is independence of the errors, without independence of observation pairs or of the relevant summands, insufficient for the variance calculation?