Despite the widespread occurrence of classification problems and the increasing collection of point process data across many disciplines, study of error probability for point process classification only emerged very recently. Here, we consider classification of renewal processes. We obtain asymptotic expressions for the Bhattacharyya bound on misclassification error probabilities for heavy-tailed renewal processes.

The Hawkes process has been applied successfully to model point process data in a number of application areas including seismology, neural coding, high-frequency finance, and genomics. But the Hawkes process does not accommodate inhibitory effects observed jointly with excitatory behavior, e.g., for regulating neural activity in the brain and binding of transcriptional factors along a genome. In this paper, we extend our previous work on modeling inhibitory effects in a scalar point process to develop a nonlinear Hawkes model for the vector case. We develop a maximum likelihood procedure for the vector point process and demonstrate the algorithm using some neural and genomic data. A comparison with the multivariate linear Hawkes model shows the superiority of the proposed model.

The first order vector autoregression with a symmetric transition matrix (sym-VAR(1)) occurs widely in applications such as electrical networks, robotics and mechanical systems. However, there is so far almost no work on estimating a sym- VAR(1) from time series data that guarantees the estimated system matrix is Schur stable. Here, we develop, for the first time, a closed-form estimator that guarantees symmetry and stability.

Kernel ridge regression is widely used but the theory of its performance has never been fully developed. While there are results on convergence there are few on bias and variance. Here we find expressions for local bias and variance for the important case of the exponential quadratic kernel. Using these new expressions, we explain when quadratic exponential kernel ridge regression can work well and when it will fail.

Subspace tracking has a long and still growing list of applications with recent interest stimulated by the advent of ’streaming’ data. Despite the considerable literature, all of it motivated by the desire to track time varying subspaces, there seems to be almost no analysis of that time varying problem. Rather existing analyses assume the subspace is fixed. Here, for the first time we develop a true subspace tracking analysis in continuous time. It turns out to be a very hard problem and so we tackle the problem of a rank 1 subspace.

While clustering in general is a heavily worked area, clustering of auto-correlated time series (CATS) has received relatively little attention. Here, we develop a convex clustering algorithm suited to auto-correlated time series and compare it with a state of the art method. We find the proposed algorithm is able to more accurately identify the true clusters.

Point processes are finding increasing applications in neuroscience, genomics, and social media. But basic modelling properties are little studied. Here we consider a periodic time-varying Poisson model and develop the asymptotic Cramer-Rao bound. We also develop, for the first time, a maximum likelihood algorithm for parameter estimation.

In a range of network system identification problems point processes are the nodal signals. In many of these applications a fundamental question involves characterizing the internodal dynamics. Here we address that question using mutual information flows and establish a connection with Granger causality. We show how to compute the mutual information flow between the two components of a bivariate Hawkes process. We illustrate with some simulations and data analysis.

Stochastic differential equations (SDEs) evolving in Stiefel manifold have numerous applications in Science and Engineering. While numerical schemes for ordinary differential equations (ODEs) in Stiefel manifolds are reasonably well established, much less has been done for numerical SDEs schemes in Stiefel manifolds. A crucial challenge is to ensure that the trajectory remains on the manifold. But many existing SDE numerical schemes fail to do this. Here we achieve this by extending the so-called 'tangent space parameterization’ (TaSP) for ODEs to SDEs. In so doing we discover a previously missed constraint. We give simulations to illustrate the new numerical scheme.

Attention has turned in recently to clustering of collections networks where each network is characterized by a sparse inverse covariance matrix. Here we consider a mixture model of sparse Gaussian graphical models and develop an exact EM algorithm that improves considerably on a previous approximation. We provide numerical examples to show its superior performance.

Many applications in science and engineering involve tracking the state of a stochastic differential equation (SDE) evolving in a Lie group. This has been tackled by particle filtering although some existing schemes fail to satisfy geometric constraints. Further the conventional particle filter suffers from particle deprivation. Here we overcome these problems by managing the geometry with the Cayley transform and particle depletion with optimal transport. Simulations show the superiority over 'regular' geometry aware particle filtering.

Time-series clustering involves grouping homogeneous time series together based on certain similarity measures. The mixture AR model (MxAR) has already been developed for time series clustering, as has an associated EM algorithm. However, this EM clustering algorithm fails to perform satisfactorily in large-scale applications due to its high computational complexity. This paper proposes a new algorithm, k-ARs, which is a limiting version of the existing EM algorithm. It shows remarkably good computational performance when applied to large-scale clustering problems as illustrated on some benchmark simulations motivated by some real applications.

While there is now a significant literature on sparse inverse covariance estimation, all that literature, with only a couple of exceptions, has dealt only with univariate (or scalar) networks where each node carries a univariate signal. However in many, perhaps most, applications, each node may carry multivariate signals representing multi-attribute data, possibly of different dimensions. Modelling such multivariate (or vector) networks requires fitting block-sparse inverse covariance matrices. Here we achieve maximal block sparsity by maximizing a block-l0-sparse penalized likelihood. There is only one previous algorithm that already does this, but it does not scale. Here we address key computational bottlenecks and develop a new algorithm which is much faster and has massively reduced requirements on matrix conditioning. A benchmark study shows a computational speed-up by many orders of magnitude.

Individual-level resting-state networks (RSNs) based on resting-state fMRI (rs-fMRI) are of great interest due to evidence that network dysfunction may underlie some diseases. Most current rs-fMRI analyses use linear correlation. Since correlation is a bivariate measure of association, it discards most of the information contained in the spatial variation of the thousands of hemodynamic signals within the voxels in a given brain region. Subject-specific functional RSNs using typical rs-fMRI data, are therefore dominated by indirect connections and loss of spatial information and can only deliver reliable connectivity after group averaging. While bivariate partial correlation can rule out indirect connections, it results in connectivity that is too sparse due to lack of sensitivity. We have developed a method that uses all the spatial variation information in a given parcel by employing a multivariate information-theoretic association measure based on canonical correlations. Our method, multivariate conditional mutual information (mvCMI) reliably constructs single-subject connectivity estimates showing mostly direct connections. Averaging across subjects is not needed. The method is applied to Human Connectome Project data and compared to diffusion MRI. The results are far superior to those obtained by correlation and partial correlation.

We sketch some aspects of LMS particularly focussing on some puzzles, problems and potentials. We explain the 'independence heuristic' by means of averaging theory for which we give a simple expose'. We suggest an explanation for the 'urban myth' that white noise only performance formulae can be used as surrogates for the correct performance formulae which involve autocorrelations.We then extend the discussion to more recent network versions such as diffusion LMS. We show that the sharing aspect of network LMS algorithms induces a two time scale structure (something not yet widely known) and exhibit its consequences. We comment on the recent upsurge of interest in `online learning' in machine learning and its failure to reference the adaptive signal processing literature. Finally we speculate on future developments of LMS.

One of the main approaches to system identification of networks of time series or signals is conditional independence graphical (CIG) modeling. In the Gaussian case, the conditional dependence structure of the nodal time series is determined by the location of zeros in the precision matrix (inverse covariance matrix). And this determines the graph structure of the network. Despite the many applications of CIG models, the theory and algorithms have so far only dealt with networks of univariate or scalar signals. But in most applications the nodes carry multivariate or vector signals. Here we extend CIG modeling to handle such data by posing a group ℓ0sparse penalised block precision matrix estimation problem. We develop a double cyclic descent algorithm to solve it. And we compare the method with a group ℓ1penalised alternative in simulations.

Increasingly in the `big data' era, point process network data is appearing e.g. in social media but also genomics, high frequency finance and neurosicence. And the most basic problem is that of identifying network structure. As with all network system identification problems this is enabled by the use of sparsity. Here, for the first time, we develop a vector l0 sparsity approach for identification of a network of interacting Hawkes processes. We apply the new method to genomic data on a transcriptional regulatory network in embryonic stem cells and the results are compared with previous work and experimental findings.

Networks of interacting Hawkes processes have emerged as useful models in neuroscience, geophysics, high frequency finance, and social network analysis. The Hawkes process is of fundamental importance, being a point process analog of an autoregression. Here we develop a fixed gain adaptive (aka online) distributed estimator for the parameters of a Hawkes process model. The stochastic intensity is modeled by a causal Laguerre basis expansion. The natural recursive structure of this basis is exploited to derive a new two time scale adaptive algorithm based on exponentially weighted least squares which preserves non-negativity constraints. Simulations illustrate the results.

There is intense interest in fMRI research on whole-brain functional connectivity, and however, two fundamental issues are still unresolved: the impact of spatiotemporal data resolution (spatial parcellation and temporal sampling) and the impact of the network construction method on the reliability of functional brain networks. In particular, the impact of spatiotemporal data resolution on the resulting connectivity findings has not been sufficiently investigated. In fact, a number of studies have already observed that functional networks often give different conclusions across different parcellation scales. If the interpretations from functional networks are inconsistent across spatiotemporal scales, then the whole validity of the functional network paradigm is called into question. This paper investigates the consistency of resting state network structure when using different temporal sampling or spatial parcellation, or different methods for constructing the networks. To pursue this, we develop a novel network comparison framework based on persistent homology from a topological data analysis. We use the new network comparison tools to characterize the spatial and temporal scales under which consistent functional networks can be constructed. The methods are illustrated on Human Connectome Project data, showing that the DISCOH2 network construction method outperforms other approaches at most data spatiotemporal resolutions.

In recent years, driven by scientific and clinical concerns, there has been an increased interest in the analysis of functional brain networks. The goal of these analyses is to better understand how brain regions interact, how this depends upon experimental conditions and behavioral measures and how anomalies (disease) can be recognized. In this paper, we provide, first, a brief review of some of the main existing methods of functional brain network analysis. But rather than compare them, as a traditional review would do, instead, we draw attention to their significant limitations and blind spots. Then, second, relevant experts, sketch a number of emerging methods, which can break through these limitations. In particular we discuss five such methods. The first two, stochastic block models and exponential random graph models, provide an inferential basis for network analysis lacking in the exploratory graph analysis methods. The other three addresses: network comparison via persistent homology, time-varying connectivity that distinguishes sample fluctuations from neural fluctuations, and network system identification that draws inferential strength from temporal autocorrelation.

In many engineering applications the state of a dynamical system is modelled by a Stochastic Differential Equation (SDE) evolving in a “curved” (non-Euclidean) space such as the Stiefel manifold - the set of n × p real matrices with orthonormal columns, (n ≥ p). Due to the advances in computing power, the problem of state estimation can be efficiently addressed by the Particle Filter (PF). However, PF algorithms have to be completely reworked to handle the geometry, and the very few papers that properly deal with either the geometry or the stochastics of the problem are in the mathematics literature and are not accessible to an engineering audience. With this in mind and motivated by deterministic schemes on the Stiefel manifold, we give a direct accessible derivation of a novel PF algorithm for state estimation on the Stiefel manifold such that the resulting estimators always remain on the manifold. Our method can be applied to ANY dynamical system (SDE) on the Stiefel manifold. We do not rely on differential geometry or advanced stochastic calculus. Simulation examples are provided.

The Hawkes process is the workhorse of dynamic point process modelling - the point process version of the autoregression. It has been applied, for example, in high frequency finance, electricity price spike modelling and gene regulatory network modelling. But, in all these and other applications, it is assumed the parameters are time invariant. However, it is becoming clear that in many applications the parameters vary with time. Here, we develop for the first time, a very simple local likelihood approach to estimation of time-variant Hawkes processes. The new algorithm is tested on simulations and then applied to data from the Australian electricity market.

In many important engineering applications the state dynamics of a system are modelled by Stochastic Differential Equations (SDEs) evolving in non-Euclidean spaces such as matrix Lie groups. Due to the advances in computing power, the problem of state estimation can be efficiently addressed by the particle filtering method. This requires dealing with both the geometry and the stochastics of the problem. However, the very few papers that properly deal with either are in the mathematics literature and not accessible. The engineering literature is also small but plagued with problems. With this in mind, we give a direct accessible derivation of the particle filter algorithm for state estimation in matrix Lie groups. We do not rely on differential geometry or advanced stochastic calculus. Simulation examples are provided.

There has been significant interest in sparse inverse covariance estimation in areas such as statistics, machine learning, and signal processing. In this problem, the sparse inverse of a covariance matrix of a multivariate normal distribution is estimated. A Penalised Log-Likelihood (PLL) optimisation problem is solved to obtain the matrix estimator, where the penalty is responsible for inducing sparsity. The most natural sparsity promoting penalty is the non-convex l0function. Due to speed and memory limitations, the existing algorithms for dealing with the non-convex l0PLL problem are unable to be used in high dimensional settings. Here we address this issue by presenting a new block iterative approach for this problem, which can handle large-scale data sizes. Simulations demonstrate that our approach outperforms existing methods for this problem.

There is an achilles heel underlying the consensus literature. It has been known for some time that measurement noise causes the explosive growth of the consensus mode. Here we state the noise problem; characterise the behaviour of the noisy consensus system including its drift to ∞ critique the existing remedies and develop a new remedy.

Non-negative matrix factorization (NMF) has found use in fields such as remote sensing and computer vision where the signals of interest are usually non-negative. Data dimensions in these applications can be huge and traditional algorithms break down due to unachievable memory demands. One is then compelled to consider distributed algorithms. In this paper, we develop for the first time a distributed version of NMF using the alternating direction method of multipliers (ADMM) algorithm and dyadic cyclic descent. The algorithm is compared to well established variants of NMF using simulated data, and is also evaluated using real remote sensing hyperspectral data.

The recent interest in the dynamics of networks and the advent, across a range of applications, of measuring modalities that operate on different temporal scales have put the spotlight on some significant gaps in the theory of multivariate time series. Fundamental to the description of network dynamics is the direction of interaction between nodes, accompanied by a measure of the strength of such interactions. Granger causality and its associated frequency domain strength measures (GEMs) (due to Geweke) provide a framework for the formulation and analysis of these issues. In pursuing this setup, three significant unresolved issues emerge. First, computing GEMs involves computing submodels of vector time series models, for which reliable methods do not exist. Second, the impact of filtering on GEMs has never been definitively established. Third, the impact of downsampling on GEMs has never been established. In this work, using state-space methods, we resolve all these issues and illustrate the results with some simulations. Our analysis is motivated by some problems in (fMRI) brain imaging, to which we apply it, but it is of general applicability.

Recently, the availability of high-dimensional point process data in a growing number of application areas is driving the demand for analysis tools for such data. An extremely challenging yet important problem is the inference of causal relationships from network data or network topology identification. This problem has received little attention in the literature until very recently. Here we develop, perhaps for the first time, a distributed optimization algorithm for large-scale dynamic networks of interacting Hawkes processes. Genomic data are analyzed to construct a transcriptional regulatory network in embryonic stem cells.

We consider estimation of the angular velocity of a satellite undergoing noisy process dynamics. Modelling the attitude as a diffusion process evolving in SO(3), we propose an offline nonparametric estimator of the angular velocity over an interval based on noise free observations of the attitude. This is an ill-posed problem, and Tikhonov regularization is employed along with a Fourier-Galerkin method to solve the resulting stochastic Euler equation. A geometry preserving numerical scheme is proposed to simulate the attitude dynamics. Simulations are given which provide heuristics for choosing the regularization parameter in different SNR scenarios.

Three emerging applications are driving a renewed interest in vector point processes: neural coding, high frequency finance and genomics. This pressure has revealed a gross lack of models and system identification methods. In particular in at least the first two applications coincidences can occur i.e. more than one event can occur at the same time. Yet the models in common use exclude this possibility. In this paper we develop a class of time-varying vector Poisson models that allow coincident events and develop for the first time an hypothesis test for no coincidences. We show simulation results and an application to high frequency finance data.

High dimensional data is often modeled as a linear combination of a sparse component, a low-rank component, and noise. An example is a video sequence of a busy scene where the background is the low-rank part and the foreground, e.g. moving pedestrians, is the sparse part. Sparse and low rank (SLR) matrix decomposition is a recent method that estimates those components. In this paper we develop an l0based SLR method and an associated tuning parameter selection method based on the extended Bayesian information criterion (EBIC) method. In simulations the new algorithm is compared with state of the art algorithms from the literature.

Principal component analysis (PCA) is one of the most widely used methods in multivariate signal processing. An important problem is to select the number of principal components (PCs). In this paper we develop an automatic method for selecting the number of PCs based on Stein's unbiased risk estimator (SURE). In simulations the new method outperforms state of the art cross-validation methods.

We develop a new approach to functional brain connectivity analysis, which deals with four fundamental aspects of connectivity not previously jointly treated. These are: temporal correlation, spurious spatial correlation, sparsity, and network construction using trajectory (as opposed to marginal) Mutual Information. We call the new method Sparse Conditional Trajectory Mutual Information (SCoTMI). We demonstrate SCoTMI on simulated and real fMRI data, showing that SCoTMI gives more accurate and more repeatable detection of network links than competing network estimation methods.

In this paper, we describe a new method for solving the magnetoencephalography inverse problem: temporal vector ℓ0-penalized least squares (TV-L0LS). The method calculates maximally sparse current dipole magnitudes and directions via spatial ℓ0 regularization on a cortically-distributed source grid, while constraining the solution to be smooth with respect to time. We demonstrate the utility of this method on real and simulated data by comparison to existing methods.

Principal component analysis (PCA) is a widely used signal processing technique. Instead of performing PCA in the data space, we consider the problem of sparse PCA in a potentially higher dimensional latent space. To do so, we zero-out groups of variables using vector £o regularization. The estimation is based on the maximization of the penalized log-likelihood, for which we develop an efficient coupled expectation-maximization (EM) - minorization-maximization (MM) algorithm. For the special case when the latent- and observation space are identical, our method corresponds to an existing vector £o PCA method, which we verify using simulations. The proposed method can also be utilized for penalized linear regression and we use simulations to demonstrate superior estimation performance. As an example of a practical application, we use our method to localize cortical activity from magnetoencephalography (MEG) data.

Recently, there has been a growing interest in dynamic networks for understanding interactions and information flows. A fundamental problem is the identification of the links or the network topology. In comparison with its time series counterpart, the problem has received little attention in the point process literature. But with high-dimensional point process data becoming available in a number of application areas such as communication networks and neural coding, topology identification has become crucial for understanding the information flows. Here we discuss for the first time topology identification of a dynamic network of interacting Hawkes processes. Cortical recordings from cats are used to identify the interaction of neurons in the primary visual cortex.

Sparse component analysis (SCA) is a widely used method for solving the blind source separation problem. We develop a new cyclic descent algorithm for SCA based on a dyadic expansion. To select the associated tuning parameter a method based on the Bayesian information criterion is developed. In simulations the new algorithm is compared with state of the art algorithms from the literature.

Any regularization method requires the selection of a penalty parameter and many model selection criteria have been developed based on various discrepancy measures. Most of the attention has been focused on prediction mean squared error. In this paper we develop a model selection criterion based on regression parameter mean squared error via SURE (Stein's unbiased risk estimator). We then apply this to the l1penalized least squares problem with grouped variables on over-determined systems. Simulation results based on topology identification of a sparse network are presented to illustrate and compare with alternative model selection criteria.

Recently, a lot of attention has been given to penalized least squares problem formulations for sparse signal reconstruction in the presence of noise. The penalty is responsible for inducing sparsity, where the common choice used is the convex l1norm. While an l0penalty generates maximum sparsity it has been avoided due to lack of convexity. With the hope of gaining improved sparsity but more computational tractability there has been recent interest in the lqpenalty. In this paper we provide a novel cyclic descent algorithm for optimizing the lqpenalized least squares problem when 00, l1and lq, 0 <; q <; 1.

In recent years there has been a surge in the demand for analysis tools for multivariate point process data driven by work in neural coding and high frequency finance. In both these areas data volumes have become huge but few dimension reduction methods have been developed. Here we introduce a reduced rank model for the multivariate point process and provide a maximum likelihood estimator which we compute by an NMF type algorithm. However, the dependence on the point process history in the model implies our algorithm does not fit the traditional framework. The method is illustrated with a simulation and some data from cortical recordings from cats.

While noise free matrix completion has a considerable history recent interest has centered on a noisy version of the problem. We consider a nuclear norm penalised least squares formulation. And by applying the SURE method, we develop for the first time, an automatic procedure for selecting the penalty parameter. We illustrate its use with some simulation results.

Finding low rank nonnegative decomposition of multivariate data has many important applications in signal processing. A standard method is the nonnegative matrix factorization (NMF). In recent years, many algorithm have been proposed for NMF. However, an important problem that has not received as much attention is the selection of the rank of NMF. In this paper we develop a method for selecting the rank of NMF based on the Stein's unbiased risk estimator (SURE). In simulations we compare the method against crossvalidation. In addition we apply the method for selecting the rank of NMF for high dimensional hyperspectral data.

There has been a fast-growing demand for analysis tools for multivariate point-process data driven by work in neural coding and, more recently, high-frequency finance. Here we develop a true or exact (as opposed to one based on time binning) principal components analysis for preliminary processing of multivariate point processes. We provide a maximum likelihood estimator, an algorithm for maximization involving steepest ascent on two Stiefel manifolds, and novel constrained asymptotic analysis. The method is illustrated with a simulation and compared with a binning approach.

Multivariate regression is one of the most widely applied multivariate statistical methods with many uses across a range of disciplines. But the number of parameters increases exponentially with dimension and reduced-rank regression (RRR) is a well known approach to dimension reduction. But traditional RRR applies only to an overdetermined system. For increasingly common undetermined systems this issue can be managed by regularization, e.g., with a quadratic penalty. A significant problem is then the choice of the two tuning parameters: one discrete i.e., the rank; the other continuous i.e., the Tikhonov penalty parameter. In this paper we resolve this problem via Stein's unbiased risk estimator (SURE). We compare SURE to cross-validation and apply it on both simulated and real data sets.

Rank minimization problems, which consist of finding a matrix of minimum rank subject to linear constraints, have been proposed in many areas of engineering and science. A specific problem is the matrix completion problem in which a low rank data matrix is recovered from incomplete samples of its entries by solving a rank penalized least squares problem. The rank penalty is in fact the l0norm of the matrix singular values. A convex relaxation of this penalty is the commonly used l1norm of the matrix singular values. In this paper we bridge the gap between these two penalties and propose a simple method for solving the lq, q ∈ (0, 1), penalized least squares problem for matrix completion. We illustrate with simulations comparing our method to others in terms of solution quality.

The scalar sparse under-determined linear regression problem has had a rapid development with the multivariate version being of more recent interest. In this paper we pose a vector l0penalized multivariate regression problem to generate coefficient vectors with shared sparsity profile and then solve the problem with a new cyclic descent algorithm. We give optimality conditions and also discuss penalty parameter selection. Finally we present simulation results that compare our algorithm with alternatives.

The most common approach to attitude estimation involves direct estimation of the rotation matrix or its associated quaternion. This results in a difficult constrained estimation problem. Here instead we develop a new online algorithm that directly estimates the angular velocity which is unconstrained. The rotation matrix or quaternion is then easily obtained from the kinematics. We sketch an averaging analysis of stability of the new algorithm.

In this paper we present a novel model based sparse principal component analysis method based on the l0penalty. We develop an estimation method based on the generalized EM algorithm and iterative hard thresholding and an associated model selection method based on Bayesian information criterion (BIC). The method is compared to a previous sparse PCA method using both simulated data and DNA microarray data.

The standard modeling framework in functional magnetic resonance imaging (fMRI) is predicated on assumptions of linearity, time invariance and stationarity. These assumptions are rarely checked because doing so requires specialized software, although failure to do so can lead to bias and mistaken inference. Identifying model violations is an essential but largely neglected step in standard fMRI data analysis. Using Lagrange multiplier testing methods we have developed simple and efficient procedures for detecting model violations such as nonlinearity, nonstationarity and validity of the common double gamma specification for hemodynamic response. These procedures are computationally cheap and can easily be added to a conventional analysis. The test statistic is calculated at each voxel and displayed as a spatial anomaly map which shows regions where a model is violated. The methodology is illustrated with a large number of real data examples.

Graphical models are well established in providing compact conditional probability descriptions of complex multivariable interactions. In the Gaussian case, graphical models are determined by zeros in the precision or concentration matrix, i.e. the inverse of the covariance matrix. Hence, there has been much recent interest in sparse precision matrices in areas such as statistics, machine learning, computer vision, pattern recognition and signal processing. In this paper we propose a simple new algorithm for constructing a sparse estimator for the precision matrix from multivariate data where the sparsity is enforced by an l0penalty. We compare and test the quality of our method on a synthetic graphical model.

A standard fMRI experiment is structured around the assumption that onset of relevant neural activity occurs almost immediately after external stimulus. Introducing deliberate lengthy delays in cognition, e.g. in problem solving tasks, necessitates adjusting the analysis model to account for any possible temporal latency. This is usually a computationally intensive task. Using Lagrange Multipliers, we develop a fast method to detect such delay in the Time-to-Onset of BOLD signal. This is conceptually different to measuring the Time-to-Peak of the BOLD Hemodynamic Response. Additionally, our method does not require prior knowledge of the length of the delay.

Reduced rank regression (RRR) has found application in various fields of signal processing. In this paper we propose a novel extension of the RRR model which we call sparse variable reduced rank regression (svRRR). By using a vector l1penalty we remove variables completely from the RRR. The proposed estimation algorithm involves optimization on the Stiefel manifold and we illustrate it both on a simulated and a real functional magnetic resonance imaging (fMRI) data set.

Many problems in signal processing involve finding sparse solutions to linear systems of equations. The usual way of achieving this involves minimizing a mixed penalty function composed of a quadratic l2term and a sparse inducing l1term. Some existing algorithms for minimization include cyclic descent, gradient projection and iterative fixed point methods. Cojugate gradient is well known as a fast algorithm for linear quadratic problems. Here we develop a nonlinear conjugate gradient algorithm for the l1penalized least squares problem. This new method uses no line search and is found to be very stable. Description of its performance is provided as well as simulations to demonstrate convergence and comparison to another algorithm.

Since voxel time courses in functional magnetic resonance imaging (fMRI) are mostly produced from complex-valued data by taking the magnitudes, they obey Rician distributions, which can be approximated as Gaussian distributions only when signal-to-noise ratios (SNRs) are high. In this paper, we derive the asymptotic power of our recently developed activation detection statistic for Rician fMRI. The analysis shows that the asymptotic power is dependent only on the ratios of signal parameters to noise parameter of Rician distributed voxel time series, and allows us to better understand the nature of low SNRs in fMRI data analysis. Based on the power analysis, a more general and descriptive definition of SNR is provided than classical one.

Random Matrix Theory has generated tremendous interest in recent years, partly from powerful results developed for multi-user detection theory but also for growing applications in statistics, signal processing and econometrics. However the current theory has emphasized white noise data. In this paper we present for the first time some results applicable to temporally stationary processes.

In neural coding ultra fine electrodes are used to record spike trains from the brains (neurons) of awake animals such as rats,cats and monkeys. But analog singals such as EEG and local field potentials are also recorded. The aim of such experiments is to study the dynamics of the brain with potential development of true neural prosthetics. One problem that arises in extracting information from these recordings is joint modeling of both point processes and analog signals. Here we develop for the first time a test for independence between an analog signal and a point process.

The group Lasso is an extension of the Lasso or l1-penalised least squares procedure. It forces simultaneous zeroing of groups of variables and has already been applied to sparse component analysis and ill-conditioned inverse problems. In this paper we address the unresolved problem of threshold or penalty parameter selection.

Functional Magnetic Resonance (fMRI) data is most often analyzed using linear regression type methods that consider each voxel separately or by using exploratory methods such as Principal Component Analysis (PCA) or Independent Component Analysis (ICA). In this paper we introduce a model, which we call XnPCA, that combines regression with PCA. Unlike the linear regression methods XnPCA allows for non-stationary noise. Additionally, since XnPCA is based on the maximum likelihood framework the Bayesian information criterion (BIC) can be used for model selection and comparison. We compare XnPCA to a regression model commonly used in fMRI research using real data from a combined visual-motor experiment.

Sparse principal component analysis combines the idea of sparsity with principal component analysis (PCA). There are two kinds of sparse PCA; sparse loading PCA (slPCA) which keeps all the variables but zeroes out some of their loadings; and sparse variable PCA (svPCA) which removes whole variables by simultaneously zeroing out all the loadings on some variables. In this paper we propose a model based svPCA method based on the l0penalty. We compare the detection performance of the proposed method with other subset selection method using a simulated data set. Additionally, we apply the method on a real high dimensional functional magnetic resonance imaging (fMRI) data set.

Recently there has developed considerable interest in using sparseness with PCA. Almost all previous methods concentrate on zeroing out some loadings. Here we develop a new approach which zeros out whole variables automatically. We formulate a vector l1penalized PCA criterion and optimize it by steepest descent along geodesic on a Grassman manifold. This ensures that each step obeys PCA orthogonality as well as an invariance property of the criterion. We show in simulations that it outperforms a previous svPCA algorithm and apply it to a real high dimensional functional Magnetic Resonance Imaging (fMRI) data.

Space-time separability has been assumed and applied to most of data analysis methods in functional magnetic resonance imaging (FMRI). We developed a procedure for testing space and time separability in the framework of the parametric cepstrum. In this paper, the asymptotic power of the proposed space-time separability test is analyzed. The analysis shows two important properties of the proposed test. The asymptotic power function involves cepstral coefficients only in the non-separable region (parameters of interest). And the non-centrality parameter of the asymptotic power function is a scaled Euclidean metric between the logarithms of a non- separable power spectral density (PSD) and a separable PSD.

The CLEAN algorithm is one of the best known signal processing algorithms in (radio-)astronomy. It is essentially a deconvolution procedure and is used e.g. to reconstruct a sparse (star) brightness distribution from noisy ('dirty' - hence 'cleaning') observed data. In this paper we make a connection for the first time between CLEAN and l1-denoising. We show CLEAN does a crude version of l1-denoising and develop a modified algorithm with much improved behaviour.

Principal component analysis (PCA) is probably the best known method for dimensionality reduction. Perhaps the most important problem in PCA is to determine the number of principal components in a given data set, and in effect separate signal from noise in the data set. Many methods have been proposed to deal with this problem but almost all of them fail in the important practical case when the number of observations is comparable to the number of variables, i.e., the realm of random matrix theory (RMT). In this paper, we propose to use Stein's unbiased risk estimator (SURE) to estimate, with some assistance from RMT, the number of principal components. The method is applied on simulated data and compared to BIC and the Laplace method.

A main purpose of data analysis in functional magnetic resonance imaging (fMRI) is to determine which regions of the brain are activated by pre-specified temporal stimuli. In recent work, under the assumption of known spectra, we developed a detection statistic based on a spatially and temporally correlated noise model. In this paper, we implement the developed test statistic, which includes spatial and temporal whitening operators. For the estimation of spatial and temporal correlations, we use the parametric cepstral modeling, which allows dramatic reduction of computation in the model fitting and very simple methods to obtain spatiotemporal whitening operators. Model comparison and selection are discussed as well. We apply the developed techniques to a human dataset.

A goal in functional Magnetic Resonance Imaging (fMRI) data analysis is determining whether a certain region of brain is activated by presented temporal stimuli. Since the fMRI data is a sequence of images, spatiotemporal models are needed and spatially and temporally correlated noise plays a crucial role in the models. Until very recently, most attention has focussed on temporally correlated but spatially independent models. And spatial correlation has been dealt with in an ad hoc fashion. We develop, for the first time, a properly formulated true spatiotemporal detection statistic based on a spatially and temporally correlated noise model. Additionally, we develop a theoretical performance analysis method for comparing different test statistics through Asymptotic Relative Efficiency (ARE) for the first time in fMRI. We perform simulations for the comparison of new test statistic with a standard statistic as well.

A number of signal estimation problems are arising where a relatively low dimensional state is to be estimated from a high dimensional observation sequence. In previous work we have shown this leads to considerable simplification in the structure of optimal state estimators even in non-linear problems. In these and other state estimation problems there is a growing interest in the computation of mutual information between unobserved state and observed sequence. Here we show that the mutual information computation can be likewise considerably simplified

Multivariate methods such as principal component analysis (PCA) and independent component analysis (ICA) have been found to be useful in functional magnetic resonance imaging (fMRI) research. They are often able to decompose the fMRI data so that the researcher can associate their components to some biological processes of interest such as the brain response resulting from a stimulus. In this paper we develop a new smooth version of the PCA derived from a maximum likelihood framework. We are thus led to an unusual use of AIC, BIC namely to choose two (rather than one) parameters simultaneously; the number of principal components and the degree of smoothness. The algorithm is applied to real fMRI data

PCA has found use as an exploratory technique for fMRI analysis. However underlying it is an implicit model that while allowing temporal non-stationary covariance assumes the same covariance structure for all voxels. Here we relax this assumption for the first time by developing a version of PCA that allows the covariance structure to vary spatially. The new method is applied to real data and provides interesting new insight.

Support vector machines have become important in classification, biometrics, machine learning and pattern recognition. However, successful application requires selection of various tuning parameters such as kernel parameters and penalty or margin parameters. We apply a new technique for this problem which provides very simple structure for the automatic selector.

Recent availability of high-speed analog to digital converters (ADC) has enabled the development of digital radiometers for Earth remote sensing. In the unprotected C band, undesired radio frequency interferences (RFI) have been observed and mitigation efforts are currently being carried out in the research community. They mainly consist of detection and filtering in the time or frequency domain. Because RFI can be very small and mistaken as geophysical signals, an additional approach could rely on the inversion of a model describing internal noise, RFI and geophysical data processed through the radiometer. Such an approach demands a thorough understanding of receiver internal noise statistics. As a first step, we propose to model the system by taking advantage of the time series record available at the digital back end of a receiver. The experimental setup includes an X band benchtop radiometer serving as a standard analog gain chain, and various inputs fed into the radiometer front end. The digital output of an eight-bit ADC following the analog chain was recorded using a logic analyzer. An exploratory data analysis enabled us to select a priori valid data. Statistical signal processing techniques encompassed loglog relationships, classical spectral estimation (smoothed periodogram and windowed autocovariance), and autoregressive modeling.

In order to construct spatial activation plots from functional magnetic resonance imaging (fMRI) data, a complex spatio-temporal modeling problem must be solved. A crucial part of this process is the estimation of the hemodynamic response (HR) function, an impulse response relating the stimulus signal to the measured noisy response. The estimation of the HR is complicated by the presence of low frequency colored noise. The standard approach to modeling the HR is to use simple parametric models, although FIR models have been used. We offer two contributions. First, we pursue a nonparametric approach using orthonormal causal Laguerre polynomials which have become popular in the system identification literature. It also happens that the shape of the basis elements is similar to that of a typical HR. We thus expect to achieve a compact, and so bias reduced, and low noise representation of the HR. Additionally, we develop a procedure for providing confidence intervals for the whole HR function. This feature is completely lacking in all previous work.

Optical flow estimation algorithms such as the Lukas-Kanade method and Horn and Schunk method require selection of a tuning parameter. In the former case, it is a neighbourhood size, in the latter, a penalty parameter. Selection of these tuning parameters is difficult in general but has a profound effect on the results. Therefore, automatic methods of selection are of great interest. In previous work, we developed selection methods for the above algorithms. Now we develop a selection procedure for a robust version of the Lukas-Kanade method. This is a non-trivial task since the robust algorithm is nonlinear.

It is implicit in traditional discussions of linear or nonlinear state estimation filters that there is no relation specified between the dimension of the state and the observation vector dimension. If anything though, the state would often be thought to have higher dimension. But increasingly in practice problems are arising where the reverse is the case. In this paper we show that state estimation filters, such as the Kalman filter undergo a remarkable simplification in structure and computation when the observation dimension is much larger than the state dimension. Both linear and nonlinear cases (including point processes) are discussed.

In order to construct spatial activation plots from functional magnetic resonance imaging (fMRI) data, a complex spatio-temporal modeling problem must be solved. A crucial part of this process is the estimation of the hemodynamic response (HR) function, an impulse response relating the stimulus signal to the measured noisy response. The estimation of the HR is complicated by the presence of low frequency colored noise. The standard approach to modeling the HR is to use simple parametric models, although FIR models have been used. We pursue a nonparametric approach using orthonormal causal Laguerre polynomials which have become popular in the system identification literature. It also happens that the shape of the basis elements is similar to that of a typical HR. We thus expect to achieve a compact and so bias reduced and low noise representation of the HR. This is not the case in FIR modeling, because a low FIR order is unable to cover the whole length of the HR over its region of support while a high FIR order results in overestimation of signal and underestimation of noise leading to misleading interpretations.

Neural spike train decoding algorithms and techniques to compute Shannon mutual information are important methods for analyzing how neural systems represent biological signals. Decoding algorithms are also one of several strategies being used to design controls for brain-machine interfaces. Developing optimal strategies to design decoding algorithms and compute mutual information are therefore important problems in computational neuroscience. We present a general recursive filter decoding algorithm based on a point process model of individual neuron spiking activity and a linear stochastic state-space model of the biological signal. We derive from the algorithm new instantaneous estimates of the entropy, entropy rate, and the mutual information between the signal and the ensemble spiking activity. We assess the accuracy of the algorithm by computing, along with the decoding error, the true coverage probability of the approximate 0.95 confidence regions for the individual signal estimates. We illustrate the new algorithm by reanalyzing the position and ensemble neural spiking activity of CA1 hippocampal neurons from two rats foraging in an open circular environment. We compare the performance of this algorithm with a linear filter constructed by the widely used reverse correlation method. The median decoding error for Animal 1 (2) during 10 minutes of open foraging was 5.9 (5.5) cm, the median entropy was 6.9 (7.0) bits, the median information was 9.4 (9.4) bits, and the true coverage probability for 0.95 confidence regions was 0.67 (0.75) using 34 (32) neurons. These findings improve significantly on our previous results and suggest an integrated approach to dynamically reading neural codes, measuring their properties, and quantifying the accuracy with which encoded information is extracted.

Neural receptive fields are dynamic in that with experience, neurons change their spiking responses to relevant stimuli. To understand how neural systems adapt their representations of biological information, analyses of receptive field plasticity from experimental measurements are crucial. Adaptive signal processing, the well-established engineering discipline for characterizing the temporal evolution of system parameters, suggests a framework for studying the plasticity of receptive fields. We use the Bayes' rule Chapman-Kolmogorov paradigm with a linear state equation and point process observation models to derive adaptive filters appropriate for estimation from neural spike trains. We derive point process filter analogues of the Kalman filter, recursive least squares, and steepest-descent algorithms and describe the properties of these new filters. We illustrate our algorithms in two simulated data examples. The first is a study of slow and rapid evolution of spatial receptive fields in hippocampal neurons. The second is an adaptive decoding study in which a signal is decoded from ensemble neural spiking activity as the receptive fields of the neurons in the ensemble evolve. Our results provide a paradigm for adaptive estimation for point process observations and suggest a practical approach for constructing filtering algorithms to track neural receptive field dynamics on a millisecond timescale.

Optical flow estimation algorithms such as the Lukas-Kanade (1981) method and Horn and Schunk (1981) method require selection of a tuning parameter. In the former case a neighbourhood size, in the latter, a penalty parameter. Selection of these tuning parameters is difficult in general but has a profound effect on the results. So automatic methods of selection are of great interest. In previous work we have developed such methods based on white noise assumptions and here we show how to adjust for the effect of spatially correlated errors. These always occur in practice and can degrade the performance of white noise based procedures.

There is a considerable literature on non-parametric estimation of signals buried in coloured noise of known or finitely parameterised spectrum. And a considerable literature on non-parametric estimation of the spectrum of an observed coloured noise. Here we consider the bivariate problem of jointly, non-parametrically estimating a signal in coloured noise as well as the coloured noise spectrum. It is straightforward enough to develop estimators for each infinite dimensional parameter but what is much less clear is how to choose jointly the two required tuning parameters. We develop, apparently for the first time, a criterion that achieves this.

Within the context of a properly formulated spatio-temporal model for functional magnetic resonance imaging (fMRI), we develop a procedure for wavelet-based activation estimation that explicitly accounts for the temporally correlated 'brain' noise. This leads to a kind of resolution level - dependent wavelet thresholding, a novelty in the spatial setting. The method is illustrated on a visual A-B fMRI experiment and is seen to give significantly better results compared with those obtainable via a more conventional smoothing step involving preprocessing with a spatially homogeneous filter.

In previous work we developed, apparently for the first time, an automatic criterion to choose when to stop the iteration in anisotropic diffusion signal reconstruction. However the method was computationally very intensive. Here, using asymptotic perturbation methods we derive an approximation to the criterion that is quite accurate but can be implemented with almost no additional computation beyond that needed to construct the signal reconstruction itself.

Anisotropic diffusion has become a valuable tool for multiscale nonlinear image analysis for example in edge detection and segmentation. We develop, apparently for the first time, an automatic criterion to choose when to stop the iteration in anisotropic diffusion signal reconstruction.

We consider signal estimation for functional MRI studies on multiple subjects. There are two major issues; alignment or registration of images across subjects, and using the multisubject information to capture covariance information; we discuss only the latter. Capturing this covariance information properly can lead to great improvements in statistical efficiency beyond what simple averaging can offer as well as compact description of group features.

Gradient-based optical flow estimation methods typically do not take into account errors in the spatial derivative estimates. The presence of these errors causes an errors-in-variables (EIV) problem. Moreover, the use of finite difference methods to calculate these derivatives ensures that the errors are strongly correlated between pixels. Total least squares (TLS) has often been used to address this EIV problem. However, its application in this context is flawed as TLS implicitly assumes that the errors between neighborhood pixels are independent. In this paper, a new optical flow estimation method (EIVM) is formulated to properly treat the EIV problem in optical flow. EIVM is based on Sprent's (1966) procedure which allows the incorporation of a general EIV model in the estimation process. In EIVM, the neighborhood size acts as a smoothing parameter. Due to the weights in the EIVM objective function, the effect of changing the neighborhood size is more complex than in other local model methods such as Lucas and Kanade (1981). These weights, which are functions of the flow estimate, can alter the effective size and orientation of the neighborhood. In this paper, we also present a data-driven method for choosing the neighborhood size based on Stein's unbiased risk estimators (SURE).

In the last half decade, fast methods of magnetic resonance imaging have led to the possibility, for the first time, of non-invasive dynamic brain imaging. This has led to an explosion of work in the Neurosciences. From a signal processing viewpoint the problems are those of nonlinear spatio-temporal system identification. In this paper, we develop new methods of identification using novel spatial regularization. We also develop a new model comparison technique and use that to compare our method with existing techniques on some experimental data.

We consider total variation denoising in the presence of coloured noise. We show that noise can sometimes but not always cause major problems with reconstruction. We suggest that signal to noise ratio is not the crucial factor, rather the spectral content of the noise is crucial. We show how to modify the reconstruction algorithm to deal with coloured noise as well as how to do automatic tuning parameter selection.

Local methods have long been used to estimate optical flow by fitting measurements in a small neighbourhood to a simple model. What is less well known are procedures to choose the neighbourhood size, weights and model order. In this paper, we show that the choice of these local model tuning variables can have a significant effect on the flow estimate. The optimal choice of these variables will depend on the image content, the noise level and the type of motion in the sequence. Hence, the development of a data-driven selection method is important research goal. This paper presents such a procedure based on Stein's unbiased risk estimators (SURE).

Using information-theoretic methods we develop simple results quantifying a lower bound for minimax estimation, a kind of infinite-dimensional Cramer-Rao lower bound, for signal estimation in possibly nonlinear, ill-conditioned, inverse problems. Our results reduce calculation to a geometric computation based on a modulus of continuity and make explicit connections with results in the literature on deterministic ill-conditioned inverse problems. Several applications are discussed.

We apply a general procedure of the author to choose penalty parameters in total variation denoising. This is an automatic method of tuning parameter choice for total variation denoising. The method is computationally much simpler than cross-validation.

Motion fields of real image sequences are typically piecewise smooth with discontinuities at object boundaries. Further, because of the aperture problem, only the normal component of the optical flow can be observed. Direct application of wavelet shrinkage to the normal component flow may yield undesirable results due to the correlated noise affecting the wavelet coefficients. In this paper, we present a new technique for estimating optical flow based on L/sub 1/ regularisation. The resulting flow estimate typically has many zero wavelet coefficients, but unlike wavelet shrinkage some of the remaining coefficients are allowed to "grow". To highlight this difference, we have named our new technique wavelet zeroing. Additionally, we present a data-driven, sequence adaptive method for optimally choosing the penalty parameter.

We consider a model comparison based on a new model selection criterion. We treat fMRI as a spatio-temporal system identification problem and compare our model fitting method based on spatial regularization methods with the so-called statistical parametric map technique currently popular in the fMRI literature. We illustrate our results with data from the brain obtained during a combined visual and motor experiment.

Although still in practice, the use of total least squares (TLS) in optical flow estimation is unreliable. The TLS implicitly assumes that the error terms affecting the partial derivatives of the image intensities are independent. The usual methods for estimating the partial derivatives ensures that the errors are strongly correlated. Due to this correlation, an alternative method is required to treat the resulting errors-in-variables (EN) problem. We propose a new method for estimating optical flow based on Sprent's (1966, 1969) procedure. This method incorporates a general EIV model and provides a far simpler computational procedure than found in previous solutions.

This paper brings together two research issues in optical flow estimation: the selection of the smoothing parameter and the proper treatment of the correlated noise in the partial derivatives of the image intensity. In previous work, the authors have developed an optical flow estimation method which properly treats the errors-in-variables problem cause by the correlated noise in the spatio-temporal derivatives of the image intensity. This method requires the specification of a neighbourhood size to use when optical flow calculations are made. In this paper we present a data-driven method for choosing the optimal neighbourhood size.

The authors consider spatio-temporal modelling of functional MRI data from an inverse problems point of view. Most of the modelling to date has been on a pixel by pixel basis with no acknowledgment given to spatial smoothness when it exists. The authors discuss regularization methods that address that issue and illustrate their results with experimental data.

In optical flow estimation, an additional constraint to the constant brightness assumption is required to uniquely determine both components of the flow. Typically, these constraints impose a smoothness requirement on the flow estimate. Since the smoothness constraint may be inconsistent with the brightness constraint, a smoothing parameter is introduced to control the tradeoff between satisfying the requirements of both constraints. Previously, there have only been heuristic discussions on how to choose the smoothing parameter. We show that the choice of the smoothing parameter can have a significant effect on the flow estimate and present a data-driven method based on estimated risk to select the smoothing parameter in the Horn and Schunck (1981) algorithm.

Discusses modelling of functional magnetic resonance imaging (fMRI) data as an ill-conditioned inverse problem. The model involved is non-linear and requires both temporal and spatial regularization.

Regularisation methods for the solution of inverse problems are well known although the theoretical study of their performance especially in image processing contexts is not well advanced. What is also much less resolved is smoothing or penalty parameter estimation. We describe a general procedure for estimation of auxiliary finite dimensional parameters in ill-conditioned inverse problems. The method is applicable to nonlinear problems, involves no approximations but offers computational advantages over cross validation and maximum likelihood.

Presents a new variational approach to the problem of computed tomography reconstruction from sparse data. The author uses a Tikhonov regularisation (quite different from that of Louis (1985)) which deals without approximation with discrete or nonuniform grids.

We present a new variational approach to the problem of computed tomography reconstruction from sparse data. We use a Tikhonov regularisation (quite different from that of Louis [1985]) which deals without approximation with discrete or nonuniform grids. Our algorithm requires calculation of a Green's function on a finite region and we show how this can be done very efficiently computationally using the numerical/analytic boundary element method (BEM).

The effects of a nonunity attenuation factor on the performance of an adaptive time delay estimation algorithm are studied using both deterministic and stochastic averaging. Analysis shows that failing to estimate such an attenuation factor may degrade the performance of the algorithm. When the adaptive delay estimator is augmented with an attenuation factor estimator, larger adaptive gain is desired to estimate the time variant delay. When there is noise in the signal, first order stochastic averaging analysis reveals that, if the attenuation factor is estimated, the estimator is biased. It is nevertheless true that the time delay estimator is unbiased.>

A method is presented for parametric modeling of stationary random fields. The class of parametric models considered allows the most general elliptic field, and by linear constraints can include such special cases as isotropic, quarter plane, and separable fields. The technique, based on the cepstrum, has the great advantage of requiring only the use of fast Fourier transforms in the fitting process. Thus, unlike the fitting of two-dimensional autoregressions, no iteration is necessary. Other advantages are that any (Wiener) filters constructed from the fitted spectrum are guaranteed to be stable, and that the spectrum is guaranteed to be positive. Statistical tests for determining various special types of field from data are developed. The choice of model order is discussed as well.