Steven T. Smith

The Geometry of Algorithms with Orthogonality Constraints

In this paper we develop new Newton and conjugate gradient algorithms on the Grassmann and Stiefel manifolds. These manifolds represent the constraints that arise in such areas as the symmetric eigenvalue problem, nonlinear eigenvalue problems, electronic structures computations, and signal processing. In addition to the new algorithms, we show how the geometrical framework gives penetrating new insights allowing us to create, understand, and compare algorithms. The theory proposed here provides a taxonomy for numerical linear algebra algorithms that provide a top level mathematical view of previously unrelated algorithms. It is our hope that developers of new algorithms and perturbation theories will benefit from the theory, methods, and examples in this paper.

Covariance, subspace, and intrinsic Crame/spl acute/r-Rao bounds

Crame/spl acute/r-Rao bounds on estimation accuracy are established for estimation problems on arbitrary manifolds in which no set of intrinsic coordinates exists. The frequently encountered examples of estimating either an unknown subspace or a covariance matrix are examined in detail. The set of subspaces, called the Grassmann manifold, and the set of covariance (positive-definite Hermitian) matrices have no fixed coordinate system associated with them and do not possess a vector space structure, both of which are required for deriving classical Crame/spl acute/r-Rao bounds. Intrinsic versions of the Crame/spl acute/r-Rao bound on manifolds utilizing an arbitrary affine connection with arbitrary geodesics are derived for both biased and unbiased estimators. In the example of covariance matrix estimation, closed-form expressions for both the intrinsic and flat bounds are derived and compared with the root-mean-square error (RMSE) of the sample covariance matrix (SCM) estimator for varying sample support K. The accuracy bound on unbiased covariance matrix estimators is shown to be about (10/log 10)n/K/sup 1/2/ dB, where n is the matrix order. Remarkably, it is shown that from an intrinsic perspective, the SCM is a biased and inefficient estimator and that the bias term reveals the dependency of estimation accuracy on sample support observed in theory and practice. The RMSE of the standard method of estimating subspaces using the singular value decomposition (SVD) is compared with the intrinsic subspace Crame/spl acute/r-Rao bound derived in closed form by varying both the signal-to-noise ratio (SNR) of the unknown p-dimensional subspace and the sample support. In the simplest case, the Crame/spl acute/r-Rao bound on subspace estimation accuracy is shown to be about (p(n-p))/sup 1/2-1/2/SNR/sup -1/2/ rad for p-dimensional subspaces. It is seen that the SVD-based method yields accuracies very close to the Crame/spl acute/r-Rao bound, establishing that the principal invariant subspace of a random sample provides an excellent estimator of an unknown subspace. The analysis approach developed is directly applicable to many other estimation problems on manifolds encountered in signal processing and elsewhere, such as estimating rotation matrices in computer vision and estimating subspace basis vectors in blind source separation.

Statistical resolution limits and the complexified Crame/spl acute/r-Rao bound

Array resolution limits and accuracy bounds on the multitude of signal parameters (e.g., azimuth, elevation, Doppler, range, cross-range, depth, frequency, chirp, polarization, amplitude, phase, etc.) estimated by array processing algorithms are essential tools in the evaluation of system performance. The case in which the complex amplitudes of the signals are unknown is of particular practical interest. A computationally efficient formulation of these bounds (from the perspective of derivations and analysis) is presented for the case of deterministic and unknown signal amplitudes. A new derivation is given using the unknown complex signal parameters and their complex conjugates. The new formula is readily applicable to obtaining either symbolic or numerical solutions to estimation bounds for a very wide class of problems encountered in adaptive sensor array processing. This formula is shown to yield several of the standard Crame/spl acute/r-Rao results for array processing, along with new results of fundamental interest. Specifically, a new closed-form expression for the statistical resolution limit of an aperture for any asymptotically unbiased superresolution algorithm (e.g., MUSIC, ESPRIT) is provided. The statistical resolution limit is defined as the source separation that equals its own Crame/spl acute/r-Rao bound, providing an algorithm-independent bound on the resolution of any high-resolution method. It is shown that the statistical resolution limit of an array or coherent integration window is about 1.2/spl middot/SNR/sup -1/4/ relative to the Fourier resolution limit of 2/spl pi//N radians (large number N of array elements). That is, the highest achievable resolution is proportional to the reciprocal of the fourth root of the signal-to-noise ratio (SNR), in contrast to the square-root (SNR/sup -1/2/) dependence of standard accuracy bounds. These theoretical results are consistent with previously published bounds for specific superresolution algorithms derived by other methods. It is also shown that the potential resolution improvement obtained by separating two collinear arrays (synthetic ultra-wideband), each with a fixed aperture B wavelengths by M wavelengths (assumed large), is approximately (M/B)/sup 1/2/, in contrast to the resolution improvement of M/B for a full aperture. Exact closed-form results for these problems with their asymptotic approximations are presented.

Optimum phase-only adaptive nulling

Optimum weighting of adaptive antenna arrays is accomplished by computing the weight vector that maximizes the signal-to-interference-plus-noise ratio (SINR). The optimal weight vector is in general complex with each weight having different magnitudes and phases. However, the design of some antenna arrays facilitates phase-only weighting or phase-only adjustments, whereupon it is desirable to compute the constrained optimal weight vector whose components have fixed magnitudes but variable phases. This constrained optimization problem may be posed as the problem of maximizing the SINR on the space of phase-only vectors. This paper addresses the problem of computing optimal phase-only adaptive weight vectors by exploiting several properties of phasor and matrix algebra. Two new algorithms (the phase-only conjugate gradient and phase-only Newtons method) are introduced. The convergence properties. SINR performance, sidelobe level performance, and nulling performance of these algorithms are demonstrated using simulations and experimental data.

On conjugate gradient-like methods for eigen-like problems

Numerical analysts, physicists, and signal processing engineers have proposed algorithms that might be called conjugate gradient for problems associated with the computation of eigenvalues. There are many variations, mostly one eigenvalue at a time though sometimes block algorithms are proposed. Is there a correct “conjugate gradient” algorithm for the eigenvalue problem? How are the algorithms related amongst themselves and with other related algorithms such as Lanczos, the Newton method, and the Rayleigh quotient?

Structured covariance estimation for space-time adaptive processing

Adaptive algorithms require a good estimate of the interference covariance matrix. In situations with limited sample support such an estimate is not available unless there is structure to be exploited. In applications such as radar space-time adaptive processing (STAP) the underlying covariance matrix is structured (e.g., block Toeplitz), and it is possible to exploit this structure to arrive at improved covariance estimates. Several structured covariance estimators have been proposed for this purpose. The efficacy of several of these are analyzed in this paper in the context of a variety of STAP algorithms. The SINR losses resulting from the different methods are compared. An example illustrating the superior performance resulting from a new maximum likelihood algorithm (based upon the expectation-maximization algorithm) is demonstrated using simulation and experimental data.

Bayesian Discovery of Threat Networks

A novel unified Bayesian framework for network detection is developed, under which a detection algorithm is derived based on random walks on graphs. The algorithm detects threat networks using partial observations of their activity, and is proved to be optimum in the Neyman-Pearson sense. The algorithm is defined by a graph, at least one observation, and a diffusion model for threat. A link to well-known spectral detection methods is provided, and the equivalence of the random walk and harmonic solutions to the Bayesian formulation is proven. A general diffusion model is introduced that utilizes spatio-temporal relationships between vertices, and is used for a specific space-time formulation that leads to significant performance improvements on coordinated covert networks. This performance is demonstrated using a new hybrid mixed-membership blockmodel introduced to simulate random covert networks with realistic properties.Network detection is an important capability in many areas of applied research in which data can be represented as a graph of entities and relationships. Oftentimes the object of interest is a relatively small subgraph in an enormous, potentially uninteresting background. This aspect characterizes network detection as a “big data” problem. Graph partitioning and network discovery have been major research areas over the last ten years, driven by interest in internet search, cyber security, social networks, and criminal or terrorist activities. The specific problem of network discovery is addressed as a special case of graph partitioning in which membership in a small subgraph of interest must be determined. Algebraic graph theory is used as the basis to analyze and compare different network detection methods. A new Bayesian network detection framework is introduced that partitions the graph based on prior information and direct observations. The new approach, called space-time threat propagation, is proved to maximize the probability of detection and is therefore optimum in the Neyman-Pearson sense. This optimality criterion is compared to spectral community detection approaches which divide the global graph into subsets or communities with optimal connectivity properties. We also explore a new generative stochastic model for covert networks and analyze using receiver operating characteristics the detection performance of both classes of optimal detection techniques.

Space-time clutter covariance matrix computation and interference subspace tracking

A fast method of computing ideal space-time clutter covariance matrices is described and the methods application to the study of clutter interference subspaces is considered. These areas are of principal importance in the study of space-time adaptive processing (STAP) for airborne pulse Doppler radar arrays. Formulae for clutter covariance matrices assuming various stochastic models are first derived, then a fast algorithm is developed in the case of uncorrelated patch-to-patch clutter with arbitrary illumination and reflectance. Finally, the clutter interference subspace as a function of scan angle for the case of a rotating sensor array is studied.

Subspace tracking with full rank updates

A new method of trading the principal invariant subspaces of a time-varying covariance matrix is proposed. The method addresses the case encountered frequently in applications where the covariance is updated by a full rank matrix at each time step; it is not assumed that the covariance changes by rank-one updates. In contrast to subspace tracking algorithms that exploit the algebraic structure that rank-one updates provide, the proposed algorithm uses the geometric structure of full rank updates. The main idea is to determine the time derivative of the subspace using the subspace tracking equation (introduced here), perform an optimization line search in this direction, then finish with a (truncated) Newton, method comparable to Rayleigh quotient iteration. The algorithm is performed on the constraint surface (Grassmann manifold) of matrices with orthonormal columns. The convergence rate and cost of this method are compared with other subspace tracking and standard eigenvalue decomposition (EVD) algorithms using the example problem of tracking the time-varying clutter interference of a rotating sensor array. Simulations indicate that the proposed method is cheaper than a full EVD, but its superlinear convergence rate and higher overhead make it about 25% more costly than a linearly convergent subspace iteration method. Nevertheless, the proposed methods generality makes it appropriate for nonlinear eigenvalue problems and other time-varying problems where linear EVD algorithms cannot be applied.

Linear and nonlinear conjugate gradient methods for adaptive processing

Both fully adaptive and partially adaptive processing algorithms can be implemented using linear and nonlinear conjugate gradient (CG) methods. This approach can be computationally attractive when a family of continuous time-varying adaptive problems must be solved. An iterative approach that starts with the previous solution can converge quickly to the new solution. This paper considers linear and nonlinear CG methods for adaptive processing. The example used is space-time adaptive processing (STAP) and clutter interference subspace tracking within the context of an airborne radar with a rotating array. It is seen that the linear CG method for adaptive filtering in this context can be several times less costly than sample matrix inversion if the interference subspace has sufficient eigenvalue structure, and that a nonlinear CG algorithm is capable of tracking the principal clutter interference subspace.