Polynomial Representations of High-Dimensional Observations of Random Processes

Abstract

The paper investigates the problem of performing a correlation analysis when the number of observations is large. In such a case, it is often necessary to combine random observations to achieve dimensionality reduction of the problem. A novel class of statistical measures is obtained by approximating the Taylor expansion of a general multivariate scalar symmetric function by a univariate polynomial in the variable given as a simple sum of the original random variables. The mean value of the polynomial is then a weighted sum of statistical central sum-moments with the weights being application dependent. Computing the sum-moments is computationally efficient and amenable to mathematical analysis, provided that the distribution of the sum of random variables can be obtained. Among several auxiliary results also obtained, the first order sum-moments corresponding to sample means are used to reduce the numerical complexity of linear regression by partitioning the data into disjoint subsets. Illustrative examples provided assume the first and the second order Markov processes.

1. Introduction

The interest in developing and improving statistical methods and models is driven by the ever increasing volumes and variety of data. Making sense of data often requires one to uncover the existing patterns and relationships in data. One of the most commonly used universal tools for data exploration is evaluation of cross-correlation among different data sets and of autocorrelation of data sequence with itself. The correlation values can be utilized to carry out sensitivity analysis, to forecast future values, to visualize complex relationships among subsystems, and to evaluate other important spatio-temporal statistical properties. However, it has been well established that correlations do not imply causalization.

The correlations measure linear dependencies between pairs of random variables which is implicitly exploited in linear regression. The correlation values rely on estimated or empirical statistical moments such as sample mean and sample variance. If the mean values are removed from data, the correlations are referred to as covariances. The pairwise correlations can be represented as a fully connected graph with edges parameterized by time shifts and possibly by amplitude adjustments between the corresponding data sequences. However, the graph may become excessively complex when the number of random variables considered is large, for example, when analyzing multiple long sequences of random observations.

The problem of extending the notion of correlations and covariances to more than two random variables has been considered previously in literature. In particular, the multirelation has been defined as a measure of linearity for multiple random variables in [1]. This measure is based on a geometric analysis of linear regression. Under the assumption of multivariate Student-t distribution, the statistical significance of the multirelation coefficient is defined in terms of eigenvalues of the correlation matrix in [2]. An univariate correlation measure for multiple random variables is defined in [3] to be a sum of elements on the main diagonal of the covariance matrix. Linear regression is again utilized in [4] to define the multiple correlation coefficient. It is derived from the coefficient of determination of linear regression as a proportion of the variance in the dependent variable which is predictable from the independent variables.

The distance correlation measure between two random vectors based on the hypothesis testing of independence of random variables is proposed in [5]. Different operations on random variables including ordering, weighting, and a nonlinear transformation are assumed in [6] to define a family of Minkowski distance measures for random vectors. A maximal correlation measure for random vectors is introduced in [7]. It generalizes the concept of maximal information coefficient while nonlinear transformations are also used to allow assessment of nonlinear correlations. Similarly to [3], the sample cross-correlation between multiple realizations of two random vectors is shown in [8] to be proportional to the sum of diagonal elements of the product of the corresponding correlation matrices. The reference [9] investigates two-dimensional general and central moments for random matrices which are invariant to similarity transformations. Two new distance measures for random vectors are defined in [10] assuming the joint characteristic function. These measures can be also used to test for statistical independence of random vectors. The most recent paper [11] derives a linear correlation measure for multiple random variables using determinants of the correlation submatrices.

This brief literature survey indicates that there are still no commonly accepted correlation measures for multiple random variables. The measures considered in the literature are either rigorously derived but mathematically rather complicated, or there are many modifications of the existing, simpler, but well understood measures. In this paper, it is shown that by constraining the complexity of multivariate Taylor series to reduce the number of its parameters or degrees-of-freedom, the Taylor series can be rewritten as a finite degree univariate polynomial. The independent variable of this polynomial is a simple sum of the random variables considered. The polynomial coefficients are real-valued constants, and they are application dependent. The polynomial defines a many-to-one transformation of multiple random variables to another scalar random variable. The mean value of the polynomial represents a broad class of polynomial measures which can be used for any number of random variables. The mean value of each polynomial element corresponds to a general or central moment of the sum of random variables. Therefore, these moments are referred to here as sum-moments. In case of multiple random vectors, similarly to computing the correlation or covariance matrix by first concatenating the vectors into one long vector, the sum-moments can be readily defined and computed for such a concatenated random vector. The main advantages of assuming sum-moments to study statistical properties of multiple random variables or multiple random vectors are the clarity in understanding their statistical significance, and mathematical simplicity of their definition. Moreover, as long as the distribution of the sum of random variables can be found, the closed-form expression for the sum-moments can be obtained.

Before introducing the polynomial representations of random vectors in Section 4 together with central and general sum-moments, a number of auxiliary results are presented in Section 2 and Section 3. In particular, Section 2 summarizes the key results and concepts from the literature concerning stationary random processes, their parameter estimation via linear regression and method of moments, and how to generate the 1st and the 2nd order Markov processes is also given. Section 3 extends the results from Section 2 by deriving additional results which are used in Section 4 and in Section 5 such as a low complexity approximation of linear regression and a procedure to generate multiple Gaussian processes with defined autocorrelation and cross-correlation. The main results of the paper are obtained in Section 4 including defining a class of polynomial statistical measures and sum-moments for multiple random variables and random processes. Other related concepts involving sums of random variables are also reviewed. Section 5 provides several examples to illustrate and evaluate the results obtained in the paper. The paper is concluded in Section 6.

2. Background

This section reviews key concepts and results from the literature which are used to develop new results in the subsequent sections. Specifically, the following concepts are briefly summarized: stationarity of random processes, estimation of general and central moments and of correlation and covariance using the method of moments, definition of cosine similarity and Minkowski distance, parameter estimation via linear regression, generation of the 1st and the 2nd order Markov processes, and selected properties of polynomials and multivariate functions are also given. Note that more straightforward proofs for some lemmas are only indicated and not fully and rigorously elaborated.

2.1. Random Processes

Consider a real-valued one-dimensional random process $x\left(t\right)\in \mathcal{R}$ over a continuous time $t\in \mathcal{R}$. The process is observed at N discrete time instances $t_1<t_2<\dots <t_N$ corresponding to N random variables $X_i=x\left(t_i\right)$, $i=1,\dots ,N$. The random variables $\mathit{X}=\{X_1,\dots ,X_N\}\in {\mathcal{R}}^N$ are completely statistically described by their joint density function $f_x\left(\mathit{X}\right)$, so that $f_x\left(\mathit{X}\right)\ge 0$ and ${\int }_{{\mathcal{R}}^N}{f}_x\left(\mathit{X}\right)d\mathit{X}=1$. The process $x\left(t\right)$ is further assumed to be K-th order stationary [12].

Definition 1.

A random process $x\left(t\right)$ is K-th order stationary, if,

$$f_x(X_1,\dots ,X_K;t_1,\dots ,t_K)=f_x(X_1,\dots ,X_K;t_1+t_0,\dots ,t_K+t_0)\forall t_0\in \mathcal{R}.$$

Lemma 1.

The K-th order stationary process is also $(K-k)$-th order stationary, $k=1,2,\dots ,K-1$, for any subset $\{X_{i_1},\dots ,X_{i_{K-k}}\}\subseteq \{X_1,\dots ,X_N\}$.

Proof. 

The unwanted random variables can be integrated out from the joint density. □

Unless otherwise stated, in the sequel, all observations of random processes are assumed to be stationary.

The expectation, $E\left[X\right]={\int }_{\mathcal{R}}x{f}_x\left(x\right)dx$, of a random variable X is a measure of its mean value. A linear correlation between two random variables $x\left(t_1\right)$ and $x\left(t_2\right)$ is defined as [12],

$$R_x(t_1,t_2)=corr\left[x\left(t_1\right)x\left(t_2\right)\right]=R_x(t_2-t_1)=R_x\left(\tau \right).$$

The (auto-) covariance measures a linear dependency between two zero-mean random variables, i.e.,

$$C_x(t_1,t_2)=cov\left[x\left(t_1\right),x\left(t_2\right)\right]=E\left[\left(x\left(t_1\right)-E\left[x\left(t_1\right)\right]\right)\left(x\left(t_2\right)-E\left[x\left(t_2\right)\right]\right)\right]=C_x(t_2-t_1).$$

It can be shown that the maximum of $C_x\left(\tau \right)$, $\tau =t_2-t_1$, occurs for $\tau =0$, corresponding to the variance of a stationary process $x\left(t\right)$, i.e.,

$$C_x\left(0\right)=var\left[x\left(t\right)\right]=E\left[{\left(x\left(t\right)-E\left[x\left(t\right)\right]\right)}^2\right]\forall t\in \mathcal{R}.$$

The covariance $C_x\left(\tau \right)$ can be normalized, so that, $-1\le C_x\left(\tau \right)/C_x\left(0\right)\le 1$. Furthermore, for real-valued regular processes, the covariance has an even symmetry, i.e., $C_x\left(\tau \right)=C_x(-\tau )$,12].

For $N>2$, it is convenient to define a random vector, $\mathit{X}={[X_1,\dots ,X_N]}^T$, where ${(·)}^T$ denotes the transpose, and the corresponding mean vector, $\overline{\mathit{X}}=E\left[\mathit{X}\right]$. Then, the covariance matrix,

$${\mathit{C}}_x=E\left[(\mathit{X}-\overline{\mathit{X}}){(\mathit{X}-\overline{\mathit{X}})}^T\right]\subseteq {\mathcal{R}}^{N\times N}$$

has as its elements the pairwise covariances, ${\left[{\mathit{C}}_x\right]}_{ij}=C_x(t_j-t_i)$, $i,j=1,2,\dots ,N$.

Assume now the case of two K-th order stationary random processes $x_1\left(t\right)$ and $x_2\left(t\right)$, and the discrete time observations,

$$\begin{array}{cc}\hfill X_{1i}& =x\left(t_{1i}\right),i=1,2,\dots ,N_1\hfill \\ \hfill X_{2i}& =x\left(t_{2i}\right),i=1,2,\dots ,N_2.\hfill \end{array}$$

Using the set notation, ${\left\{X_i\right\}}_{1:K}=\{X_1,X_2,\dots ,X_K\}$, the jointly or mutually stationary processes imply time-shift invariance of their joint density function.

Definition 2.

Random processes $x_1\left(t\right)$ and $x_2\left(t\right)$ are K-th order jointly stationary, if,

$$f_x({\left\{X_{1i}\right\}}_{1:K},{\left\{X_{2i}\right\}}_{1:K};{\left\{t_{1i}\right\}}_{1:K},{\left\{t_{2i}\right\}}_{1:K})=f_x({\left\{X_{1i}\right\}}_{1:K},{\left\{X_{2i}\right\}}_{1:K};{\{t_{1i}+t_0\}}_{1:K},{\{t_{2i}+t_0\}}_{1:K})\forall t_0\in \mathcal{R}$$

is satisfied for all subsets, ${\left\{X_{1i}\right\}}_{1:K}\subseteq {\left\{X_{1i}\right\}}_{1:N_1}$, $K\le N_1$, and ${\left\{X_{2i}\right\}}_{1:K}\subseteq {\left\{X_{2i}\right\}}_{1:N_2}$, $K\le N_2$.

Lemma 2.

The K-th order joint stationarity implies the joint stationarity of all orders smaller than K.

Proof. 

The claim follows from marginalization of the joint density function to remove unwanted variables. □

The cross-covariance of random variables $X_1=x_1\left(t_1\right)$ and $X_2=x_2\left(t_2\right)$ being discrete time observations of jointly stationary random processes $x_1\left(t\right)$ and $x_2\left(t\right)$, respectively, is defined as,

$$C_{x_1{x}_2}(t_1,t_2)=cov\left[X_1,X_2\right]=cov\left[\left(x\left(t_1\right)-E\left[x\left(t_1\right)\right]\right)\left(x\left(t_2\right)-E\left[x\left(t_2\right)\right]\right)\right]=C_{x_1{x}_2}(t_2-t_1).$$

The cross-covariance can be again normalized, so it is bounded as,

$$-1\le C_{x_1{x}_2}\left(\tau \right)/\sqrt{var\left[x_1\left(t\right)\right]var\left[x_2\left(t\right)\right]}\le 1.$$

Note that, unlike for autocovariance, the maximum of $C_{x_1{x}_2}\left(\tau \right)$ can occur for any value of the argument $\tau$.

The covariance matrix for the random vectors ${\mathit{X}}_1={[X_{1i},\dots ,X_{1N_1}]}^T$ and ${\mathit{X}}_2={[X_{2i},\dots ,X_{2N_2}]}^T$ having the means, ${\overline{\mathit{X}}}_1=E\left[{\mathit{X}}_1\right]$ and ${\overline{\mathit{X}}}_2=E\left[{\mathit{X}}_2\right]$, respectively, is computed as,

$${\mathit{C}}_{x_1{x}_2}=E\left[({\mathit{X}}_1-{\overline{\mathit{X}}}_1){({\mathit{X}}_2-{\overline{\mathit{X}}}_2)}^T\right]\subseteq {\mathcal{R}}^{N_1\times N_2}.$$

Its elements are the covariances, ${\left[{\mathit{C}}_{x_1{x}_2}\right]}_{ij}=C_{x_1{x}_2}(t_{1i}-t_{2j})$, $i=1,\dots ,N_1$, $j=1,\dots ,N_2$.

In addition to the first order (mean value) and the second order (covariance) statistical moments, higher order statistics of a random variable X are given by the general and the central moments defined, respectively, as [13],

$$g_m\left(X\right)=E\left[{\left|X\right|}^m\right]and{\mu }_m\left(X\right)=E\left[|X-E\left[X\right]{|}^m\right],m=1,2,\dots$$

(1)

where $\left|X\right|$ denotes the absolute value of scalar variable X. The positive integer-valued moments (1) facilitate mathematically tractable integration, and prevent producing complex numbers, if the argument is negative. Note also that the absolute value in (1) is necessary if X is complex valued, or if m is odd, in order to make the moments to be real-valued and convex. The central moment can be normalized by the variance as,

$${\mu }_m\left(X\right)=E\left[|X-E\left[X\right]{|}^m\right]/E{\left[{(X-E\left[X\right])}^2\right]}^{m/2}.$$

The cosine similarity between two equal-length vectors, ${\mathit{X}}_1={[X_{11},\dots ,X_{1N}]}^T$ and ${\mathit{X}}_2={[X_{21},\dots ,X_{2N}]}^T$, is defined as [14],

$$S_{\cos }({\mathit{X}}_1,{\mathit{X}}_2)=\frac{{\sum }_{i=1}^N{X}_{1i}{X}_{2i}}{\sqrt{{\sum }_{i=1}^N{X}_{1i}^2}\sqrt{{\sum }_{i=1}^N{X}_{2i}^2}}.$$

The Minkowski distance between two equal-length vectors ${\mathit{X}}_1$ and ${\mathit{X}}_2$ is defined as the $l_m$-norm [6], i.e.,

$$S_{\mathrm{mnk}}({\mathit{X}}_1,{\mathit{X}}_2)={∥{\mathit{X}}_1-{\mathit{X}}_2∥}_m=(\sum _{i=1}^N{\left|X_{1i}-X_{2i}{|}^m\right)}^{1/m},m=1,2,\dots$$

2.2. Estimation Methods

Statistical moments can be empirically estimated from measured data using sample moments. Such inference strategy is referred to as the method of moments [15]. In particular, under the ergodicity assumption, a natural estimator of the mean value of a random variable X from its measurements $X_i$ is the sample moment,

$$\overline{X}=E\left[X\right]\approx \frac{1}{N}\sum _{i=1}^N{X}_i.$$

The sample mean estimator is unbiased and consistent [15]. More generally, the sample mean estimator of the first moment of a transformed random variable $h\left(X\right)$ is,

$$E\left[h\left(X\right)\right]\approx \frac{1}{N-d}\sum _{i=1}^{N}h\left(X_i\right)$$

where d is the number of degrees of freedom used in the transformation h, i.e., the number of other parameters which must be estimated. For example, the variance of X is estimated as,

$$var\left[X\right]=E\left[|X-\overline{X}{|}^2\right]\approx \frac{1}{N-1}\sum _{i=1}^N{(X_i-\widehat{\overline{X}})}^2$$

where $\widehat{\overline{X}}$ is the estimate of the mean value $\overline{X}$ of X.

Assuming random sequences ${\left\{X_{1i}\right\}}_{1:N_1}$ and ${\left\{X_{2i}\right\}}_{1:N_2}$, their autocovariance and cross-covariance, respectively, are estimated as [15],

$$\begin{array}{cc}\hfill C_x\left(k\right)=E\left[X_i{X}_{i+k}\right]& \approx \frac{1}{N-k}\sum _{i=1}^{N-k}{X}_i{X}_{i+k},k\ll N.\hfill \\ \hfill C_{x_1{x}_2}\left(k\right)=E\left[X_{1i}{X}_{2(i+k)}\right]& \approx \frac{1}{N-k}\sum _{i=1}^{N-k}{X}_{1i}{X}_{2(i+k)},k\ll N=min(N_1,N_2).\hfill \end{array}$$

Since these estimators are consistent, the condition $k\ll N$ is necessary to combine a sufficient number of samples and achieve an acceptable estimation accuracy.

The parameters of a random process can be estimated by fitting a suitable data model to the measurements $X_i$. Denote such data model as, $E\left[X_i\right]=M_i\left(\mathit{P}\right)$, $i=1,2,\dots ,N$. Assuming the least-squares (LS) criterion, the vector of parameters, $\mathit{P}={[P_1,\dots ,P_D]}^T\subseteq {\mathcal{R}}^D$, is estimated as,

$$\widehat{\mathit{P}}=\underset{\mathit{P}}{argmin}\sum _{i=1}^N{(X_i-M_i\left(\mathit{P}\right))}^2.$$

For continuous parameters, the minimum is obtained using the derivatives, i.e., let $\frac{d}{d\mathit{P}}{M}_i\left(\mathit{P}\right)={\dot{M}}_i\left(\mathit{P}\right)$, and,

$$\frac{d}{d\mathit{P}}\sum _{i=1}^N{(X_i-M_i\left(\mathit{P}\right))}^2\stackrel{!}{=}0\iff \sum _{i=1}^N{\dot{M}}_i\left(\widehat{\mathit{P}}\right)M_i\left(\widehat{\mathit{P}}\right)=\sum _{i=1}^N{\dot{M}}_i\left(\widehat{\mathit{P}}\right)X_i.$$

For a linear data model, $M_i\left(\mathit{P}\right)={\mathit{w}}_i^T\mathit{P}$ where ${\mathit{w}}_i={[w_{1i},\dots ,w_{Di}]}^T$ are known coefficients, the LS estimate can be obtained in the closed-form, i.e.,

$$\widehat{\mathit{P}}={\left(\sum _{i=1}^N{\mathit{w}}_i{\mathit{w}}_i^T\right)}^{-1}\sum _{i=1}^N{\mathit{w}}_i{X}_i$$

(2)

where ${(·)}^{-1}$ denotes the matrix inverse.

2.3. Generating Random Processes

The task is to generate a discrete time stationary random process with a given probability density and a given autocovariance. The usual strategy is to generate a correlated Gaussian process followed by a nonlinear memoryless transformation. For instance, the autoregressive (AR) process described by the second order difference equation with constant coefficients $a_1$, $a_2$, and $b>0$ [12], i.e.,

$$x\left(n\right)+a_{1}x(n-1)+a_{2}x(n-2)=bu\left(n\right)$$

generates the 2nd order Markov process from a zero-mean white (i.e., uncorrelated) process $u\left(n\right)$. For $a_1=-2\sqrt{\frac{1-a}{1+a}}$ and $a_2=\frac{1-a}{1+a}$, $0<a<1$, this process has the autocovariance,

$$\begin{array}{cc}\hfill C_{2\mathrm{MP}}\left(k\right)=& \underset{{\sigma }^2}{\underbrace{b^2\frac{{(1+a)}^2}{4a^3}}}\underset{\doteq {(1-a)}^{\left|k\right|}}{\underbrace{{\left(\frac{1-a}{1+a}\right)}^{\left|k\right|/2}}}(1+a|k\left|\right)\hfill \\ \hfill \doteq & {\sigma }^2{(1-a)}^{\left|k\right|}(1+a|k\left|\right)={\sigma }^2{e}^{-\left|k\right|(-ln(1-a\left)\right)}(1+a|k\left|\right).\hfill \end{array}$$

(3)

On the other hand, for $a_1=a>0$, and $a_2=0$, the AR process,

$$x\left(n\right)+ax(n-1)=bu\left(n\right)$$

generates the 1st order Markov process with the autocovariance,

$$C_{1\mathrm{MP}}\left(k\right)=\underset{{\sigma }^2}{\underbrace{\frac{b^2}{1-a^2}}}{a}^{\left|k\right|}={\sigma }^2{e}^{-\alpha \left|k\right|},a=e^{-\alpha }.$$

(4)

Lemma 3.

[16] The stationary random process $x\left(k\right)$ with autocovariance $C_x\left(k\right)$ is transformed by a linear time-invariant system with real-valued impulse response $h\left(k\right)$ into another stationary random process $y\left(k\right)={\sum }_{i=-\infty }^{\infty }h\left(i\right)x(k-i)\equiv h\left(k\right)⊛x\left(k\right)$ with autocovariance, $C_y\left(k\right)=h\left(k\right)⊛h(-k)⊛C_x\left(k\right)$. The symbol, ⊛, denotes (discrete time) convolution.

Proof. 

By definition, the output covariance, $C_y(k,l)=E\left[(y\left(k\right)-E\left[y\left(k\right)\right])(y\left(l\right)-E\left[y\left(l\right)\right])\right]$. Substituting $y\left(k\right)={\sum }_{i=-\infty }^{\infty }h\left(i\right)x(k-i)$, and rearranging, the covariance, $C_y(k-l)={\sum }_i{\sum }_{m}h\left(m\right)C_x(i-m)h(k-i)=h\left(k\right)⊛h(-k)⊛C_x\left(k\right)$. □

Lemma 4.

A stationary random process at the output of a linear or nonlinear time-invariant system remains stationary.

Proof. 

For any multivariate function $h\left(\mathit{x}\right)\in \mathcal{R}$ and any $t_0\in \mathcal{R}$, the expectation,

$$E\left[h\left(\mathit{x}\right)\right]={\int }_{{\mathcal{R}}^N}h\left(\mathit{x}\right)f_X(\mathit{x};t_1,\dots ,t_N)d\mathit{x}={\int }_{{\mathcal{R}}^N}h\left(\mathit{x}\right)f_X(\mathit{x};t_1-t_0,\dots ,t_N-t_0)d\mathit{x}$$

assuming Definition 1 and provided that dimension N of observations is at most equal to the order of stationarity K. □

For shorter sequences, the linear transformation, $\mathit{X}=\mathit{T}\mathit{U}$, can be used to generate a normally distributed vector $\mathit{X}\in {\mathcal{R}}^N$ having the covariance matrix, ${\mathit{C}}_x=\mathit{T}{\mathit{T}}^T$, from uncorrelated Gaussian vector $\mathit{U}\in {\mathcal{R}}^N$. The mean, $E\left[\mathit{X}\right]=\mathit{T}E\left[\mathit{U}\right]$. For longer sequences, a linear time-invariant filter can be equivalently used as indicated by Lemma 3.

2.4. Polynomials and Multivariate Functions

Lemma 5.

A univariate twice differentiable function $p\left(x\right)$ is convex, if and only if, $\frac{d^2}{d{x}^2}p\left(x\right)=\ddot{p}\left(x\right)>0$ for $\forall x\in \mathcal{R}$. More generally, a twice differentiable multivariate function $f:{\mathcal{R}}^N\mapsto \mathcal{R}$ is convex, if and only if, the domain of f is convex, and its Hessian is positive semidefinite, i.e., ${\nabla }^{2}f\ge 0$ for $\forall x\in domf$.

Proof. 

See [17] (Sec. 3.1). □

Consequently, convex polynomials can be generated as follows.

Lemma 6.

Let $\mathit{Q}\in {\mathcal{R}}^{m\times m}$ be a positive semidefinite matrix, and assume a polynomial, $\ddot{p}\left(x\right)={\sum }_{i=0}^{2m-2}{p}_i{x}^i$ for $\forall x\in \mathcal{R}$ where $p_i={\sum }_{k+l=i}{\mathit{Q}}_{kl}$. Then, for any $q_0,q_1\in \mathcal{R}$, the polynomial $p\left(x\right)$ of degree $2m$,

$$p\left(x\right)=q_0+q_{1}x+\sum _{i=0}^{2m-2}\frac{p_i}{(i+1)(i+2)}{x}^{2+i}$$

is convex.

Proof. 

Let $\mathit{x}={[x^0,x^1,\dots ,x^{m-1}]}^T$. Then, ${\mathit{x}}^T\mathit{Q}\mathit{x}={\sum }_{i=0}^{2m-2}{p}_i{x}^i=\ddot{p}\left(x\right)\ge 0$ for $\forall \mathit{x}$, since $\mathit{Q}$ is positive semidefinite. Using Lemma 5 concludes the proof. □

For a non-negative integer $m\in \left\{0\right\}\cup {\mathbb{N}}_{+}=\{0,1,2,\dots \}$, assume the following notations to simplify the subsequent mathematical expressions [18]:

$$\begin{array}{ccc}\mathit{n}& =& \{n_1,n_2,\dots ,n_N\}\in {\{{\mathbb{N}}_{+}\cup 0\}}^N\hfill \\ \mathit{x}& =& \{x_1,x_2,\dots ,x_N\}\in {\mathcal{R}}^N\hfill \\ {∥\mathit{x}∥}_p& =& {\left(x_1^p+x_2^p+\cdots +x_N^p\right)}^{1/p}\hfill \\ {∥\mathit{x}∥}_1& =& |x_1|+|x_2|+\cdots +|x_N|\hfill \\ {\left|\mathit{x}\right|}_1& =& x_1+x_2+\cdots +x_N\hfill \\ \mathit{h}& =& \{h_1,h_2,\dots ,h_N\}\in {\mathcal{R}}^N\hfill \\ {\mathit{h}}^{\mathit{n}}& =& h_1^{n_1}{h}_2^{n_2}\cdots h_N^{n_N}\hfill \\ {\partial }^{\mathit{n}}f\left(\mathit{x}\right)& =& {\partial }_1^{n_1}{\partial }_2^{n_2}\cdots {\partial }_N^{n_N}f\left(\mathit{x}\right)=\frac{{\partial }^{{\left|\mathit{n}\right|}_1}}{\partial x_1^{n_1}{x}_2^{n_2}\cdots x_N^{n_N}}f\left(\mathit{x}\right)\hfill \\ m!& =& \prod _{i=1}^{m}i\hfill \\ \mathit{n}!& =& n_1!n_2!\cdots n_N!\hfill \end{array}$$

Note that ${\left|\mathit{x}\right|}_1$ denotes the sum of elements of $\mathit{x}$ whereas ${∥\mathit{x}∥}_1$ is the sum of absolute values of its elements.

Lemma 7.

The m-th power of a finite sum can be expanded as [19],

$${\left|\mathit{x}\right|}_1^m=\sum _{{\left|\mathit{n}\right|}_1=m}\frac{m!}{\mathit{n}!}{\mathit{x}}^{\mathit{n}}.$$

Proof. 

See [18]. □

Theorem 1.

The multivariate Taylor’s expansion of a $(m+1)$-order differentiable function $f:{\mathcal{R}}^N\mapsto \mathcal{R}$ about the point $\mathit{x}\in {\mathcal{R}}^N$ is written as [19],

$$f(\mathit{x}+\mathit{h})=f\left(\mathit{x}\right)+\sum _{{\left|\mathit{n}\right|}_1\le m}\frac{{\partial }^{\mathit{n}}f\left(\mathit{x}\right)}{\mathit{n}!}{\mathit{h}}^{\mathit{n}}+\sum _{{\left|\mathit{n}\right|}_1=m+1}\frac{{\partial }^{\mathit{n}}f(\mathit{x}+t\mathit{h})}{\mathit{n}!}{\mathit{h}}^{\mathit{n}}$$

for some $t\in (0,1)$.

Proof. 

See [18,19]. □

Definition 3.

A multivariate function $f\left(\mathit{x}\right)=f(x_1,\dots ,x_N)$ is said to be symmetric, if, for any permutation of its arguments $\mathit{x}$ denoted as ${\mathit{x}}^{\prime }$, $f\left(\mathit{x}\right)=f\left({\mathit{x}}^{\prime }\right)$.

3. Background Extensions

In this section, additional results are obtained which are used in the next section to introduce statistical sum-moments for random vectors. In particular, the mean cosine similarity, the mean Minkowski distance as well as the higher central moments for random vectors are defined. A polynomial approximation of univariate functions is shown to be a linear regression problem. A numerically efficient solution of the LS problem is derived. Finally, a procedure to generate multiple Gaussian processes with defined autocovariance and cross-covariance is devised.

Recall the second moment of a random variable X, i.e.,

$${\mu }_2=E\left[{(X-c)}^2\right],c\in \mathcal{R}.$$

It is straightforward to show that ${\mu }_2$ is minimized for $c=E\left[X\right]$, giving the variance of X. On the other hand, ${\mu }_2=0$, if and only if $c=E\left[X\right]\pm \sqrt{-var\left[X\right]}$.

For random vectors, both cosine similarity and Minkowski distance are random variables. Assuming that the vectors are jointly stationary, and their elements are identically distributed, the mean cosine similarity can be defined as,

$${\overline{S}}_{\cos }({\mathit{X}}_1,{\mathit{X}}_2)=\frac{{\sum }_{i=1}^{N}E\left[X_{1i}{X}_{2i}\right]}{\sqrt{{\sum }_{i=1}^{N}E\left[X_{1i}^2\right]-{\overline{\mathit{X}}}_{1i}^2}\sqrt{{\sum }_{i=1}^{N}E\left[X_{2i}^2\right]-{\overline{\mathit{X}}}_{2i}^2}}=\frac{1}{N}\sum _{i=1}^N\frac{E\left[X_{1i}{X}_{2i}\right]}{{\sigma }_1{\sigma }_2}=\frac{1}{N}\sum _{i=1}^N{\rho }_{1i,2i}$$

where ${\sigma }_1^2=var\left[X_{1i}\right]$, ${\sigma }_2^2=var\left[X_{2i}\right]$, and, ${\rho }_{1i,2i}$ denotes the Pearson correlation coefficient of the i-th elements of the vectors ${\mathit{X}}_1$ and ${\mathit{X}}_2$. It should be noted that this definition of the mean cosine similarity does not account for other correlations, $E\left[X_{1i}{X}_{2j}\right]$, $i\ne j$.

The mean Minkowski distance for random vectors can be defined as,

$${\overline{S}}_{\mathrm{mnk}}({\mathit{X}}_1,{\mathit{X}}_2)={\left(\sum _{i=1}^{N}E\left[{\left|X_{1i}-X_{2i}\right|}^m\right]\right)}^{1/m},m=1,2,\dots$$

(5)

Recognizing the m-th general moment in (5), the m-th power of the mean Minkowski distance can be normalized as,

$${\tilde{S}}_{\mathrm{mnk}}^m({\mathit{X}}_1,{\mathit{X}}_2)=\sum _{i=1}^N\frac{E\left[{\left|X_{1i}-X_{2i}\right|}^m\right]}{N{\left(E\left[{\left|X_{1i}-X_{2i}\right|}^2\right]\right)}^{m/2}}={\overline{\mu }}_m({\mathit{X}}_1-{\mathit{X}}_2),m=1,2,\dots$$

where the average Minkowski distance between two random vectors is,

$${\overline{\mu }}_m({\mathit{X}}_1-{\mathit{X}}_2)=\frac{1}{N}\sum _{i=1}^N{\mu }_m(X_{1i}-X_{2i}),m=1,2,\dots$$

Furthermore, note that for $m=2$,

$$\frac{1}{2}\left(E\left[{∥{\mathit{X}}_1-{\mathit{X}}_2∥}_2^2\right]+E\left[{∥{\mathit{X}}_1+{\mathit{X}}_2∥}_2^2\right]\right)=E\left[{∥{\mathit{X}}_1∥}_2^2\right]+E\left[{∥{\mathit{X}}_2∥}_2^2\right].$$

Assuming positive integers $\mathit{m}=\{m_1,m_2,\dots ,m_N\}$, the higher order joint central moments of random vector $\mathit{X}={\left\{X_i\right\}}_{1:N}$ can be defined as,

$${\mu }_{m_1,\dots ,m_N}(X_1,\dots ,X_N)=E\left[\prod _{i=1}^N{(X_i-{\overline{X}}_i)}^{m_i}\right]$$

(6)

or, using more compact index notation as,

$${\mu }_{\mathit{m}}\left(\mathit{X}\right)=E\left[{(\mathit{X}-\overline{\mathit{X}})}^{\mathit{m}}\right].$$

3.1. Linear LS Estimation

The linear LS estimation can be used to fit a degree $(D-1)$ polynomial to N samples of a random process $x\left(t\right)$ at discrete time instances $t_1,t_2,\dots ,t_N$. Hence, consider the polynomial data model,

$$E\left[x\left(t\right)\right]\approx M(t;\mathit{P})=\sum _{k=1}^D{P}_k{t}^{k-1}.$$

Denoting ${\mathit{w}}_i={[t_i^0,t_i^1,\dots ,t_i^{D-1}]}^T$, the linear LS solution (2) gives the estimates,

$$\widehat{\mathit{P}}={\left(\sum _{i=1}^N\left[\begin{array}{cccc}{t}_i^0& t_i^1& \cdots & t_i^{D-1}\\ t_i^1& t_i^2& \cdots & t_i^D\\ ⋮& ⋮& \ddots & ⋮\\ t_i^{D-1}& t_i^D& \cdots & t_i^{2D-2}\end{array}\right]\right)}^{-1}\sum _{i=1}^{N}x\left(t_i\right)\left[\begin{array}{c}{t}_i^0\\ t_i^1\\ ⋮\\ t_i^{D-1}\end{array}\right].$$

Assuming $D=2$ parameters, i.e., the linear LS regression for a straight line, the parameters $P_1$ and $P_2$ to be estimated must satisfy the following equality:

$$\frac{d}{d{P}_1}\sum _{i=1}^N{\left(X_i-w_{1i}{P}_1-w_{2i}{P}_2\right)}^2=-2\sum _{i=1}^N{w}_{1i}{X}_i+2\sum _{i=1}^N{w}_{1i}^2{P}_1+2\sum _{i=1}^N{w}_{1i}{w}_{2i}{P}_2\stackrel{!}{=}0.$$

Denoting the weighted averages, $\overline{X}={\sum }_{i=1}^N{w}_{1i}{X}_i$, ${∥{\mathit{w}}_1∥}_2^2={\sum }_{i=1}^N{w}_{1i}^2$, and ${\overline{w}}_{12}={\sum }_{i=1}^N{w}_{1i}{w}_{2i}$, a necessary but not sufficient condition for the linear LS estimation of parameters $P_1$ and $P_2$ is,

$$\overline{X}={∥{\mathit{w}}_1∥}_2^2{P}_1+{\overline{w}}_{12}{P}_2.$$

(7)

In the LS terminology, the values ${∥{\mathit{w}}_1∥}_2^2$ and ${\overline{w}}_{12}$ represent independent variables whereas $\overline{X}$ is a dependent variable.

Note that all N measurements are used in (7). However, if N is sufficiently large, the data points could be split into two disjoint sets of $N_1$ and $N_2$ elements, respectively, $N_1+N_2=N$. For convenience, denote the sums,

$$\begin{array}{ccc}{\overline{X}}_1=\frac{1}{a_1}\sum _{i\in {\mathcal{I}}_1}{w}_{1i}{X}_i,& {\overline{W}}_{11}=\frac{1}{a_1}\sum _{i\in {\mathcal{I}}_1}{w}_{1i}^2,& {\overline{W}}_{12}=\frac{1}{a_1}\sum _{i\in {\mathcal{I}}_1}{w}_{1i}{w}_{2i},\\ {\overline{X}}_2=\frac{1}{a_2}\sum _{i\in {\mathcal{I}}_2}{w}_{1i}{X}_i,& {\overline{W}}_{21}=\frac{1}{a_2}\sum _{i\in {\mathcal{I}}_2}{w}_{1i}^2,& {\overline{W}}_{22}=\frac{1}{a_2}\sum _{i\in {\mathcal{I}}_2}{w}_{1i}{w}_{2i},\end{array}$$

where $a_1,a_2>1$ are some constants (to be determined later), and ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$ are two disjoint index sets, such that ${\mathcal{I}}_1\cup {\mathcal{I}}_2=\{1,2,\dots ,N\}$, and the cardinality, $|{\mathcal{I}}_1|=N_1$ and $|{\mathcal{I}}_2|=N_2$. Note also that,

$$\begin{array}{ccc}{∥{\mathit{w}}_1∥}_2^2& =& a_1{\overline{W}}_{11}+a_2{\overline{W}}_{21}\\ {\overline{w}}_{12}& =& a_1{\overline{W}}_{12}+a_2{\overline{W}}_{22}\\ \overline{X}& =& a_1{\overline{X}}_1+a_2{\overline{X}}_2.\end{array}$$

(8)

Consequently, the approximate LS estimates of parameters $P_1$ and $P_2$ are readily computed as,

$$\left[\begin{array}{c}{\widehat{P}}_1\\ {\widehat{P}}_2\end{array}\right]={\left[\begin{array}{cc}{\overline{W}}_{11}& {\overline{W}}_{12}\\ {\overline{W}}_{21}& {\overline{W}}_{22}\end{array}\right]}^{-1}\left[\begin{array}{c}{\overline{X}}_1\\ {\overline{X}}_2\end{array}\right].$$

(9)

There are $2^N$ possibilities how to split N data points into two disjoint subsets indexed by ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$. More importantly, the estimates (9) do not guarantee the minimum LS fit, i.e., achieving the minimum squared error (MSE). However, the complexity of performing the LS fit is greatly reduced by splitting the data, since independently of the value of $N\gg 1$, only a $2\times 2$ matrix needs to be inverted. The optimum LS fit and the approximate LS fit (9) are depicted in Figure 1. The points $A_1$ and $A_2$ in Figure 1 correspond to the data subsets indexed by ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$, respectively. The mid-point, $B=a_1{A}_1+a_2{A}_2$, follows from (8). Note that B is always located at the intersection of the optimum and the approximate lines $L{S}_{\mathrm{opt}}$ and $L{S}_{\mathrm{apr}}$, respectively. The vertical arrows at points $A_1$ and $A_2$ in Figure 1 indicate that the dependent values ${\overline{X}}_1$ and ${\overline{X}}_2$ are random variables.

The larger the variation of the gradient of the line $L{S}_{\mathrm{apr}}$ in Figure 1, the larger the uncertainty and the probability that the line $L{S}_{\mathrm{apr}}$ deviates from the optimum regression line $L{S}_{\mathrm{opt}}$. Since the line $L{S}_{\mathrm{apr}}$ is defined by points ${\overline{X}}_1$ and ${\overline{X}}_2$, and always, $B\in L{S}_{\mathrm{apr}}$, the spread of random variables ${\overline{X}}_1$ and ${\overline{X}}_2$ about their mean values affect the likelihood that $L{S}_{\mathrm{apr}}$ deviates from $L{S}_{\mathrm{opt}}$. In particular, given $0<p<1$, there exists $\xi >0$, such that the probability of the gradient $G_{\mathrm{apr}}$ of $L{S}_{\mathrm{apr}}$ to be within the given bounds is at least,

$$Pr\left(b_{\mathrm{l}}\left(\xi \right)<G_{\mathrm{apr}}<b_{\mathrm{u}}\left(\xi \right)\right)\ge p$$

(10)

where

$$\begin{array}{ccc}\hfill b_{\mathrm{l}}\left(\xi \right)& =& \frac{\left(E\left[{\overline{X}}_2\right]-\xi \sqrt{var\left[{\overline{X}}_2\right]}\right)-\left(E\left[{\overline{X}}_1\right]+\xi \sqrt{var\left[{\overline{X}}_1\right]}\right)}{{\overline{W}}_{22}-{\overline{W}}_{12}}\hfill \\ \hfill b_{\mathrm{u}}\left(\xi \right)& =& \frac{\left(E\left[{\overline{X}}_2\right]+\xi \sqrt{var\left[{\overline{X}}_2\right]}\right)-\left(E\left[{\overline{X}}_1\right]-\xi \sqrt{var\left[{\overline{X}}_1\right]}\right)}{{\overline{W}}_{22}-{\overline{W}}_{12}}.\hfill \end{array}$$

For stationary measurements (cf. (9)), the means and the variances in (10) are equal to,

$$\begin{array}{cc}E\left[{\overline{X}}_1\right]=E\left[\overline{X}\right]N_1/\left(a_{1}N\right)\hfill & var\left[{\overline{X}}_1\right]=var\left[\overline{X}\right]N_1/\left(a_1^{2}N\right)\hfill \\ E\left[{\overline{X}}_2\right]=E\left[\overline{X}\right]N_2/\left(a_{2}N\right)\hfill & var\left[{\overline{X}}_2\right]=var\left[\overline{X}\right]N_2/\left(a_2^{2}N\right).\hfill \end{array}$$

The uncertainty in computing $G_{\mathrm{apr}}$ from data, and thus, also the probability of line $L{S}_{\mathrm{apr}}$ deviating from line $L{S}_{\mathrm{opt}}$ is inversely proportional to the width of the interval in (10), i.e.,

$$b_{\mathrm{u}}\left(\xi \right)-b_{\mathrm{l}}\left(\xi \right)=\frac{2\xi (\sqrt{var\left[{\overline{X}}_1\right]}+\sqrt{var\left[{\overline{X}}_2\right]})}{{\overline{W}}_{22}-{\overline{W}}_{12}}\propto \frac{\frac{\sqrt{N_1}}{a_1}+\frac{\sqrt{N_2}}{a_2}}{{\overline{W}}_{22}-{\overline{W}}_{12}}.$$

(11)

Consequently, the numerator of (11) must be minimized, and the denominator maximized.

In order to minimize the numerator in (11), it is convenient to choose, $a_1=N_1^{\nu }$ and $a_2=N_2^{\nu }$, where $\nu \in {\mathcal{R}}^{+}$ is a constant to be optimized. It is straightforward to show that the expression, $(N_1^{1/2-\nu }+N_2^{1/2-\nu })$ is convex, i.e., it has a unique global minimum for $\forall \nu >1/2$. Hence, a necessary condition to reduce the approximation error is that $a_1>\sqrt{N_1}$ and $a_2>\sqrt{N_2}$. For numerical convenience, let $a_1=N_1$ and $a_2=N_2$. Then, the dependent and independent variables assumed in (9) become arithmetic averages, and the optimum index subsets have the cardinality,

$$\left\{\begin{array}{cc}|{\mathcal{I}}_1|=|{\mathcal{I}}_2|=N/2& N-even\\ |{\mathcal{I}}_1|=|{\mathcal{I}}_2|\pm 1=(N\pm 1)/2& N-odd.\end{array}\right.$$

(12)

In order to maximize the denominator in (11), assume that the independent variables $({∥{\mathit{w}}_1∥}_2^2,{\overline{w}}_{12})$ in (9) are sorted by $w_{li}$, i.e., let $w_{11}<w_{12}<\dots <w_{1N}$. This ordering and the condition (12) suggests that the disjoint index sets ${\mathcal{I}}_1$ and ${\mathcal{I}}_2$ maximizing the difference, ${\overline{W}}_{22}-{\overline{W}}_{12}$, are,

$$\left\{\begin{array}{cc}{\mathcal{I}}_1=\{1,2,\dots ,N/2\},{\mathcal{I}}_2=\{N/2+1,\dots ,N\}& N-\mathrm{even}\\ {\mathcal{I}}_1=\{1,2,\dots ,(N-1)/2\},{\mathcal{I}}_2=\{(N+1)/2,\dots ,N\}& N-\mathrm{odd}\\ \mathrm{or},{\mathcal{I}}_1=\{1,2,\dots ,(N+1)/2\},{\mathcal{I}}_2=\{(N+3)/2,\dots ,N\}.\end{array}\right.$$

Such a partitioning corresponds to splitting the data into two equal (N-even) or approximately equal (N-odd) sized subsets by the median (2-quantile) index point.

In summary, the approximate linear LS regression can be efficiently computed with a good accuracy by splitting the data into multiple disjoint subsets, calculating the average data points in each of these subsets, and then solving the set of linear equations with the same (or smaller) number of unknown parameters. The data splitting should exploit ordering of data points by one of the independent variables. It can be expected that the accuracy of approximate LS regression is going to improve with the number of data points N. Numerical evaluation of the approximate LS regression is considered in Section 5.

3.2. Generating Pairwise-Correlated Gaussian Processes

How to generate a single correlated Gaussian process is well established in the literature, and it has been described in Section 2.3. Moreover, the linear transformation to generate correlated Gaussian variables from uncorrelated Gaussian variables does not have to be square. A sufficient condition on the rank of linear transformation to obtain a positive definite covariance matrix is given by the following lemma.

Lemma 8.

The matrix, ${\mathit{T}}^T\mathit{T}$, is positive definite, provided that the matrix, $\mathit{T}\in {\mathcal{R}}^{N_1\times N_2}$, has rank $N_2$.

Proof. 

The rank $N_2$ of $\mathit{T}$ implies that $\mathit{T}$ consists of $N_2$ linearly independent columns and that $N_1\ge N_2$. The matrix ${\mathit{T}}^T\mathit{T}$ is positive definite, provided that ${\mathit{U}}^T{\mathit{T}}^T\mathit{T}\mathit{U}={∥\mathit{T}\mathit{U}∥}_2^2>0$ for $\forall \mathit{U}\in {\mathcal{R}}^{N_2}$ where ${∥·∥}_2$ denotes the Euclidean norm of a vector. Since the columns of $\mathit{T}$ are linearly independent, ${∥\mathit{T}\mathit{U}∥}_2=0$, if and only if $\mathit{U}=\mathit{0}$. □

Corollary 1.

The matrix, $\mathit{T}{\mathit{T}}^T$, is positive definite, provided that the matrix, $\mathit{T}\in {\mathcal{R}}^{N_1\times N_2}$, has rank $N_1$.

Corollary 2.

A rank $N_1$ linear transformation $\mathit{T}$ of uncorrelated Gaussian vector $\mathit{U}\in {\mathcal{R}}^{N_2}$ generates $N_1\le N_2$ correlated Gaussian variables having the (positive definite) covariance matrix, $\mathit{T}{\mathit{T}}^T$.

Furthermore, it is often necessary to generate multiple mutually correlated Gaussian processes with given autocorrelation as well as cross-correlation.

Lemma 9.

The linear transformation,

$$\left[\begin{array}{c}{\mathit{x}}_1\\ {\mathit{x}}_2\end{array}\right]=\left[\begin{array}{cc}{\mathit{T}}_1& \mathit{K}\\ \mathit{0}& {\mathit{T}}_2\end{array}\right]\left[\begin{array}{c}{\mathit{u}}_1\\ {\mathit{u}}_2\end{array}\right]$$

(13)

generates a pair of correlated Gaussian vectors ${\mathit{x}}_1\in {\mathcal{R}}^{N_1}$ and ${\mathit{x}}_2\in {\mathcal{R}}^{N_2}$ from uncorrelated zero-mean Gaussian vectors ${\mathit{u}}_1,{\mathit{u}}_2\in {\mathcal{R}}^N$ where $\mathit{0}$ denotes a zero matrix, and according to Corollary 2, it is necessary that, $max(N_1,N_2)\le N$. The corresponding (auto-) correlation and cross-correlation matrices are,

$$E\left[{\mathit{u}}_1{\mathit{u}}_1^T\right]={\sigma }_1^2\mathit{I},E\left[{\mathit{u}}_2{\mathit{u}}_2^T\right]={\sigma }_2^2\mathit{I},E\left[{\mathit{u}}_1{\mathit{u}}_2^T\right]=\mathit{0}$$

$${\mathit{C}}_{x_1}=E\left[{\mathit{x}}_1{\mathit{x}}_1^T\right]={\mathit{T}}_1{\mathit{T}}_1^T+\mathit{K}{\mathit{K}}^T,{\mathit{C}}_{x_2}=E\left[{\mathit{x}}_2{\mathit{x}}_2^T\right]={\mathit{T}}_2{\mathit{T}}_2^T,{\mathit{C}}_{x_1{x}_2}=E\left[{\mathit{x}}_1{\mathit{x}}_2^T\right]=\mathit{K}{\mathit{T}}_2^T={\mathit{C}}_{x_2{x}_1}^T$$

where $\mathit{I}$ denotes an identity matrix.

Proof. 

The proof is straightforward by substituting (13) to definitions of ${\mathit{C}}_{x_1}$, ${\mathit{C}}_{x_2}$ and ${\mathit{C}}_{x_1{x}_2}$. □

The following corollary details the procedure described in Lemma 9.

Corollary 3.

Given ${\mathit{T}}_2$, calculate the (auto-) correlation matrix, ${\mathit{C}}_{x_2}=\mathit{K}{\mathit{T}}_2^T$, or vice versa. Then, given the cross-correlation matrix, ${\mathit{C}}_{x_1{x}_2}$, compute $\mathit{K}={\mathit{C}}_{x_1{x}_2}{\mathit{T}}_2{\mathit{C}}_{x_2}^{-T}$. Finally, given ${\mathit{T}}_1$, calculate the (auto-) correlation matrix, ${\mathit{C}}_{x_1}={\mathit{T}}_1{\mathit{T}}_1^T+\mathit{K}{\mathit{K}}^T$, or, obtain ${\mathit{T}}_1$ by solving the matrix equation, ${\mathit{T}}_1{\mathit{T}}_1^T={\mathit{C}}_{x_1}-\mathit{K}{\mathit{K}}^T$. Note that the matrix equation, ${\mathit{T}}^T\mathit{T}=\mathit{C}$, can be solved for $\mathit{T}$ using the singular value decomposition, $\mathit{C}=\mathit{U}\Lambda {\mathit{U}}^T$ where $\mathit{U}$ is a unitary matrix and Λ is a diagonal matrix of eigenvalues. Then, $\mathit{T}=\mathit{U}\sqrt{\Lambda }$.

4. Polynomial Statistics and Sum-Moments for Vectors of Random Variables

The main objective of this section is to define a universal function to effectively measure the statistics of random vectors and random processes observed at multiple discrete time instances. The measure function should: (1) be universally applicable for an arbitrary number of random variables and random vectors, (2) be symmetric, so that all random variables are considered equally, (3) lead to mathematically tractable expressions, and (4) be convex to allow defining convex optimization problems.

Let $f:{\mathcal{R}}^N\mapsto \mathcal{R}$ denote such a mapping or transformation of N random variables $X_i$ to a scalar random variable Y, i.e.,

$$Y=f(X_1,X_2,\dots ,X_N)=f\left(\mathit{x}\right).$$

(14)

In order to satisfy the symmetry requirement, the random variables $X_i$ can be first combined as,

$$Y=f(X_1\circ X_2\circ \cdots \circ X_N)$$

using a binary commutative operator, ∘, such as addition or multiplication. In case of addition, it is important that the function f is nonlinear. The nonlinearity can be also used to limit the extreme values of the combined variables.

For a random process $x\left(t\right)$, the random variables are defined as, $X_i=x\left(t_i\right)$. Define a vector of discrete time instances, $\mathit{t}=(t_1,\dots ,t_N)$, and assume the index notation, $x\left(\mathit{t}\right)=\{X_1=x\left(t_1\right),\dots ,X_N=x\left(t_N\right)\}$. Then the mapping (14) can be rewritten as,

$$Y=f\left(x\left(\mathit{t}\right)\right)=\mathcal{F}\left(\mathit{t}\right)=\mathcal{F}(t_1,t_2,\dots ,t_N).$$

The mean value is the most important statistical property of the random variable Y. In addition, if the process $x\left(t\right)$ is K-th order stationary, the dimension of the mean mapping for N observations is reduced by one.

Lemma 10.

The mean of $Y=f\left(x\right(\mathit{t}\left)\right)$ for N discrete time observations of a K-th order stationary random process has dimension $(N-1)$, i.e., if $N\le K$,

$$\overline{Y}=E\left[f\left(x\right(\mathit{t})\right]=\mathcal{C}(t_2-t_1,t_3-t_1,\dots ,t_N-t_1)=\mathcal{C}({\tau }_1,{\tau }_2,\dots ,{\tau }_{N-1})=\mathcal{C}\left(\mathit{\tau }\right)$$

where ${\tau }_i=t_{i+1}-t_1$, $i=1,2,\dots ,N-1$.

Proof. 

Assuming Lemma 1 and the first sample $X_1=x\left(t_1\right)$ to be a reference, the joint probability density of N process observations becomes, $f_X(\mathit{x};0,t_2-t_1,\dots ,t_N-t_1)\equiv {\tilde{f}}_X(\mathit{x};t_2-t_1,\dots ,t_N-t_1)$, so the corresponding statistical moments have dimension $(N-1)$. □

In optimization problems, it is useful to consider the gradient of $\overline{Y}$, i.e.,

$$\nabla \overline{Y}=\left(\frac{\partial \overline{Y}}{\partial {\tau }_1},\frac{\partial \overline{Y}}{\partial {\tau }_2},\dots ,\frac{\partial \overline{Y}}{\partial {\tau }_{N-1}}\right)=\frac{\partial \overline{Y}}{\partial \mathcal{T}}\left(\frac{\partial \mathcal{T}}{\partial {\tau }_1},\frac{\partial \mathcal{T}}{\partial {\tau }_2},\dots ,\frac{\partial \mathcal{T}}{\partial {\tau }_{N-1}}\right)$$

where $\mathcal{T}$ is a application dependent measure of the vector $\mathit{\tau }$ such as the norms, $\mathcal{T}={∥\mathit{\tau }∥}_1=|{\tau }_1|+\cdots +|{\tau }_{N-1}|$, or $\mathcal{T}={∥\mathit{\tau }∥}_{\infty }=max({\tau }_1,\dots ,{\tau }_{N-1})$.

In general, assuming the Taylor’s series defined in Theorem 1 on p. 8, a multivariate function of a random vector $\mathit{x}$ can be expanded about the mean $\overline{\mathit{x}}=E\left[\mathit{x}\right]$ as,

$$f(\overline{\mathit{x}}+\mathit{h})\approx f\left(\overline{\mathit{x}}\right)+\sum _{l=1}^m\sum _{{\left|\mathit{n}\right|}_1=l}\frac{{\partial }^{\mathit{n}}f\left(\mathit{x}\right)}{\mathit{n}!}{\mathit{h}}^{\mathit{n}}.$$

(15)

Thus, the value of $f(\overline{\mathit{x}}+\mathit{h})$ is a weighted sum of ${\mathit{h}}^{\mathit{n}}$ plus an offset $f\left(\overline{\mathit{x}}\right)$. More importantly, if the partial derivatives ${\partial }^{\mathit{n}}f\left(\mathit{x}\right)$ are replaced with the coefficients $(l!p_l)$ which are independent of $\mathit{n}$, the number of parameters in (15) is greatly reduced. Moreover, instead of precisely determining the values of $p_l\in \mathcal{R}$ to obtain the best possible approximation of the original function f, it is useful as well as sufficient to constrain the Taylor expansion (15) to the class of functions that are exactly constructed as,

$$f(\overline{\mathit{x}}+\mathit{h})=f\left(\overline{\mathit{x}}\right)+\sum _{l=1}^m{p}_l\sum _{{\left|\mathit{n}\right|}_1=l}\frac{l!}{\mathit{n}!}{\mathit{h}}^{\mathit{n}}=\sum _{l=0}^m{p}_l{(h_1+\cdots +h_N)}^l$$

(16)

where $p_0=f\left(\overline{\mathit{x}}\right)$. The function expansion (16) represents a m-th degree polynomial in variable ${\left|\mathit{h}\right|}_1$. The coefficients $p_l$ of this polynomial can be set using Lemma 6, so the polynomial is convex.

The key realization is that the polynomial functions (16) have all the desired properties specified at the beginning of this section.

Claim 1.

A multivariate function $f:{\mathcal{R}}^N\mapsto \mathcal{R}$ having the desirable properties to measure the statistics of a random vector $\mathit{x}$ is the m-th degree polynomial,

$$Y=f\left(\mathit{x}\right)=\sum _{l=0}^m{p}_l{\left|\sum _{i=1}^N(X_i-E\left[X_i\right])\right|}^l$$

(17)

where $p_0=E\left[\mathit{x}\right]$, and the values of m and of the coefficients $p_l\in \mathcal{R}$ are determined by the application requirements.

The polynomial Function (17) can be used for any number of observations N. It is symmetric, so all observations are treated equally. Moreover, it is prone to integration and employing the expectation operator. The convexity can be achieved by Lemma 6.

4.1. Related Concepts

Assume the scalar function f defined in (17) for a K-th order stationarity random process $x\left(t\right)$. Define the auxiliary random variable,

$$Z\left(a\right)=\frac{1}{a}\sum _{i=1}^N(X_i-{\overline{X}}_i)$$

where a is a normalization constant, $a\ne 0$. The expression (17) can be then rewritten as, $Y\left(a\right)={\sum }_{l=0}^m{p}_l{Z}^l\left(a\right)$. For $a=1$, $Z\left(1\right)$ has a zero mean, and the variance, $var\left[Z\left(1\right)\right]={\sum }_{i,j}E\left[(X_i-{\overline{X}}_i)(X_j-{\overline{X}}_j)\right]$. For $a=N$, $Z\left(N\right)$ represents a sample mean, and its variance is equal to, $var\left[Z\left(N\right)\right]=var\left[Z\left(1\right)\right]/N^2$. For $a=\sqrt{N}$, the variance of $Z\left(\sqrt{N}\right)$ is normalized by the number of dimensions, i.e., $var\left[Z\left(\sqrt{N}\right)\right]=var\left[Z\left(1\right)\right]/N$.

For $p_l=s^l/l!$, in the limit of large m, Equation (17) gives,

$$Y=\underset{m\to \infty }{lim}\sum _{l=0}^m\frac{s^l}{l!}{\left(Z\left(a\right)\right)}^l=e^{sZ\left(a\right)}.$$

The mean, $E\left[Y\right]=E\left[e^{sZ\left(a\right)}\right]$, is the moment generating function of the random variable $Z\left(a\right)$.

In data processing, the sample mean is intended to be an estimate of the true population mean, i.e., $N\gg 1$ is required. Here, $Z\left(a\right)$ is calculated over a finite number of vector or process dimensions, so it is a random variable for $\forall a\in \mathcal{R}\backslash \left\{0\right\}$. The variable $Z\left(N\right)$ should be then referred to as an arithmetic average or a center of gravity of the random vector $\mathit{X}$ in the Euclidean space ${\mathcal{R}}^N$, i.e.,

$$Z\left(N\right)\triangleq \overline{X}=\frac{1}{N}\sum _{i=1}^N{X}_i\in \mathcal{R}.$$

(18)

Note that (18) is not an $l_1$-norm, since the variables $X_i$ are not summed with absolute values.

If the random variables $X_i$ are independent, the distribution of $Z\left(a\right)$ is given by convolution of their marginal distributions. For correlated observations, if the characteristic function, $\tilde{f}(s/a)=E\left[e^{j·{\left|\mathit{X}\right|}_{1}s/a}\right]$, of the sum ${\left|\mathit{X}\right|}_1$ can be obtained, the distribution of $Z\left(a\right)={\left|\mathit{X}\right|}_1/a$ is computed as the inverse transform,

$$f_Z\left(Z\left(a\right)\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{-jsZ\left(a\right)}\tilde{f}\left(\frac{s}{a}\right)ds.$$

Many other properties involving the sums of random variables can be obtained. For instance, if the random variables $X_i$ are independent and have zero mean, then,

$$E\left[Z^m\left(1\right)\right]-E\left[{(Z\left(1\right)-X_N)}^m\right]=\left\{\begin{array}{cc}E\left[X_N^m\right]& m=2,3\\ E\left[X_N^4\right]+6\sum _{i=1}^{N-1}E\left[X_N^2\right]E\left[X_i^2\right]& m=4\\ E\left[X_N^5\right]+10\sum _{i=1}^{N-1}E\left[X_N^2\right]E\left[X_i^3\right]+E\left[X_N^3\right]E\left[X_i^2\right]& m=5.\end{array}\right.$$

2 Considering Claim 1, an important statistic for a random vector can be defined as the m-th central sum-moment.

Definition 4.

The m-th central sum-moment of random vector $\mathit{X}\in {\mathcal{R}}^N$ is computed as,

$$\sum {\mu }_m\left(\mathit{X}\right)={\mu }_m{\left(\right|\mathit{X}|}_1)=E\left[{\left|\sum _{i=1}^N(X_i-{\overline{X}}_i)\right|}^m\right],m=1,2,\dots$$

Lemma 11.

The second central sum-moment of random vector is equal to the sum of all elements of its covariance matrix, i.e.,

$$\sum {\mu }_2\left(\mathit{X}\right)={\mu }_2\left(\right|\mathit{X}-\overline{\mathit{X}}{|}_1)=E\left[{\left(\sum _{i=1}^N(X_i-{\overline{X}}_i)\right)}^2\right]=\sum _{i,j=1}^{N}cov\left[X_i,X_j\right]={\left|{\mathit{C}}_x\right|}_1.$$

Furthermore, the second central sum-moment is also equal to the variance of $\left|\mathit{X}\right|$, i.e.,

$$\sum {\mu }_2\left(\mathit{X}\right)=var\left[{\left|\mathit{X}\right|}_1\right]=var\left[\sum _{i=1}^N(X_i-{\overline{X}}_i)\right].$$

Proof. 

The first equality is shown by expanding the expectation, and substituting for elements of the covariance matrix, ${\left[{\mathit{C}}_x\right]}_{i,j}=cov\left[X_i,X_j\right]$. The second expression follows by noting that ${\sum }_{i=1}^N(X_i-{\overline{X}}_i)$ has zero mean. □

In the literature, there are other measures involving sums of random variables. For instance, in Mean-Field Theory, the model dimensionality is reduced by representing N-dimensional vectors by their center of gravity [20]. The central point of a vector is also used in the first order approximation of multivariate functions in [21] and in the model overall variance in [22].

In Measure Theory [23], the total variation (TV) of a real-valued univariate function, $x:(t_0,t_N)\mapsto \mathcal{R}$, is defined as the supremum over all possible partitions $\mathcal{P}:t_0<t_1<\cdots <t_N$ of the interval $(t_0,t_N)$, i.e.,

$$\mathrm{TV}\left(x\right)=\underset{\mathcal{P}}{sup}\sum _{i=0}^{N-1}\left|x\right(t_{i+1}-x\left(t_i\right)|.$$

The TV concept can be adopted for observations $X_i=x\left(t_i\right)$ of a stationary random process $x\left(t\right)$ at discrete time instances, ${\left\{t_i\right\}}_{0:N}$. A mathematically tractable mean TV measure can be defined as,

$${\overline{\mathrm{TV}}}^2\left(x\right)=\underset{\mathcal{P}}{sup}E\left[\sum _{i=0}^{N-1}{|X_{i+1}-X_i|}^2\right]=\underset{\mathcal{P}}{sup}2N(E\left[X_i^2\right]-cov\left[X_{i+1},X_i\right]).$$

Jensen’s inequality for a random vector assuming equal weights can be stated as [17],

$$E\left[{\left|\frac{1}{N}\sum _{i=1}^N(X_i-{\overline{X}}_i)\right|}^m\right]\le \frac{1}{N}\sum _{i=1}^{N}E\left[|X_i-{\overline{X}}_i{|}^m\right].$$

Alternatively, exchanging the expectation and summation, Jensen’s inequality becomes,

$$\sum _{i=1}^N|E\left[X_i-{\overline{X}}_i\right]{|}^m\le E\left[\sum _{i=1}^N{|X_i-{\overline{X}}_i|}^m\right].$$

(19)

Furthermore, if the right-hand side of (19) is to be minimized and $m=2$, the inequality in (19) changes to equality. In particular, consider the minimum mean square error (MMSE) estimation of a vector of random parameters $\mathit{P}={\left\{P_i\right\}}_{1:N}$ from measurements $\mathit{X}$. Denoting ${\overline{P}}_i\left(\mathit{X}\right)=E\left[P_i|\mathit{X}\right]$, conditioned on $\mathit{X}$, the MMSE estimator $\widehat{\mathit{P}}\left(\mathit{X}\right)$ minimizes [15],

$$\begin{array}{ccc}\hfill \underset{\widehat{P}|\mathit{X}}{min}E\left[\sum _{i=1}^N{({\widehat{P}}_i\left(\mathit{X}\right)-P_i)}^2\right]& =& \underset{\widehat{P}|\mathit{X}}{min}E\left[\sum _{i=1}^N{\left(({\widehat{P}}_i\left(\mathit{X}\right)-{\overline{P}}_i\left(\mathit{X}\right))-(P_i-{\overline{P}}_i\left(\mathit{X}\right))\right)}^2\right]\hfill \\ & =& \underset{\widehat{P}|\mathit{X}}{min}\left\{\sum _{i=1}^N{({\widehat{P}}_i\left(\mathit{X}\right)-{\overline{P}}_i\left(\mathit{X}\right))}^2+E\left[\sum _{i=1}^n{(P_i-{\overline{P}}_i\left(\mathit{X}\right))}^2\right]\right\}\hfill \\ & =& \underset{\widehat{P}|\mathit{X}}{min}\sum _{i=1}^N{(E\left[{\widehat{P}}_i\left(\mathit{X}\right)-P_i|\mathit{X}\right])}^2=\underset{\widehat{P}|\mathit{X}}{min}\sum _{i=1}^N{({\widehat{P}}_i\left(\mathit{X}\right)-E\left[P_i|\mathit{X}\right])}^2\hfill \end{array}$$

(20)

where the expectations are over the conditional distribution $f_{\mathit{P}|\mathit{X}}$.

In signal processing, a length N moving average (MA) filter transforms the input sequence $X_i$ into an output sequence $Y_i$ by discrete-time convolution ⊛, i.e.,

$$Y_i=\sum _{j=0}^{N-1}{X}_{i-j}=\left[\underset{{\mathit{1}}_N}{\underbrace{11\dots 1}}\right]⊛X_i.$$

The (auto-) correlations of the input and output sequences are related by Lemma 3 on p. 6, i.e.,

$$C_Y\left(i\right)={\mathit{1}}_N⊛{\mathit{1}}_N⊛C_x\left(i\right)=\sum _{j=-N+1}^{N-1}(N-|j\left|\right)C_x(i-j).$$

(21)

Note that if the input process is stationary, then the input and output processes are jointly stationary.

4.2. Multiple Random Processes

The major complication with observing, evaluating, and processing multiple random processes is how to achieve their time alignment and amplitude normalization (scaling). Focusing here on the time alignment problem only, denote the discrete time observation instances of L random processes as,

$${\mathit{t}}_l=\{t_{l1}<t_{l2}<\dots <t_{l{N}_l}\},l=1,2,\dots ,L.$$

Assume that the first time instance $t_{l1}$ of every process serves as a reference. Then, there are $(L-1)$ uncertainties in time alignment of L processes, i.e.,

$${\Delta }_l=(t_{l1}-t_{11})\in \mathcal{R},l=2,3,\dots ,L.$$

The $(L-1)$ values ${\Delta }_l$ are unknown parameters which must be estimated. Note also that the difference,

$${\Delta }_l-{\Delta }_k=t_{l1}-t_{k1}$$

represents an unknown time shift between the process $x_l\left(t\right)$ and $x_k\left(t\right)$.

For any multivariate stationary distribution of observations of two random processes, the corresponding cross-correlation normally attains a maximum for some time shift between these processes [12]. Hence, a usual strategy for aligning the observed sequences is to locate the maximum value of their cross-correlation. The time shifts ${\Delta }_l$, $l=1,2,\dots ,L$, are then estimated as,

$${\widehat{\Delta }}_l=\underset{\Delta }{argmax}{C}_{x_1{x}_l}(\Delta ),\Delta \in {\left\{(t_{li}-t_{11})\right\}}_{i=1,\dots ,N_l}.$$

Assuming the center values ${\overline{X}}_1$ and ${\overline{X}}_2$ in (18) as scalar representations of the vectors ${\mathit{X}}_1$ and ${\mathit{X}}_2$, their cross-covariance can be computed as,

$$cov\left[|{\mathit{X}}_1{|}_1,{\left|{\mathit{X}}_2\right|}_1\right]=N^{2}cov\left[{\overline{X}}_1,{\overline{X}}_2\right]=N^{2}E\left[({\overline{X}}_1-E\left[X_1\right])({\overline{X}}_2-E\left[X_2\right])\right].$$

(22)

The task now is how to generalize the pairwise cross-covariance (22) to the case of multiple random vectors having possibly different lengths. If all random vectors of interest are concatenated into one single vector, the m-th joint central sum-moment can be defined by utilizing Claim 1.

Definition 5.

The m-th central sum-moment for L random processes with $N_l$ observations, $l=1,2,\dots ,L$, is computed as,

$$\sum {\mu }_m({\mathit{X}}_1,\dots ,{\mathit{X}}_L)={\mu }_m\left(\right|{\mathit{X}}_1{|}_1+\dots +|{\mathit{X}}_L{|}_1)=E\left[{\left|\sum _{l=1}^L\sum _{i=1}^{N_l}(X_{li}-{\overline{X}}_{li})\right|}^m\right].$$

Lemma 12.

The second central sum-moment for L random processes with $N_l$ observations, $l=1,2,\dots ,L$, is equal to the sum of all pairwise covariances, i.e.,

$$\sum {\mu }_2({\mathit{X}}_1,\dots ,{\mathit{X}}_L)=\sum _{l,k=1}^L\sum _{i,j=1}^{N_l}cov\left[X_{li},X_{kj}\right]=\sum _{l,k=1}^L{|cov\left[{\mathit{X}}_l,{\mathit{X}}_k\right]|}_1=var\left[|{\mathit{X}}_1{|}_1+\dots +{\left|{\mathit{X}}_L\right|}_1\right].$$

Proof. 

The expression can be obtained by expanding the sum and then applying the expectation. □

Many other properties of central and noncentral sum-moments can be obtained. For example, assuming two equal-length vectors ${\mathit{X}}_1$ and ${\mathit{X}}_2$, it is straightforward to show that,

$$\begin{array}{ccc}\hfill E\left[{\left(|{\mathit{X}}_1{|}_1+{\left|{\mathit{X}}_2\right|}_1\right)}^2\right]-E\left[|{\mathit{X}}_1{|}_1^2+{\left|{\mathit{X}}_2\right|}_1^2\right]& =& 2\sum _{i,j=1}^{N}E\left[X_{1i}{X}_{2j}\right]\hfill \\ \hfill E\left[{\left(|{\mathit{X}}_1{|}_1+{\left|{\mathit{X}}_2\right|}_1\right)}^2\right]-E\left[{∥{\mathit{X}}_1-{\mathit{X}}_2∥}_2^2\right]& =& \sum _{\begin{array}{c}i,j=1\\ i\ne j\end{array}}^{N}E\left[X_{1i}{X}_{1j}+X_{2i}{X}_{2j}\right]+2\sum _{i,j=1}^{N}E\left[X_{1i}{X}_{2j}\right]+2\sum _{i=1}^{N}E\left[X_{1i}{X}_{2i}\right].\hfill \end{array}$$

5. Illustrative Examples

This section provides examples to quantitatively evaluate the results developed in the previous sections. In particular, the accuracy of approximate linear LS regression proposed in Section 3.1 is assessed to justify its lower computational complexity. The central sum-moments introduced in Section 4 are compared assuming correlated Gaussian processes. Finally, several signal processing problems involving the 1st order Markov process are investigated.

5.1. Linear Regression

Consider a classical one-dimensional linear LS regression problem with independent and identically normally distributed errors. The errors are also assumed to be independent from all other data model parameters. The data points are generated as,

$$X_i=\Delta i:Y_i=P_2{X}_i+P_1+E_i,i=1,2,\dots ,N$$

where $E_i$ are zero-mean, uncorrelated Gaussian samples having the equal variance ${\sigma }_e^2$, and $P_1$ and $P_2$ are unknown parameters to be estimated. This LS problem can be solved exactly using the expression (2), and substituting $w_{1i}=1$ and $w_{2i}=\Delta i$, $\forall i=1,2,\dots ,N$. Alternatively, to avoid inverting the $N\times N$ data matrix, the procedure devised in Section 3.1 suggests to split the data into two equal-size subsets, compute the average data point for each subset, and then solve the corresponding set of two equations with two unknowns. Specifically, the set of two equations with the unknown parameters ${\widehat{P}}_1$ and ${\widehat{P}}_2$ is,

$$\begin{array}{ccc}\frac{2}{N}\sum _{i=1}^{N/2}{Y}_i& =& {\widehat{P}}_1+{\widehat{P}}_2\frac{\Delta }{8}N(2+N)\\ \frac{2}{N}\sum _{i=N/2+1}^N{Y}_i& =& {\widehat{P}}_1+{\widehat{P}}_2\frac{\Delta }{8}N(2+3N)\end{array}$$

(23)

assuming N is even, and using, ${\sum }_{i=1}^{N/2}\Delta i=\frac{\Delta }{8}N(2+N)$, and ${\sum }_{i=N/2+1}^N\Delta i=\frac{\Delta }{8}N(2+3N)$. Denoting ${\overline{Y}}_1=\frac{2}{N}{\sum }_{i=1}^{N/2}{Y}_i$ and ${\overline{Y}}_2=\frac{2}{N}{\sum }_{i=N/2+1}^N{Y}_i$, the closed-form solution of (23) is,

$$\begin{array}{ccc}{\widehat{P}}_1& =& \frac{1}{2N}\left((2+3N){\overline{Y}}_1-(2+N){\overline{Y}}_2\right)\hfill \\ {\widehat{P}}_2& =& \frac{4}{\Delta N^2}({\overline{Y}}_2-{\overline{Y}}_1).\hfill \end{array}$$

As a numerical example, assume the true values $P_1=1.5$, $P_2=0.3$, $E\left[E_i\right]=0$, $var\left[E_i\right]=1$, and $N=40$ data points. Figure 2 shows the intervals $(\overline{T}-\sqrt{var\left[T\right]},\overline{T}+\sqrt{var\left[T\right]})$ versus the subset size $1\le N_1\le N/2$ for the random variable T defined as,

$$T=100\frac{S_{\mathrm{apr}}\left(N\right)-S_{\mathrm{opt}}\left(N\right)}{S_{\mathrm{opt}}\left(N\right)}$$

where $S_{\mathrm{apr}}\left(N\right)={\sum }_{i=1}^N{(Y_i-{\widehat{P}}_1-{\widehat{P}}_2{X}_i)}^2$ and $S_{\mathrm{opt}}\left(N\right)={\sum }_{i=1}^N{(Y_i-P_1-P_2{X}_i)}^2$ are the total MSEs. In the limit, ${lim}_{N\to \infty }(S_{\mathrm{apr}}\left(N\right)-S_{\mathrm{opt}}\left(N\right))=0$, since a sufficiently large subset of data is as good as the complete set of data. For finite N, it is likely that $S_{\mathrm{apr}}\left(N\right)>S_{\mathrm{opt}}\left(N\right)$, so the lower bounds in Figure 2 converge much faster to zero than the upper bounds.

5.2. Comparison of Central Moments

Assuming Lemma 11 and Equation (3), the second central sum-moment of the 2nd order Markov process of length N is,

$$\sum {\mu }_2\left(\mathit{X}\right)=\sum _{i,j=1}^N{C}_{2\mathrm{MP}}(i-j)={\sigma }^2\left(N+2\sum _{i=1}^{N-1}i(1+(N-i)\alpha )e^{-\alpha (N-i)}\right).$$

These moments are compared in Figure 3 for different values of the sequence length N, three values of the parameter $\alpha$, and assuming ${\sigma }^2=1$. It can be observed that the values of the second central sum-moment are increasing with N and the level of correlation, $e^{-\alpha }$.

Consider now the following three central moments of order $m=1,2,\dots$, i.e.,

$$\begin{array}{ccc}\hfill {\overline{S}}_{\mathrm{mnk}}\left(N\right)& =& \sum _{i=1}^{N}E\left[|\sqrt{N}{X}_i{|}^m\right]\hfill \\ \hfill \sum {\mu }_2\left(N\right)& =& E\left[{\left|\sum _{i=1}^N{X}_i\right|}^m\right]\hfill \\ \hfill \sum {\tilde{\mu }}_2\left(N\right)& =& E\left[{\left(\sum _{i=1}^N\left|X_i\right|\right)}^m\right].\hfill \end{array}$$

The moment, ${\overline{S}}_{\mathrm{mnk}}$, is the mean Minkowski distance; the scaling by $\sqrt{N}$ is introduced to facilitate the comparison with the other two moments, i.e., the mean sum-moment $\sum {\mu }_2$, and the mean sum-moment $\sum {\tilde{\mu }}_2$ having the samples summed as the $l_1$ norm. More importantly, assuming a correlated Gaussian process $X_i=x\left(t_i\right)$, the central sum-moment can be readily obtained in a closed form whereas obtaining the closed form expressions for the other two metrics may be mathematically intractable. In particular, by factoring the covariance matrix as, ${\mathit{C}}_x={\mathit{T}}_x{\mathit{T}}_x^T$, the correlated Gaussian vector can be expressed as, $\mathit{X}=\mathit{T}\mathit{U}$. Then the sum of elements, ${\left|\mathit{X}\right|}_1={\mathit{1}}^T\mathit{T}\mathit{U}$, and the m-th central sum-moment can be computed as,

$$\sum {\mu }_2\left(N\right)=E\left[{\left|{\mathit{1}}^T\mathit{T}\mathit{U}\right|}^m\right]={∥{\mathit{1}}^T\mathit{T}∥}_2^{m}E\left[{\left|U\right|}^m\right]={∥{\mathit{1}}^T\mathit{T}∥}_2^m\frac{2^{m/2}}{\sqrt{\pi }}\mathsf{\Gamma }\left(\frac{m+1}{2}\right)$$

where U is a zero-mean, unit-variance Gaussian random variable, and $\mathsf{\Gamma }$ denotes the gamma function.

Figure 4 shows all three moments as a function of sequence length N for three values of parameter $\alpha$ assuming the 1st order Markov process. The vertical axis in Figure 4 is scaled by $1/N$ for convenience. Note that, for uncorrelated, i.e., independent Gaussian samples, the moments $\sum {\mu }_2$ and $\sum {\tilde{\mu }}_2$ are identical. More importantly, all three moments are strictly increasing with the number of samples N and with the moment order m.

2

Finally, the second central sum-moment can be visualized in two dimensions. Consider $N=3$ observations of a zero-mean real-valued stationary random process, $X_i=x\left(t_i\right)$, $i=1,2,3$. Let $\sum {\mu }_2=E\left[{(X_1+X_2+X_3)}^2\right]={\sum }_{i,j=1}^{3}E\left[X_i{X}_j\right]$. Assuming the 1st order Markov correlation model, $E\left[X_i{X}_j\right]\propto e^{-0.5\left|\tau \right|}$, Figure 5 shows the values of the central sum-moment $\sum {\mu }_2$ versus the sample distances ${\tau }_1=(t_2-t_1)$ and ${\tau }_2=(t_3-t_1)$. Several symmetries can be observed in Figure 5. In particular, the central sum-moment $\sum {\mu }_2$ is symmetric about the axis ${\tau }_1={\tau }_2$ as well as about the axis ${\tau }_1=-{\tau }_2$. These symmetries are consequences of the following equalities:

$$\left.\begin{array}{ccc}\sum {\mu }_2({\tau }_1,{\tau }_2)& =& \sum {\mu }_2({\tau }_2,{\tau }_1)\\ \sum {\mu }_2({\tau }_1,{\tau }_2)& =& \sum {\mu }_2(-{\tau }_1,-{\tau }_2)\end{array}\right\}\forall {\tau }_1,{\tau }_2.$$

5.3. Signal Processing Problems for the 1st Order Markov Process

Consider the 1st order Markov process observed at the output of a length N MA filter. According to (21), the (auto-) covariance of the output process is,

$$C_{1\mathrm{MP}+\mathrm{MA}}(k;\alpha )=\sum _{j=-N+1}^{N-1}(N-|j\left|\right){\sigma }^2{e}^{-\alpha |k-j|}.$$

(24)

Conjecture 1.

The MA filtering of the 1st order Markov process generates nearly a 2nd order Markov process.

The parameter of the 2nd order Markov process approximating the combined (auto-) covariance (24) can be obtained using the LS regression fit, i.e.,

$$\widehat{\alpha }=\underset{\tilde{\alpha }}{argmin}\sum _k{(C_{1\mathrm{MP}+\mathrm{MA}}(k;\alpha )-C_{2\mathrm{MP}}(k;\tilde{\alpha }))}^2.$$

(25)

Substituting (3) and (24) into (25), and letting the first derivative to be equal to zero, the LS estimate $\widehat{\alpha }$ must satisfy the linear equation,

$$\sum _k\widehat{\alpha }\left|k\right|+1+W_{-1}\left(-\frac{C_{1\mathrm{MP}+\mathrm{MA}}(k;\alpha )}{e}\right)=0$$

which can be readily solved for $\widehat{\alpha }$, and $W_{-1}$ denotes the Lambert function [24].

A discrete time sequence of N elements has the (auto-) covariance constrained to $(2N-1)$ time indexes as indicated in (24). Assuming the length N MA filter, and that there are $\left(n_{\mathrm{x}}N\right)$ samples, $n_{\mathrm{x}}=1,2,\dots$, of the 1st order Markov process available, the (auto-) covariance (24) has the overall length $2N(n_{\mathrm{x}}+1)-3$ samples. Figure 6 compares the MSE,

$$\mathrm{MSE}\left(n_{\mathrm{x}}\right)=100\times \frac{{\sum }_k{(C_{1\mathrm{MP}+\mathrm{MA}}(k;\alpha )-C_{2\mathrm{MP}}(k;\widehat{\alpha }))}^2}{{\sum }_k{C}_{1\mathrm{MP}+\mathrm{MA}}^2(k;\alpha )}$$

of the LS fit of the (auto-) covariance of the 2nd order Markov process to the combined (auto-) covariance of the 1st order Markov process and the impulse response of the MA filter assuming three values of $\alpha$ and two values of $n_{\mathrm{x}}$. Given $\alpha$ and $n_{\mathrm{x}}$, Figure 6 shows that the best LS fit occurs for a certain value of the MA filter length N. It can be concluded that, in general, the 1st order Markov process changes to the 2nd order Markov process at the output of the MA filter.

2

The second problem to investigate is a linear MMSE (LMMSE) prediction of the 1st order Markov process observed at the output of a MA filter. In particular, given N samples $X_i=x\left(t_i\right)$, $i=1,2,\dots ,N$, of random process having the (auto-) covariance (24), the task is to predict its future value, $X_{N+1}=x\left(t_{N+1}\right)$, $t_{N+1}>t_N$.

In general, the impulse response $\mathit{h}$ of the LMMSE filter to estimate an unknown scalar parameter P from the measurements $\mathit{X}$ is computed as [15],

$$\mathit{h}=E\left[(\mathit{x}-\overline{\mathit{x}}){(\mathit{X}-\overline{\mathit{X}})}^T\right]E\left[P\mathit{X}\right].$$

(26)

Here, the unknown parameter $P=X_{N+1}$, and $E\left[X_{N+1}{X}_i\right]=C_{1\mathrm{MP}+\mathrm{MA}}(N+1-i)$ and $E\left[X_i{X}_j\right]=C_{1\mathrm{MP}+\mathrm{MA}}(i-j)$ in (26) which gives the length N LMMSE filter,

$$\mathit{h}=\left[\underset{N-1}{\underbrace{00\dots 0}}{C}_{1\mathrm{MP}+\mathrm{MA}}\left(1\right)\right].$$

Consequently, the predicted value, $X_{N+1}=X_N{C}_{1\mathrm{MP}+\mathrm{MA}}\left(1\right)/{\sigma }_X^2$. Note that the same procedure, but excluding the MA filter, gives the LMMSE estimate, $X_{N+1}=X_N{C}_{1\mathrm{MP}}\left(1\right)/{\sigma }_X^2$.

The last problem to consider is a time alignment of two zero-mean, jointly stationary processes. It is assumed that the normalized cross-covariance of these two processes is,

$$\frac{E\left[X_{1i}{X}_{2j}\right]}{\sqrt{E\left[X_{1i}^2\right]E\left[X_{2i}^2\right]}}=e^{-\alpha |i-j|}(1+\alpha |i-j\left|\right).$$

(27)

Denote the uncertainty in determining the difference, $(i-j)$, as $\Delta$. In order to estimate the unknown parameters $\alpha$ and $\Delta$, the left-hand side of (27) can be estimated by the method of moments, i.e., let $E\left[X_i^2\right]\approx \frac{1}{N}{\sum }_{i=1}^N{X}_i^2$, and, $E\left[X_{1i}{X}_{2j}\right]\approx \frac{1}{N}{\sum }_{i=1}^N{X}_{1i}{X}_{2|i-\Delta |}$. The cross-covariance (27) can be then rewritten as,

$$v_k=\frac{{\sum }_{i=1}^N{X}_{1i}{X}_{2|i-\Delta -k|}}{\sqrt{{\sum }_{i=1}^N{X}_{1i}^2}\sqrt{{\sum }_{i=1}^N{X}_{2i}^2}}=e^{-\alpha |\Delta +k|}(1+\alpha |\Delta +k\left|\right),k=0,1,2,\dots$$

Utilizing the Lambert function $W_{-1}$, the cross-covariance can be rewritten further as,

$$\alpha |\Delta +k|=-1-W_{-1}\left(-\frac{v_k}{e}\right)\triangleq {\tilde{v}}_k,k=0,1,2,\dots$$

(28)

Assuming, without loss of generality, that $\Delta \ge 0$, the absolute value in (28) can be ignored. Consequently, the unknown parameters $\alpha$ and $\Delta$ can be obtained as a linear LS fit to N measured and calculated values ${\tilde{v}}_k$ in the linear model (28).

6. Conclusions

The development of a novel statistical measure to enable correlation analysis for multiple random vectors resumed by summarizing background knowledge on statistical description of discrete time random processes. This was then extended with the derivation of several supporting results which were used in the following sections. Specifically, it was shown that linear regression can be effectively approximated by splitting the data into disjoint subsets and assuming only one average data point within each subset. In addition, a procedure for generating multiple Gaussian processes with the prescribed autocovariance and cross-covariance was devised. The main result of the paper was obtained by assuming the Taylor’s expansion of multivariate symmetric scalar functions, and then approximating the Taylor’s expansion by a univariate polynomial. The single polynomial variable is a simple sum of variables in the original multivariate function. The polynomial approximation represents a mapping from multiple discrete time observations of a random process to a multidimensional scalar field. The mean field value is a weighted sum of canonical central moments with increasing orders. These moments were named central sum-moments to reflect how they are defined. The sum-moments were then discussed in light of other similar concepts such as total variance, Mean Field Theory, and moving average sequence filtering. Illustrative examples were studied in the last section of the paper. In particular, the accuracy of approximate linear regression was evaluated quantitatively assuming two disjoint data subsets. Assuming the 1st and the 2nd order Markov processes, the central sum-moments were compared with the mean Minkowski distance. For Gaussian processes, the central sum-moments can be obtained in closed form. The remaining problems investigated moving average filtering of the 1st order Markov processes and its prediction using a linear MMSE filter.