Kosambi–Karhunen–Loève theorem

Since the limit in the mean of jointly Gaussian random variables is jointly Gaussian, and jointly Gaussian random (centered) variables are independent if and only if they are orthogonal, we can also conclude:

Theorem. The variables Z_i have a joint Gaussian distribution and are stochastically independent if the original process {X_t}_t is Gaussian.

In the Gaussian case, since the variables Z_i are independent, we can say more:

${\displaystyle \lim _{N\to \infty }\sum _{i=1}^{N}e_{i}(t)Z_{i}(\omega )=X_{t}(\omega )}$

almost surely.

This is a consequence of the independence of the Z_k.

In the introduction, we mentioned that the truncated Karhunen–Loeve expansion was the best approximation of the original process in the sense that it reduces the total mean-square error resulting of its truncation. Because of this property, it is often said that the KL transform optimally compacts the energy.

More specifically, given any orthonormal basis {f_k} of L²([a, b]), we may decompose the process X_t as:

${\displaystyle X_{t}(\omega )=\sum _{k=1}^{\infty }A_{k}(\omega )f_{k}(t)}$

where

${\displaystyle A_{k}(\omega )=\int _{a}^{b}X_{t}(\omega )f_{k}(t)\,dt}$

and we may approximate X_t by the finite sum

${\displaystyle {\hat {X}}_{t}(\omega )=\sum _{k=1}^{N}A_{k}(\omega )f_{k}(t)}$

for some integer N.

Claim. Of all such approximations, the KL approximation is the one that minimizes the total mean square error (provided we have arranged the eigenvalues in decreasing order).

An important observation is that since the random coefficients Z_k of the KL expansion are uncorrelated, the Bienaymé formula asserts that the variance of X_t is simply the sum of the variances of the individual components of the sum:

${\displaystyle \operatorname {var} [X_{t}]=\sum _{k=0}^{\infty }e_{k}(t)^{2}\operatorname {var} [Z_{k}]=\sum _{k=1}^{\infty }\lambda _{k}e_{k}(t)^{2}}$

Integrating over [a, b] and using the orthonormality of the e_k, we obtain that the total variance of the process is:

${\displaystyle \int _{a}^{b}\operatorname {var} [X_{t}]\,dt=\sum _{k=1}^{\infty }\lambda _{k}}$

In particular, the total variance of the N-truncated approximation is

${\displaystyle \sum _{k=1}^{N}\lambda _{k}.}$

As a result, the N-truncated expansion explains

${\displaystyle {\frac {\sum _{k=1}^{N}\lambda _{k}}{\sum _{k=1}^{\infty }\lambda _{k}}}}$

of the variance; and if we are content with an approximation that explains, say, 95% of the variance, then we just have to determine an ${\displaystyle N\in \mathbb {N} }$ such that

${\displaystyle {\frac {\sum _{k=1}^{N}\lambda _{k}}{\sum _{k=1}^{\infty }\lambda _{k}}}\geq 0.95.}$

Given a representation of ${\displaystyle X_{t}=\sum _{k=1}^{\infty }W_{k}\varphi _{k}(t)}$, for some orthonormal basis ${\displaystyle \varphi _{k}(t)}$ and random ${\displaystyle W_{k}}$, we let ${\displaystyle p_{k}=\mathbb {E} [|W_{k}|^{2}]/\mathbb {E} [|X_{t}|_{L^{2}}^{2}]}$, so that ${\displaystyle \sum _{k=1}^{\infty }p_{k}=1}$. We may then define the representation entropy to be ${\displaystyle H(\{\varphi _{k}\})=-\sum _{i}p_{k}\log(p_{k})}$. Then we have ${\displaystyle H(\{\varphi _{k}\})\geq H(\{e_{k}\})}$, for all choices of ${\displaystyle \varphi _{k}}$. That is, the KL-expansion has minimal representation entropy.

Proof:

Denote the coefficients obtained for the basis ${\displaystyle e_{k}(t)}$ as ${\displaystyle p_{k}}$, and for ${\displaystyle \varphi _{k}(t)}$ as ${\displaystyle q_{k}}$.

Choose ${\displaystyle N\geq 1}$. Note that since ${\displaystyle e_{k}}$ minimizes the mean squared error, we have that

${\displaystyle \mathbb {E} \left|\sum _{k=1}^{N}Z_{k}e_{k}(t)-X_{t}\right|_{L^{2}}^{2}\leq \mathbb {E} \left|\sum _{k=1}^{N}W_{k}\varphi _{k}(t)-X_{t}\right|_{L^{2}}^{2}}$

Expanding the right hand size, we get:

${\displaystyle \mathbb {E} \left|\sum _{k=1}^{N}W_{k}\varphi _{k}(t)-X_{t}\right|_{L^{2}}^{2}=\mathbb {E} |X_{t}^{2}|_{L^{2}}+\sum _{k=1}^{N}\sum _{\ell =1}^{N}\mathbb {E} [W_{\ell }\varphi _{\ell }(t)W_{k}^{*}\varphi _{k}^{*}(t)]_{L^{2}}-\sum _{k=1}^{N}\mathbb {E} [W_{k}\varphi _{k}X_{t}^{*}]_{L^{2}}-\sum _{k=1}^{N}\mathbb {E} [X_{t}W_{k}^{*}\varphi _{k}^{*}(t)]_{L^{2}}}$

Using the orthonormality of ${\displaystyle \varphi _{k}(t)}$, and expanding ${\displaystyle X_{t}}$ in the ${\displaystyle \varphi _{k}(t)}$ basis, we get that the right hand size is equal to:

${\displaystyle \mathbb {E} [X_{t}]_{L^{2}}^{2}-\sum _{k=1}^{N}\mathbb {E} [|W_{k}|^{2}]}$

We may perform identical analysis for the ${\displaystyle e_{k}(t)}$, and so rewrite the above inequality as:

${\displaystyle {\displaystyle \mathbb {E} [X_{t}]_{L^{2}}^{2}-\sum _{k=1}^{N}\mathbb {E} [|Z_{k}|^{2}]}\leq {\displaystyle \mathbb {E} [X_{t}]_{L^{2}}^{2}-\sum _{k=1}^{N}\mathbb {E} [|W_{k}|^{2}]}}$

Subtracting the common first term, and dividing by ${\displaystyle \mathbb {E} [|X_{t}|_{L^{2}}^{2}]}$, we obtain that:

${\displaystyle \sum _{k=1}^{N}p_{k}\geq \sum _{k=1}^{N}q_{k}}$

This implies that:

${\displaystyle -\sum _{k=1}^{\infty }p_{k}\log(p_{k})\leq -\sum _{k=1}^{\infty }q_{k}\log(q_{k})}$

To further illustrate the differences between linear and non-linear approximations, we study the decomposition of a simple non-Gaussian random vector in a Karhunen–Loève basis. Processes whose realizations have a random translation are stationary. The Karhunen–Loève basis is then a Fourier basis and we study its performance. To simplify the analysis, consider a random vector Y[n] of size N that is random shift modulo N of a deterministic signal f[n] of zero mean

${\displaystyle \sum _{n=0}^{N-1}f[n]=0}$

${\displaystyle Y[n]=f[(n-p){\bmod {N}}]}$

The random shift P is uniformly distributed on [0, N − 1]:

${\displaystyle \Pr(P=p)={\frac {1}{N}},\qquad 0\leq p<N}$

Clearly

${\displaystyle \mathbf {E} \{Y[n]\}={\frac {1}{N}}\sum _{p=0}^{N-1}f[(n-p){\bmod {N}}]=0}$

and

${\displaystyle R[n,k]=\mathbf {E} \{Y[n]Y[k]\}={\frac {1}{N}}\sum _{p=0}^{N-1}f[(n-p){\bmod {N}}]f[(k-p){\bmod {N}}]={\frac {1}{N}}f\Theta {\bar {f}}[n-k],\quad {\bar {f}}[n]=f[-n]}$

Hence

${\displaystyle R[n,k]=R_{Y}[n-k],\qquad R_{Y}[k]={\frac {1}{N}}f\Theta {\bar {f}}[k]}$

Since R_Y is N periodic, Y is a circular stationary random vector. The covariance operator is a circular convolution with R_Y and is therefore diagonalized in the discrete Fourier Karhunen–Loève basis

${\displaystyle \left\{{\frac {1}{\sqrt {N}}}e^{i2\pi mn/N}\right\}_{0\leq m<N}.}$

The power spectrum is Fourier transform of R_Y:

${\displaystyle P_{Y}[m]={\hat {R}}_{Y}[m]={\frac {1}{N}}\left|{\hat {f}}[m]\right|^{2}}$

Example: Consider an extreme case where ${\displaystyle f[n]=\delta [n]-\delta [n-1]}$. A theorem stated above guarantees that the Fourier Karhunen–Loève basis produces a smaller expected approximation error than a canonical basis of Diracs ${\displaystyle \left\{g_{m}[n]=\delta [n-m]\right\}_{0\leq m<N}}$. Indeed, we do not know a priori the abscissa of the non-zero coefficients of Y, so there is no particular Dirac that is better adapted to perform the approximation. But the Fourier vectors cover the whole support of Y and thus absorb a part of the signal energy.

${\displaystyle \mathbf {E} \left\{\left|\left\langle Y[n],{\frac {1}{\sqrt {N}}}e^{i2\pi mn/N}\right\rangle \right|^{2}\right\}=P_{Y}[m]={\frac {4}{N}}\sin ^{2}\left({\frac {\pi k}{N}}\right)}$

Selecting higher frequency Fourier coefficients yields a better mean-square approximation than choosing a priori a few Dirac vectors to perform the approximation. The situation is totally different for non-linear approximations. If ${\displaystyle f[n]=\delta [n]-\delta [n-1]}$ then the discrete Fourier basis is extremely inefficient because f and hence Y have an energy that is almost uniformly spread among all Fourier vectors. In contrast, since f has only two non-zero coefficients in the Dirac basis, a non-linear approximation of Y with M ≥ 2 gives zero error.[7]

Recall that the main implication and difficulty of the KL transformation is computing the eigenvectors of the linear operator associated to the covariance function, which are given by the solutions to the integral equation written above.

Define Σ, the covariance matrix of X, as an N × N matrix whose elements are given by:

${\displaystyle \Sigma _{ij}=\mathbf {E} [X_{i}X_{j}],\qquad \forall i,j\in \{1,\ldots ,N\}}$

Rewriting the above integral equation to suit the discrete case, we observe that it turns into:

${\displaystyle \sum _{j=1}^{N}\Sigma _{ij}e_{j}=\lambda e_{i}\quad \Leftrightarrow \quad \Sigma e=\lambda e}$

where ${\displaystyle e=(e_{1}~e_{2}~\ldots ~e_{N})^{T}}$ is an N-dimensional vector.

The integral equation thus reduces to a simple matrix eigenvalue problem, which explains why the PCA has such a broad domain of applications.

Since Σ is a positive definite symmetric matrix, it possesses a set of orthonormal eigenvectors forming a basis of ${\displaystyle \mathbb {R} ^{N}}$, and we write ${\displaystyle \{\lambda _{i},\varphi _{i}\}_{i\in \{1,\ldots ,N\}}}$ this set of eigenvalues and corresponding eigenvectors, listed in decreasing values of λ_i. Let also Φ be the orthonormal matrix consisting of these eigenvectors:

${\displaystyle {\begin{aligned}\Phi &:=\left(\varphi _{1}~\varphi _{2}~\ldots ~\varphi _{N}\right)^{T}\\\Phi ^{T}\Phi &=I\end{aligned}}}$

It remains to perform the actual KL transformation, called the principal component transform in this case. Recall that the transform was found by expanding the process with respect to the basis spanned by the eigenvectors of the covariance function. In this case, we hence have:

${\displaystyle X=\sum _{i=1}^{N}\langle \varphi _{i},X\rangle \varphi _{i}=\sum _{i=1}^{N}\varphi _{i}^{T}X\varphi _{i}}$

In a more compact form, the principal component transform of X is defined by:

${\displaystyle {\begin{cases}Y=\Phi ^{T}X\\X=\Phi Y\end{cases}}}$

The i-th component of Y is ${\displaystyle Y_{i}=\varphi _{i}^{T}X}$, the projection of X on ${\displaystyle \varphi _{i}}$ and the inverse transform X = ΦY yields the expansion of X on the space spanned by the ${\displaystyle \varphi _{i}}$:

${\displaystyle X=\sum _{i=1}^{N}Y_{i}\varphi _{i}=\sum _{i=1}^{N}\langle \varphi _{i},X\rangle \varphi _{i}}$

As in the continuous case, we may reduce the dimensionality of the problem by truncating the sum at some ${\displaystyle K\in \{1,\ldots ,N\}}$ such that

${\displaystyle {\frac {\sum _{i=1}^{K}\lambda _{i}}{\sum _{i=1}^{N}\lambda _{i}}}\geq \alpha }$

where α is the explained variance threshold we wish to set.

We can also reduce the dimensionality through the use of multilevel dominant eigenvector estimation (MDEE).[8]

There are numerous equivalent characterizations of the Wiener process which is a mathematical formalization of Brownian motion. Here we regard it as the centered standard Gaussian process W_t with covariance function

${\displaystyle K_{W}(t,s)=\operatorname {cov} (W_{t},W_{s})=\min(s,t).}$

We restrict the time domain to [a, b]=[0,1] without loss of generality.

The eigenvectors of the covariance kernel are easily determined. These are

${\displaystyle e_{k}(t)={\sqrt {2}}\sin \left(\left(k-{\tfrac {1}{2}}\right)\pi t\right)}$

and the corresponding eigenvalues are

${\displaystyle \lambda _{k}={\frac {1}{(k-{\frac {1}{2}})^{2}\pi ^{2}}}.}$

This gives the following representation of the Wiener process:

Theorem. There is a sequence {Z_i}_i of independent Gaussian random variables with mean zero and variance 1 such that

${\displaystyle W_{t}={\sqrt {2}}\sum _{k=1}^{\infty }Z_{k}{\frac {\sin \left(\left(k-{\frac {1}{2}}\right)\pi t\right)}{\left(k-{\frac {1}{2}}\right)\pi }}.}$

Note that this representation is only valid for ${\displaystyle t\in [0,1].}$ On larger intervals, the increments are not independent. As stated in the theorem, convergence is in the L² norm and uniform in t.

Similarly the Brownian bridge ${\displaystyle B_{t}=W_{t}-tW_{1}}$ which is a stochastic process with covariance function

${\displaystyle K_{B}(t,s)=\min(t,s)-ts}$

can be represented as the series

${\displaystyle B_{t}=\sum _{k=1}^{\infty }Z_{k}{\frac {{\sqrt {2}}\sin(k\pi t)}{k\pi }}}$

In communication, we usually have to decide whether a signal from a noisy channel contains valuable information. The following hypothesis testing is used for detecting continuous signal s(t) from channel output X(t), N(t) is the channel noise, which is usually assumed zero mean Gaussian process with correlation function ${\displaystyle R_{N}(t,s)=E[N(t)N(s)]}$

${\displaystyle H:X(t)=N(t),}$

${\displaystyle K:X(t)=N(t)+s(t),\quad t\in (0,T)}$

When the channel noise is white, its correlation function is

${\displaystyle R_{N}(t)={\tfrac {1}{2}}N_{0}\delta (t),}$

and it has constant power spectrum density. In physically practical channel, the noise power is finite, so:

${\displaystyle S_{N}(f)={\begin{cases}{\frac {N_{0}}{2}}&|f|<w\\0&|f|>w\end{cases}}}$

Then the noise correlation function is sinc function with zeros at ${\displaystyle {\frac {n}{2\omega }},n\in \mathbf {Z} .}$ Since are uncorrelated and gaussian, they are independent. Thus we can take samples from X(t) with time spacing

${\displaystyle \Delta t={\frac {n}{2\omega }}{\text{ within }}(0,''T'').}$

Let ${\displaystyle X_{i}=X(i\,\Delta t)}$. We have a total of ${\displaystyle n={\frac {T}{\Delta t}}=T(2\omega )=2\omega T}$ i.i.d observations ${\displaystyle \{X_{1},X_{2},\ldots ,X_{n}\}}$ to develop the likelihood-ratio test. Define signal ${\displaystyle S_{i}=S(i\,\Delta t)}$, the problem becomes,

${\displaystyle H:X_{i}=N_{i},}$

${\displaystyle K:X_{i}=N_{i}+S_{i},i=1,2,\ldots ,n.}$

The log-likelihood ratio

${\displaystyle {\mathcal {L}}({\underline {x}})=\log {\frac {\sum _{i=1}^{n}(2S_{i}x_{i}-S_{i}^{2})}{2\sigma ^{2}}}\Leftrightarrow \Delta t\sum _{i=1}^{n}S_{i}x_{i}=\sum _{i=1}^{n}S(i\,\Delta t)x(i\,\Delta t)\,\Delta t\gtrless \lambda _{\cdot }2}$

As t → 0, let:

${\displaystyle G=\int _{0}^{T}S(t)x(t)\,dt.}$

Then G is the test statistics and the Neyman–Pearson optimum detector is

${\displaystyle G({\underline {x}})>G_{0}\Rightarrow K<G_{0}\Rightarrow H.}$

As G is Gaussian, we can characterize it by finding its mean and variances. Then we get

${\displaystyle H:G\sim N\left(0,{\tfrac {1}{2}}N_{0}E\right)}$

${\displaystyle K:G\sim N\left(E,{\tfrac {1}{2}}N_{0}E\right)}$

where

${\displaystyle \mathbf {E} =\int _{0}^{T}S^{2}(t)\,dt}$

is the signal energy.

The false alarm error

${\displaystyle \alpha =\int _{G_{0}}^{\infty }N\left(0,{\tfrac {1}{2}}N_{0}E\right)\,dG\Rightarrow G_{0}={\sqrt {{\tfrac {1}{2}}N_{0}E}}\Phi ^{-1}(1-\alpha )}$

And the probability of detection:

${\displaystyle \beta =\int _{G_{0}}^{\infty }N\left(E,{\tfrac {1}{2}}N_{0}E\right)\,dG=1-\Phi \left({\frac {G_{0}-E}{\sqrt {{\tfrac {1}{2}}N_{0}E}}}\right)=\Phi \left({\sqrt {\frac {2E}{N_{0}}}}-\Phi ^{-1}(1-\alpha )\right),}$

where Φ is the cdf of standard normal, or Gaussian, variable.

When N(t) is colored (correlated in time) Gaussian noise with zero mean and covariance function ${\displaystyle R_{N}(t,s)=E[N(t)N(s)],}$ we cannot sample independent discrete observations by evenly spacing the time. Instead, we can use K–L expansion to decorrelate the noise process and get independent Gaussian observation 'samples'. The K–L expansion of N(t):

${\displaystyle N(t)=\sum _{i=1}^{\infty }N_{i}\Phi _{i}(t),\quad 0<t<T,}$

where ${\displaystyle N_{i}=\int N(t)\Phi _{i}(t)\,dt}$ and the orthonormal bases ${\displaystyle \{\Phi _{i}{t}\}}$ are generated by kernel ${\displaystyle R_{N}(t,s)}$, i.e., solution to

${\displaystyle \int _{0}^{T}R_{N}(t,s)\Phi _{i}(s)\,ds=\lambda _{i}\Phi _{i}(t),\quad \operatorname {var} [N_{i}]=\lambda _{i}.}$

Do the expansion:

${\displaystyle S(t)=\sum _{i=1}^{\infty }S_{i}\Phi _{i}(t),}$

where ${\displaystyle S_{i}=\int _{0}^{T}S(t)\Phi _{i}(t)\,dt}$, then

${\displaystyle X_{i}=\int _{0}^{T}X(t)\Phi _{i}(t)\,dt=N_{i}}$

under H and ${\displaystyle N_{i}+S_{i}}$ under K. Let ${\displaystyle {\overline {X}}=\{X_{1},X_{2},\dots \}}$, we have

${\displaystyle N_{i}}$ are independent Gaussian r.v's with variance ${\displaystyle \lambda _{i}}$

under H: ${\displaystyle \{X_{i}\}}$ are independent Gaussian r.v's.

${\displaystyle f_{H}[x(t)|0<t<T]=f_{H}({\underline {x}})=\prod _{i=1}^{\infty }{\frac {1}{\sqrt {2\pi \lambda _{i}}}}\exp \left(-{\frac {x_{i}^{2}}{2\lambda _{i}}}\right)}$