I. INTRODUCTION
Dimension reduction techniques are often applied as a data pre-processing step or as a part of the data analysis to simplify a data model. This typically involves identification of a suitable low-dimensional representation for the original high-dimensional data set or a subspace of an observation space. By working with this reduced representation, tasks such as classification or estimation can often yield more accurate and readily interpretable results, while computational costs may be significantly reduced.
In mathematical terms, the problem we consider can be stated as follows: given the N-dimensional random variable x = (x 1, x 2, … , x N)T, find a lower dimensional representation of it, y = (y 1, y 2, … , y m)T with m < N, that captures the content in the original data, according to some criterion.
Among the linear dimension reduction techniques, principal component analysis (PCA) is the best in the mean-square error sense [Reference Jolliffe1]. In various fields, it is also known as the singular value decomposition (SVD) and the Karhunen–Loeve transform. In essence, PCA seeks to reduce the dimension of the data by finding a few orthogonal linear combinations (the PCs) of the original variables. The first PC y 1 is the linear combination with the largest variance. The second PC y 2 is the linear combination with the second largest variance and orthogonal to the first PC, and so on. There are as many PCs as the number of original variables. For many data sets, the first several PCs explain most of the variance, so that the rest can be disregarded with minimal loss of information.
Here we would like to use dimension reduction techniques in a signal parameter estimation problem. We hope to obtain a similar PCA technique with a different criterion.
II. DATA DIMENSION REDUCTION
The estimation problem, which we consider here, can be phrased as follows: given x(t) = s(t, θ)+w(t), where s(t, θ), 0 ≤ t ≤ T is a signal of the known waveform to be observed in the presence of an additive Gaussian noise process w(t) with zero mean. We have an a priori probability density p(θ) in the random parameter estimation problem. In all of our discussions we assume that p(θ) is known. We need an estimation parameter θ of a signal. In general, the variable θ appears in a signal in a nonlinear manner.
We assume that the observation space corresponds to the set of N observations: x 1, x 2, … , x N. Thus, each set can be thought of as a point in an N-dimensional space and can be denoted by a column vector x = s(θ)+w, where s(θ) ∈ R N and w ∈ R N are the N-dimensional vectors of a signal and a noise correspondingly. Vector w has a nonsingular covariance matrix Rw. Then the probability density of x is
![\eqalignno{p\lpar {\bf x}\vert \theta\rpar &= \lpar \lpar 2\pi\rpar ^N \vert {\bf R}_w\vert\rpar ^{1/2} \exp \cr &\quad \times \left[-{1 \over 2}\lpar {\bf x} - {\bf s}\lpar \theta\rpar \rpar ^T {\bf R}_w^{-1}\lpar {\bf x} - {\bf s}\lpar \theta\rpar \rpar \right].&}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqn1.gif?pub-status=live)
To obtain the m-dimensional vector y with m < N, we use linear transformation y = Cx with the transformation matrix C . We need such a matrix C that guarantees minimal loss of estimation accuracy of a parameter θ, using vector y. In addition to the foregone requirement we try to represent x in a new coordinate system in which the components are statistically independent random variables: CRCT = I, where I is a diagonal identity matrix. It is convenient to write the transformation matrix as C = AR−1/2w. Here R−1/2w is a symmetric square root from R−1w(R−1/2wR−1/2w = R−1w),
AT = (a1, … , am), AAT = I. So we have
First of all, recall that the variance of any unbiased estimate of an arbitrary parameter θ is determined from the Cramer–Rao inequality [Reference Van Trees2]:
![\sigma_\theta^2\ge\left(\left[{{\partial\ln p\lpar {\bf x}\vert\theta\rpar }\over{\partial\theta}}\right]^2\right)^{-1} = \left(\left[-{{\partial^2\ln p\lpar {\bf x}\vert\theta\rpar }\over{\partial\theta^2}}\right]\right)^{-1}.](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqn3.gif?pub-status=live)
The inequalities are defined if the derivatives involved are existent and absolutely integrable. Recall that the denominator of the right-hand sides of CRI is usually referred to as the Fisher information on parameter θ:
![I\lpar \theta\rpar = {\left[{{\partial\ln p\lpar {\bf x}\vert\theta\rpar }\over{\partial\theta}}\right]}^2 = -{{\partial^2\ln p\lpar {\bf x}\vert\theta\rpar }\over{\partial\theta^2}}.](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqn4.gif?pub-status=live)
where subscript N is used to distinguish the initial dimension of the observation from a new reducing dimension m, and 
is the column vector of derivatives.
The Fisher information in the vector y [Reference Latyshev3] is
![I_m\lpar \theta\rpar = \sum_{k = 1}^m{\left[{\bf a}_k^T{\bf R}_w^{-1/2}{\bf s}^\prime\lpar \theta\rpar \right]}^2.](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqnU1.gif?pub-status=live)
The loss of the Fisher information is
![\eqalign{\Delta I\lpar \theta\rpar &= I_N\lpar \theta\rpar - I_m\lpar \theta\rpar \cr &=\lpar {\bf s}^\prime\lpar \theta\rpar \rpar ^T{\bf R}_w^{-1}{\bf s}^\prime\lpar \theta\rpar - \sum_{k = 1}^m{\left[{\bf a}_k^T{\bf R}_w^{-1/2}{\bf s}^\prime\lpar \theta\rpar \right]}^2.}](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqnU2.gif?pub-status=live)
The mean of the loss of the Fisher information is
![\Delta I = E_\theta\left\{\lpar {\bf s}^\prime\lpar \theta\rpar \rpar ^T{\bf R}_w^{-1}{\bf s}^\prime\lpar \theta\rpar - \sum_{k = 1}^m{\left[{\bf a}_k^T{\bf R}_w^{-1/2}{\bf s}^\prime\lpar \theta\rpar \right]}^2\right\}\comma \;](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqn6.gif?pub-status=live)
where E θ denotes an expectation over the random variable θ. Thus, we need the transformation matrix (Reference Van Trees2) which provides minimal value of ΔI.
Theorem: The linear transformation with the matrix ARw−1/2 provides minimal mean of the loss of Fisher information ΔI, if the column vectors a1, … , am of AT are the orthonormal eigenvectors of
corresponding to m largest eigenvalues. At the same time
where 
are the eigenvalues of B.
Proof: Rewrite (3) in the next form
![\Delta I = E_\theta\left\{\lpar {\bf s}^\prime\lpar \theta\rpar \rpar ^T{\bf R}_w^{-1}{\bf s}^\prime\lpar \theta\rpar \right\}- E_\theta \left\{\sum_{k = 1}^m{\left[{\bf a}_k^T{\bf R}_w^{-1/2}{\bf s}^\prime\lpar \theta\rpar \right]}^2\right\}.](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqnU3.gif?pub-status=live)
The first term does not depend on ak. Therefore, we find the minimal value of ΔI if the subtrahend is maximal. Denote it as H(a1, … , am). Inverting averaging with summation and taking into account the equality
![{\left[{\bf a}_k^T{\bf R}_w^{-1/2}{\bf s}^\prime\lpar \theta\rpar \right]}^2 = {\bf a}_k^T{\bf R}_w^{-1/2}{\bf s}^\prime\lpar \theta\rpar \lpar {\bf s}^\prime\lpar \theta\rpar \rpar ^T{\bf R}_w^{-1/2}{\bf a}_k\comma \;](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqnU4.gif?pub-status=live){kind=small}
we have
![H\lpar {\bf a}_1\comma \;\ldots\comma \; {\bf a}_m\rpar = \sum_{k = 1}^m{\bf a}_k^T\left[{\bf R}_w^{-1/2}E_\theta\left\{{\bf s}^\prime\lpar \theta\rpar \lpar {\bf s}^\prime\lpar \theta\rpar \rpar ^T\right\}{\bf R}_w^{-1/2}\right]{\bf a}_k.](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151022070154772-0355:S175907870900021X_eqnU5.gif?pub-status=live)
The expression in brackets is a symmetric matrix:
In compliance with the theorem about eigenvalues and eigenvectors [Reference Fukunaga4], the maximal value of H (a1, … , am) takes place if a1, … , am are the orthonormal eigenvectors of the matrix B, corresponding to m largest eigenvalues 
and
The equality 
implies
On the other hand, 
. It implies 
.
Note that we can assert that a subspace spanned by the column vectors a1, … , am is the m-dimensional subspace of an observation space with maximal Fisher information content concerning the parameter θ among any other m-dimensional subspaces.
III. NUMERICAL EXAMPLE
To illustrate the advantages of the dimension reduction technique in signal parameter estimation, let us consider the Doppler shift problem. We model the next signal
Here s(t, f D) is a signal of known form to be observed in the presence of an additive white Gaussian noise process w(t) with zero mean and variance σw2. The time delay is known and we wish to estimate the Doppler frequency shift f D only. Note that for the radar case, the standard narrowband assumption is employed here, i.e. the frequency offset f D is a narrowband approximation to the stretching or shrinking of the frequency axis due to the Doppler effect induced by the relative motion of the reflecting target.
The quantity f D is a random parameter with a uniform a priori density p 0(f D), f D ∈[−F;F], so we may use the Bayesian estimation procedure. However, it is very complicated to put it into practice. In this example, our goal is to simplify the Bayesian estimate for the cases when a cost function is the square of the error. For the problem at hand, we have to capture Fisher information concerning Doppler shift after data dimension reduction.
To illustrate our results we use the FM signal with T · Δf = 32:
A set of N observations, x 1, x 2, …, x N, consists of statistically independent values x k with variance σw2. So we have a column vector x = s(f D)+w ∈ R N with the nonsingular covariance matrix Rw = σw2I.
The Bayesian estimate uses the a posteriori density:
The estimate is the mean of the a posteriori density (or the conditional mean):
The simplification of the Bayesian estimate may be achieved using linear transformation with matrix A to provide minimal mean loss of Fisher information about Doppler shift. From the foregoing theorem the column vectors a1, … , am of AT are the orthonormal eigenvectors of the matrix
corresponding to m largest eigenvalues. Here s′ (f D) is the column vector of derivatives with respect to f D.
If we use the vector y = Ax, y ∈ R m, the Bayesian estimate is the same as in (11) and (12), where x is replaced by y.
For the simulations T = 1, N = 512, f D is a random parameter with the uniform a priori density f D ∈ [−0.4; 0.4]. The mean squared error (MSE) of the Doppler shift estimate was calculated for each of the algorithms based on 10 000 Monte Carlo trials for various signal-to-noise ratio (SNR) values. The results are plotted in Fig. 1 together with the appropriate Cramer–Rao bound (CRB). The demonstrated results correspond to the Bayesian estimate using initial 512-dimensional vector x and results based on five-dimensional and 10-dimensional vectors y.
MSE of the Bayesian estimate using initial (512-lower thin line) and reduced dimensional vectors (5, 10).
As can be seen from Fig. 1, the 512-dimensional and the 10-dimensional data sets are equivalent with respect to estimation accuracy. At the same time the computational complexity of the Bayesian algorithm is two orders of magnitude less in this example. This is due to the fact that we capture Fisher information concerning Doppler shift. Further, the effect of dimension reduction from 10 to 5 degrades accuracy noticeably.
IV. SUBSPACE-BASED MAXIMUM-LIKELIHOOD (ML) ESTIMATION
The theorem that is proved above allows us to reveal an m-dimensional subspace of the observation space containing Fisher information on the arbitrary parameter θ. If the informational contents of m-dimensional vector y and N-dimensional x are equal to one another, we may regard y as a vector of sufficient statistics in a sense that an estimate of θ based on y is characterized by the same CRB as the estimate based on the initial observation. Otherwise y may be considered as an approximate vector of sufficient statistics. The degree of approximation is determined by the sum (8).
Next, we consider an approach to the synthesis of subspace-based ML estimation algorithms providing for a higher accuracy of the results. Note that the estimate as obtained according to the ML principle in nonlinear estimation problems is asymptotically optimal for large SNR.
As an illustration, let us consider estimation of the time of arrival. Instead of (9) now we have
with
The unknown parameter is a signal delay τ with respect to the point on the time axis that is chosen as an original. It is necessary to estimate τ from the observation x(t).
We are taking into account that an estimation problem is usually solved either after a signal has been detected or simultaneously with signal detection. Therefore, below we assume that an expected signal has been detected and that its approximate position on the time axis is known. It is necessary to specify the signal's position in the presence of distortions. This assumption fits the real operation conditions of tracking equipment containing a discriminator that in each measurement cycle is tuned to a given reference value of a parameter to be measured.
Let us compare the information content of the vector y with different dimensions for this problem. We can take into account the results of calculation of Fisher information from [Reference Latyshev5, Reference Latyshev6] for signal (10) to be observed in the presence of an additive white Gaussian noise process w(t). Normalized graphs are shown in Figs. 2–4 for different dimensions of vector y and T = 1.
Fisher information in observation (1) and in the output signal of an ML filter (2).
Fisher information in observation (1) and in 2D vector (2).
Fisher information in observation (1) and in 3D vector (2).
It is noteworthy that the first eigenvector a1 coincides with the pulse response of the ML-filter which follows from the ML equation to within a constant [Reference Latyshev5]. Therefore, the solution of the ML equation corresponds to the first coordinate of the vector y. Its informational content is shown in Fig. 2. The graph in Fig. 2 can be interpreted as an informational explanation of the asymptotic optimal property of the ML estimate. The informational content of the observation x and of the one-dimensional subspace are close to each other for small τ only. In an estimation problem the estimate is close to the true value τ in case of high SNR and approximate position on the time axis is known within narrow a priori range (τmin; τmax). In general, if the SNR is not high, the 1D subspace is associated with large losses of Fisher information concerning the time delay. So the first coordinate of the vector y cannot be considered as sufficient statistic for the time delay estimation. Figure 2 shows the losses of Fisher information using the gray shading.
The rest of the coordinates of vector y yield an additional information channel that can be used to enhance the accuracy of estimation. In Figs. 3 and 4 we see informational contents of 2D and 3D subspaces, respectively. The losses are noticeably less. Therefore, we can consider 2D or 3D vectors y as an approximate vector of sufficient statistics within larger a priori range (τmin; τmax).
As long as we have an approximate vector of sufficient statistics, each of its components can be used for an independent estimation of the signal delay. In practice, this estimation procedure necessitates a multi-channel device providing for the formation of the corresponding vector components in its channels. Each of the channels contains a linear filter with an impulse response ak as the corresponding eigenvector of (7).
This approach actually deals with an m-channel device, so that parameter τ is estimated in each channel using the ML estimator (MLE) (see Fig. 5).
A multi-channel estimator of the time of arrival.
Therefore, the results in channels can be integrated in a weight adder (WA). Since an estimation procedure in each channel is nonlinear, it is difficult to exactly calculate its statistical characteristics. However, potential characteristics can be employed. As shown above, all channels are characterized by the mean of the Fisher information λk. Hence, the estimates in the channels can be characterized by the different CRBs for the error variance. The optimal integration in WA involves calculation of the average of the estimates obtained in the channels with the weights proportional to 
.
Simulation of this approach is illustrated in Fig. 6. For comparison, Fig. 6 shows the error variances obtained on the bases of 1D, 2D and 3D approximate vectors of sufficient statistics. Recall that the estimate for the case m = 1 coincides with the classic ML estimate. Each point of the plots has been obtained as a result of 10 000 estimation cycles with various white noise realizations. The initial dimension of the observation is equal to N = 1024. Here it is assumed that the approximate value of the delay is known within an a priori range (−T/4;T/4), T = 1. Thus, the domain of uncertainty is the half-duration of the signal.
MSE of the estimates based on approximate vectors of sufficient statistics.
The simulation results show that the pulse's time of arrival may be estimated more accurately than that in the case when the classical ML approach is applied (m = 1). In the example considered, this effect is the result of averaging independent estimates of the time of arrivals in independent channels. As is known from estimation theory, the averaging effectively enhances the accuracy of the determination of various signal parameters. The method proposed in this paper allows application of the aforementioned effect of averaging.
V. INDEPENDENT TIME ARRIVAL AND DOPPLER SHIFT ESTIMATIONS
Consider the widespread Gaussian observation model. Let the N-dimensional data column vector be an additive mixture of the deterministic component and the distortion vector
where signal s(τ, f D) depends on the two a priori unknown parameters. Assume that both parameters are statistically mutually independent. The parameter vector v = (τ, f D) has a bounded domain of variation, where 2D probability density p 0(τ, f D) is specified. The distortion vector w consists of independent sample values, described by the Gaussian probability distribution with the zero mean and a nonsingular covariance matrix Rw = σw2I.
Assume that it is necessary to estimate τ from the observed vector x. In such a case we can consider f D to be a nuisance parameter. To find an independent estimate of τ, it is desirable to obtain the statistics in the form of linear functions of vector x that are not affected by the nuisance parameter. We consider these statistics as invariant to variations of the nuisance parameter. It is also necessary to minimize a possible deterioration of the estimation accuracy, if τ is estimated using these statistics.
The problem formulated above may be solved on the basis of the following ideas. If the data or the statistics obtained from these data do not contain Fisher information on a certain parameter, it is impossible to estimate this parameter by the use of this statistic. This situation may be interpreted as an invariance to the changes of this parameter. On the other hand, when the observations are used to form a statistic or data set concentrating the complete Fisher information on a desired parameter, this parameter will be estimated with the same accuracy as that provided by all of the observations. Hence, it is necessary to have a tool that ensures control of the Fisher information content in the statistics obtained via linear transformations of the observations. The theorem proved above allows one to reveal the subspaces of an observation space containing Fisher information on specific scalar parameters τ or f D. We use the next two matrices, which are similar to (7):
Here, differentiation is performed with respect to the parameters indicated by the indexes τ or f D. E v denotes an expectation over the vector random variable v = (τ, f D).
In order to eliminate Fisher information on the nuisance parameter f D from an observation, it is possible to obtain a decomposition of the observed signal using the eigenvectors of Bf. Then we remove the terms with nonzero eigenvalues, and sum the remaining terms.
A drawback of this approach is that the Fisher information on desired parameter τ is beyond transformations. We need the eliminated statistics with minimum Fisher information concerning τ. Therefore, retaining the idea of the method applied to find approximate invariant statistics by means of elimination of some terms from the orthogonal decomposition, we modify the approach. Note that the considered transformation of the observed vector x to the vector of decomposition coefficients actually converts Bτ into a diagonal matrix. Denoting the matrix that consists of the eigenvectors of (15) by A, we have ABτAT = D, where D is the diagonal matrix with the eigenvalues of Bτ. It is known that two matrices Bτ and Bf can be diagonalized simultaneously. According to [Reference Fukunaga4], this can be done using the single matrix A so that
where D and AT are chosen to be the matrices of the eigenvalues and the eigenvectors of the matrix Bf−1Bτ, respectively. Simultaneously with performance of the above transformation, it is necessary to separate Fisher information on parameter τ from that on parameter f D.
Following [Reference Fukunaga4], let us introduce the auxiliary matrix
Using matrix A, we can choose a linear transformation so that Q and Bτ are diagonalized simultaneously and the following conditions are fulfilled:
where D(1) is the diagonal matrix with ordered entries 
which are the eigenvalues of matrix Bτ. The substitution (18) into (19) yields
or
Hence, it follows that matrices Bτ and Bf are diagonalized simultaneously under the above conditions. Besides, the eigenvalues of both matrices satisfy the relationship 
. Thus, when Bτ and Bf are diagonalized simultaneously, the greatest eigenvalues of the first matrix correspond to the least eigenvalues of the second matrix and vice versa. An advantage of this approach is that the orthogonal decomposition based on matrix A yields an orthogonal series, whose terms are ranked simultaneously according to the degree of decrease of Fisher information. However, these terms are arranged in reciprocal order: Fisher information on the desired parameter is concentrated in the first decomposition terms and Fisher information on the nuisance parameter is concentrated in the last terms. Naturally, in this situation, elimination of terms with the greatest content of information on the nuisance parameter guarantees both minimal loss of Fisher information on the parameter to be estimated and minimal loss in the accuracy of estimation.
In such a way, on the basis of the generalized eigenvectors of the matrix pair (Bτ, Bf) we can divide the observation space into two mutually orthogonal subspaces containing Fisher information about τ and f D, respectively. For example, to obtain an independent estimate of the time of arrival we have to remove Fisher's information about Doppler shift. For this purpose, we use the subspace with Fisher's information about τ only. Figures 7 and 8 illustrate this possibility for 7 bits M-sequences. Both figures show relief of |l(τ, f D)|, where l(τ, f D) is the logarithm of the likelihood function. Figure 7 corresponds to observation data without noise. We see here the local gap near the true values of τ and f D. Figure 8 corresponds to |l(τ, f D)| using the orthogonal projection of observation data onto the subspace with Fisher information concerning τ only. Here we see that the narrow canyon is definitely parallel to f D axes. It implies that the true value of τ may be obtained independently from f D. If τ is estimated, the next step is to estimate f D using the subspace with Fisher information concerning f D.
Logarithm of the likelihood function for 7 bits M-sequences.
Logarithm of the likelihood function in invariant to f D subspace.
VI. CONCLUSIONS
Based on the Fisher information approach, a linear dimension reduction technique has been proposed. It is similar to PCA, but it makes good use of the minimal loss of Fisher information as a criterion. As a result, we obtain subspaces of an observation space concentrating Fisher information concerning interesting parameters of a signal. The projections of an observation onto these subspaces are the vectors of the sufficient statistics. Simulation performed with the use of the test signals has shown that the accuracy of the estimation is based on the ML method and can be improved if we use these vectors. The generalized eigenvectors of the matrix pair allow obtaining mutually orthogonal subspaces for independent estimation of two different parameters.
Viacheslav Latyshev received the diploma in electrical engineering from the Moscow Aviation Institute in 1971. In 1979, he received the Ph.D. degree in electrical engineering and D.Sc. degree in 1999. From 1971 to 1979, he has held research positions at the Department of Electrical Engineering, Moscow Aviation Institute. In 1985, he was appointed Professor of Signal Processing, where he is working currently. His research interests include stochastic signal processing, fast methods for signal parameter estimation, higher order statistics.
