Search

Alternating Diffusion (AD) is a commonly applied diffusion-based sensor fusion algorithm. While it has been successfully applied to various problems, its computational burden remains a limitation. Inspired by the landmark diffusion idea considered in the Robust and Scalable Embedding via Landmark Diffusion (ROSELAND), we propose a variation of AD, called Landmark AD (LAD), which captures the essence of AD while offering superior computational efficiency. We provide a series of theoretical analyses of LAD under the manifold setup and apply it to the automatic sleep stage annotation problem with two electroencephalogram channels to demonstrate its application.

We present a theoretical foundation regarding the boundedness of the t-SNE algorithm. t-SNE employs gradient descent iteration with Kullback-Leibler (KL) divergence as the objective function, aiming to identify a set of points that closely resemble the original data points in a high-dimensional space, minimizing KL divergence. Investigating t-SNE properties such as perplexity and affinity under a weak convergence assumption on the sampled dataset, we examine the behavior of points generated by t-SNE under continuous gradient flow. Demonstrating that points generated by t-SNE remain bounded, we leverage this insight to establish the existence of a minimizer for KL divergence.

In this manuscript, we propose an efficient manifold denoiser based on landmark diffusion and optimal shrinkage under the complicated high dimensional noise and compact manifold setup. It is flexible to handle several setups, including the high ambient space dimension with a manifold embedding that occupies a subspace of high or low dimensions, and the noise could be colored and dependent. A systematic comparison with other existing algorithms on both simulated and real datasets is provided. This manuscript is mainly algorithmic and we report several existing tools and numerical results. Theoretical guarantees and more comparisons will be reported in the official paper of this manuscript.

Objective: Sleep spindles contain crucial brain dynamics information. We introduce the novel non-linear time-frequency analysis tool 'Concentration of Frequency and Time' (ConceFT) to create an interpretable automated algorithm for sleep spindle annotation in EEG data and to measure spindle instantaneous frequencies (IFs). Methods: ConceFT effectively reduces stochastic EEG influence, enhancing spindle visibility in the time-frequency representation. Our automated spindle detection algorithm, ConceFT-Spindle (ConceFT-S), is compared to A7 (non-deep learning) and SUMO (deep learning) using Dream and MASS benchmark databases. We also quantify spindle IF dynamics. Results: ConceFT-S achieves F1 scores of 0.749 in Dream and 0.786 in MASS, which is equivalent to or surpass A7 and SUMO with statistical significance. We reveal that spindle IF is generally nonlinear. Conclusion: ConceFT offers an accurate, interpretable EEG-based sleep spindle detection algorithm and enables spindle IF quantification., Comment: 22 pages, 6 figures, and 1 table

We introduce a novel ridge detection algorithm for time-frequency (TF) analysis, particularly tailored for intricate nonstationary time series encompassing multiple non-sinusoidal oscillatory components. The algorithm is rooted in the distinctive geometric patterns that emerge in the TF domain due to such non-sinusoidal oscillations. We term this method \textit{shape-adaptive mode decomposition-based multiple harmonic ridge detection} (\textsf{SAMD-MHRD}). A swift implementation is available when supplementary information is at hand. We demonstrate the practical utility of \textsf{SAMD-MHRD} through its application to a real-world challenge. We employ it to devise a cutting-edge walking activity detection algorithm, leveraging accelerometer signals from an inertial measurement unit across diverse body locations of a moving subject.

Dealing with time series with missing values, including those afflicted by low quality or over-saturation, presents a significant signal processing challenge. The task of recovering these missing values, known as imputation, has led to the development of several algorithms. However, we have observed that the efficacy of these algorithms tends to diminish when the time series exhibit non-stationary oscillatory behavior. In this paper, we introduce a novel algorithm, coined Harmonic Level Interpolation (HaLI), which enhances the performance of existing imputation algorithms for oscillatory time series. After running any chosen imputation algorithm, HaLI leverages the harmonic decomposition based on the adaptive nonharmonic model of the initial imputation to improve the imputation accuracy for oscillatory time series. Experimental assessments conducted on synthetic and real signals consistently highlight that HaLI enhances the performance of existing imputation algorithms. The algorithm is made publicly available as a readily employable Matlab code for other researchers to use.

Objective: Breathing pattern variability (BPV), as a universal physiological feature, encodes rich health information. We aim to show that, a high-quality automatic sleep stage scoring based on a proper quantification of BPV extracting from the single airflow signal can be achieved. Methods: Topological data analysis (TDA) is applied to characterize BPV from the intrinsically nonstationary airflow signal, where the extracted features are used to train an automatic sleep stage scoring model using the XGBoost learner. The noise and artifacts commonly present in the airflow signal are recycled to enhance the performance of the trained system. The state-of-the-art approach is implemented for a comparison. Results: When applied to 30 whole night polysomnogram signals with standard annotations, the leave-one-subject-out cross-validation shows that the proposed features (overall accuracy 78.8\%$\pm$8.7\% and Cohen's kappa 0.56$\pm 0.15$) outperforms those considered in the state-of-the-art work (overall accuracy 75.0\%$\pm$9.6\% and Cohen's kappa 0.50$\pm 0.15$) when applied to automatically score wake, rapid eyeball movement (REM) and non-REM (NREM). The TDA features are shown to contain complementary information to the traditional features commonly used in the literature via examining the feature importance. The respiratory quality index is found to be essential in the trained system. Conclusion: The proposed TDA-assisted automatic annotation system can accurately distinguish wake, REM and NREM from the airflow signal. Significance: Since only one single airflow channel is needed and BPV is universal, the result suggests that the TDA-assisted signal processing has potential to be applied to other biomedical signals and homecare problems other than the sleep stage annotation.

A common method for estimating the Hessian operator from random samples on a low-dimensional manifold involves locally fitting a quadratic polynomial. Although widely used, it is unclear if this estimator introduces bias, especially in complex manifolds with boundaries and nonuniform sampling. Rigorous theoretical guarantees of its asymptotic behavior have been lacking. We show that, under mild conditions, this estimator asymptotically converges to the Hessian operator, with nonuniform sampling and curvature effects proving negligible, even near boundaries. Our analysis framework simplifies the intensive computations required for direct analysis.

The moving average of the complex modulus of the analytic wavelet transform provides a robust time-scale representation for signals to small time shifts and deformation. In this work, we derive the Wiener chaos expansion of this representation for stationary Gaussian processes by the Malliavin calculus and combinatorial techniques. The expansion allows us to obtain a lower bound for the Wasserstein distance between the time-scale representations of two long-range dependent Gaussian processes in terms of Hurst indices. Moreover, we apply the expansion to establish an upper bound for the smooth Wasserstein distance and the Kolmogorov distance between the distributions of a random vector derived from the time-scale representation and its normal counterpart. It is worth mentioning that the expansion consists of infinite Wiener chaos and the projection coefficients converge to zero slowly as the order of the Wiener chaos increases, and we provide a rational-decay upper bound for these distribution distances, the rate of which depends on the nonlinear transformation of the amplitude of the complex wavelet coefficients.

The scattering transform is a multilayered, wavelet-based transform initially introduced as a model of convolutional neural networks (CNNs) that has played a foundational role in our understanding of these networks' stability and invariance properties. Subsequently, there has been widespread interest in extending the success of CNNs to data sets with non-Euclidean structure, such as graphs and manifolds, leading to the emerging field of geometric deep learning. In order to improve our understanding of the architectures used in this new field, several papers have proposed generalizations of the scattering transform for non-Euclidean data structures such as undirected graphs and compact Riemannian manifolds without boundary. In this paper, we introduce a general, unified model for geometric scattering on measure spaces. Our proposed framework includes previous work on geometric scattering as special cases but also applies to more general settings such as directed graphs, signed graphs, and manifolds with boundary. We propose a new criterion that identifies to which groups a useful representation should be invariant and show that this criterion is sufficient to guarantee that the scattering transform has desirable stability and invariance properties. Additionally, we consider finite measure spaces that are obtained from randomly sampling an unknown manifold. We propose two methods for constructing a data-driven graph on which the associated graph scattering transform approximates the scattering transform on the underlying manifold. Moreover, we use a diffusion-maps based approach to prove quantitative estimates on the rate of convergence of one of these approximations as the number of sample points tends to infinity. Lastly, we showcase the utility of our method on spherical images, directed graphs, and on high-dimensional single-cell data.

Nonsinusoidal oscillatory signals are everywhere. In practice, the nonsinusoidal oscillatory pattern, modeled as a 1-periodic wave-shape function (WSF), might vary from cycle to cycle. When there are finite different WSFs, $s_1,\ldots,s_K$, so that the WSF jumps from one to another suddenly, the different WSFs and jumps encode useful information. We present an iterative warping and clustering algorithm to estimate $s_1,\ldots,s_K$ from a nonstationary oscillatory signal with time-varying amplitude and frequency, and hence the change points of the WSFs. The algorithm is a novel combination of time-frequency analysis, singular value decomposition entropy and vector spectral clustering. We demonstrate the efficiency of the proposed algorithm with simulated and real signals, including the voice signal, arterial blood pressure, electrocardiogram and accelerometer signal. Moreover, we provide a mathematical justification of the algorithm under the assumption that the amplitude and frequency of the signal are slowly time-varying and there are finite change points that model sudden changes from one wave-shape function to another one., Comment: 39 pages, 11 figures

We develop a data-driven optimal shrinkage algorithm for matrix denoising in the presence of high-dimensional noise with a separable covariance structure; that is, the noise is colored and dependent across samples. The algorithm, coined {\em extended OptShrink} (eOptShrink) depends on the asymptotic behavior of singular values and singular vectors of the random matrix associated with the noisy data. Based on the developed theory, including the sticking property of non-outlier singular values and delocalization of the non-outlier singular vectors associated with weak signals with a convergence rate, and the spectral behavior of outlier singular values and vectors, we develop three estimators, each of these has its own interest. First, we design a novel rank estimator, based on which we provide an estimator for the spectral distribution of the pure noise matrix, and hence the optimal shrinker called eOptShrink. In this algorithm we do not need to estimate the separable covariance structure of the noise. A theoretical guarantee of these estimators with a convergence rate is given. On the application side, in addition to a series of numerical simulations with a comparison with various state-of-the-art optimal shrinkage algorithms, we apply eOptShrink to extract maternal and fetal electrocardiograms from the single channel trans-abdominal maternal electrocardiogram., Comment: arXiv admin note: text overlap with arXiv:1905.13060 by other authors

The manifold scattering transform is a deep feature extractor for data defined on a Riemannian manifold. It is one of the first examples of extending convolutional neural network-like operators to general manifolds. The initial work on this model focused primarily on its theoretical stability and invariance properties but did not provide methods for its numerical implementation except in the case of two-dimensional surfaces with predefined meshes. In this work, we present practical schemes, based on the theory of diffusion maps, for implementing the manifold scattering transform to datasets arising in naturalistic systems, such as single cell genetics, where the data is a high-dimensional point cloud modeled as lying on a low-dimensional manifold. We show that our methods are effective for signal classification and manifold classification tasks., Comment: Accepted for publication in the TAG in DS Workshop at ICML. For subsequent theoretical guarantees, please see Section 6 of arXiv:2208.08561

Objective: Commercial and research-grade wearable devices have become increasingly popular over the past decade. Information extracted from devices using accelerometers is frequently summarized as ``number of steps" (commercial devices) or ``activity counts" (research-grade devices). Raw accelerometry data that can be easily extracted from accelerometers used in research, for instance ActiGraph GT3X+, are frequently discarded. Approach: Our primary goal is proposing an innovative use of the {\em de-shape synchrosqueezing transform} to analyze the raw accelerometry data recorded from a single sensor installed in different body locations, particularly the wrist, to extract {\em gait cadence} when a subject is walking. The proposed methodology is tested on data collected in a semi-controlled experiment with 32 participants walking on a one-kilometer predefined course. Walking was executed on a flat surface as well as on the stairs (up and down). Main Results: The cadences of walking on a flat surface, ascending stairs, and descending stairs, determined from the wrist sensor, are 1.98$\pm$0.15 Hz, 1.99$\pm$0.26 Hz, and 2.03$\pm$0.26 Hz respectively. The cadences are 1.98$\pm$0.14 Hz, 1.97$\pm$0.25 Hz, and 2.02$\pm$0.23 Hz, respectively if determined from the hip sensor, 1.98$\pm$0.14 Hz, 1.93$\pm$0.22 Hz and 2.06$\pm$0.24 Hz, respectively if determined from the left ankle sensor, and 1.98$\pm$0.14 Hz, 1.97$\pm$0.22 Hz, and 2.04$\pm$0.24 Hz, respectively if determined from the right ankle sensor. The difference is statistically significant indicating that the cadence is fastest while descending stairs and slowest when ascending stairs. Also, the standard deviation when the sensor is on the wrist is larger. These findings are in line with our expectations. Conclusion: We show that our proposed algorithm can extract the cadence with high accuracy, even when the sensor is placed on the wrist.

Multivariate time-series have become abundant in recent years, as many data-acquisition systems record information through multiple sensors simultaneously. In this paper, we assume the variables pertain to some geometry and present an operator-based approach for spatiotemporal analysis. Our approach combines three components that are often considered separately: (i) manifold learning for building operators representing the geometry of the variables, (ii) Riemannian geometry of symmetric positive-definite matrices for multiscale composition of operators corresponding to different time samples, and (iii) spectral analysis of the composite operators for extracting different dynamic modes. We propose a method that is analogous to the classical wavelet analysis, which we term Riemannian multi-resolution analysis (RMRA). We provide some theoretical results on the spectral analysis of the composite operators, and we demonstrate the proposed method on simulations and on real data., Comment: 48 pages, 13 figures

Diverse disciplines are interested in how the coordination of interacting agents' movements, emotions, and physiology over time impacts social behavior. Here, we describe a new multivariate procedure for automating the investigation of this kind of behaviorally-relevant "interactional synchrony", and introduce a novel interactional synchrony measure based on features of dynamic time warping (DTW) paths. We demonstrate that our DTW path-based measure of interactional synchrony between facial action units of two people interacting freely in a natural social interaction can be used to predict how much trust they will display in a subsequent Trust Game. We also show that our approach outperforms univariate head movement models, models that consider participants' facial action units independently, and models that use previously proposed synchrony or similarity measures. The insights of this work can be applied to any research question that aims to quantify the temporal coordination of multiple signals over time, but has immediate applications in psychology, medicine, and robotics.

We explore the finite dimensional distributions of the second-order scattering transform of a class of non-Gaussian processes when all the scaling parameters go to infinity simultaneously. For frequently used wavelets, we find a coupling rule for the scale parameters of the wavelet transform within the first and second layers such that the limit exists. The coupling rule is explicitly expressed in terms of the Hurst index of the long range dependent inputs. More importantly, the spectral density function of the limiting process can be explicitly expressed in terms of the Hurst index of the long-range dependent input process and the Fourier transform of the mother wavelet. To obtain these results, we first show that the scattering transform of a class of non-Gaussian processes can be approximated by the second-order scattering transform of Gaussian processes when the scale parameters are sufficiently large. Secondly, for each approximation process, which is still a non-Gaussian random process, we show that its doubling scaling limit exists by the moment method. The long-range dependent non-Gaussian model considered in this work is realized by real data, and assumptions are supported by numerical results., Comment: 38 pages, 8 figures

Analysis of signals with oscillatory modes with crossover instantaneous frequencies is a challenging problem in time series analysis. One way to handle this problem is lifting the 2-dimensional time-frequency representation to a 3-dimensional representation, called time-frequency-chirp rate (TFC) representation, by adding one extra chirp rate parameter so that crossover frequencies are disentangled in higher dimension. The chirplet transform is an algorithm for this lifting idea, which leads to a TFC representation. However, in practice, we found that it has a strong ``blurring'' effect in the chirp rate axis, which limits its application in real-world data. Moreover, to our knowledge, we have limited mathematical understanding of the chirplet transform in the literature. Motivated by the need for the real-world data analysis, in this paper, we propose the synchrosqueezed chirplet transform (SCT) that enhances the TFC representation given by the chirplet transform. The resulting concentrated TFC representation has high contrast so that one can better distinguish different modes with crossover instantaneous frequencies. The basic idea is to use the phase information in the chirplet transform to determine a reassignment rule that sharpens the TFC representation determined by the chirplet transform. We also analyze the chirplet transform and provide theoretical guarantees of SCT., Comment: Fixed typos and added Remark 2

In many real world oscillatory signals, the fundamental component of a signal $f(t)$ might be weak or does not exist. This makes it difficult to estimate the instantaneous frequency of the signal. Traditionally, researchers apply the rectification trick, working with $|f(t)|$ or $\mbox{ReLu}(f(t))$ instead, to enhance the fundamental component. This raises an interesting question: what type of nonlinear function $g:\mathbb{R} \rightarrow \mathbb{R}$ has the property that $g(f(t))$ has a more pronounced fundamental frequency? $g(t) = |t|$ and $g(t) = \mbox{ReLu}(t)$ seem to work well in practice; we propose a variant of $g(t) = 1/(1-|t|)$ and provide a theoretical guarantee. Several simulated signals and real signals are analyzed to demonstrate the performance of the proposed solution.

We study the behavior of two kernel based sensor fusion algorithms, nonparametric canonical correlation analysis (NCCA) and alternating diffusion (AD), under the nonnull setting that the clean datasets collected from two sensors are modeled by a common low dimensional manifold embedded in a high dimensional Euclidean space and the datasets are corrupted by high dimensional noise. We establish the asymptotic limits and convergence rates for the eigenvalues of the associated kernel matrices assuming that the sample dimension and sample size are comparably large, where NCCA and AD are conducted using the Gaussian kernel. It turns out that both the asymptotic limits and convergence rates depend on the signal-to-noise ratio (SNR) of each sensor and selected bandwidths. On one hand, we show that if NCCA and AD are directly applied to the noisy point clouds without any sanity check, it may generate artificial information that misleads scientists' interpretation. On the other hand, we prove that if the bandwidths are selected adequately, both NCCA and AD can be made robust to high dimensional noise when the SNRs are relatively large.

We consider the eigenvalue problem $Ax = \lambda x$ where $A \in \mathbb{R}^{n \times n}$ and the eigenvalue is also real $\lambda \in \mathbb{R}$. If we are given $A$, $\lambda$ and, additionally, the absolute value of the entries of $x$ (the vector $(|x_i|)_{i=1}^n$), is there a fast way to recover $x$? In particular, can this be done quicker than computing $x$ from scratch? This may be understood as a special case of the phase retrieval problem. We present a randomized algorithm which provably converges in expectation whenever $\lambda$ is a simple eigenvalue. The problem should become easier when $|\lambda|$ is large and we discuss another algorithm for that case as well.

Background: We investigated clinical information underneath the beat-to-beat fluctuation of the arterial blood pressure (ABP) waveform morphology. We proposed the Dynamical Diffusion Map algorithm (DDMap) to quantify the variability of morphology. The underlying physiology could be the compensatory mechanisms involving complex interactions between various physiological mechanisms to regulate the cardiovascular system. As a liver transplant surgery contains distinct periods, we investigated its clinical behavior in different surgical steps. Methods: Our study used DDmap algorithm, based on unsupervised manifold learning, to obtain a quantitative index for the beat-to-beat variability of morphology. We examined the correlation between the variability of ABP morphology and disease acuity as indicated by Model for End-Stage Liver Disease (MELD) scores, the postoperative laboratory data, and 4 early allograft failure (EAF) scores. Results: Among the 85 enrolled patients, the variability of morphology obtained during the presurgical phase was best correlated with MELD-Na scores. The neohepatic phase variability of morphology was associated with EAF scores as well as postoperative bilirubin levels, international normalized ratio, aspartate aminotransferase levels, and platelet count. Furthermore, variability of morphology presents more associations with the above clinical conditions than the common BP measures and their BP variability indices. Conclusions: The variability of morphology obtained during the presurgical phase is indicative of patient acuity, whereas those during the neohepatic phase are indicative of short-term surgical outcomes., Comment: 5 figures and 1 table

We analyze the scattering transform with the quadratic nonlinearity (STQN) of Gaussian processes without depth limitation. STQN is a nonlinear transform that involves a sequential interlacing convolution and nonlinear operators, which is motivated to model the deep convolutional neural network. We prove that with a proper normalization, the output of STQN converges to a chi-square process with one degree of freedom in the finite dimensional distribution sense, and we provide a total variation distance control of this convergence at each time that converges to zero at an exponential rate. To show these, we derive a recursive formula to represent the intricate nonlinearity of STQN by a linear combination of Wiener chaos, and then apply the Malliavin calculus and Stein's method to achieve the goal., Comment: 32 pages, 1 figure

What is the influence of chronic maternal prenatal stress (PS) on fetal iron homeostasis? In a prospective case-control study in 164 pregnant women, we show that cord blood transferrin saturation is lower in male stressed neonates. The total effect of PS exposure on fetal ferritin revealed a decrease of 15.4% compared with controls. Electrocardiogram-based Fetal Stress Index (FSI) identified affected fetuses non-invasively during the third trimester of gestation. FSI-based timely detection of fetuses affected by PS can support early individualized iron supplementation and neurodevelopmental follow-up to prevent long-term sequelae due to PS-exacerbated impairment of the iron homeostasis., Comment: FELICITy study; ClinicalTrials.gov registration number is NCT03389178

Time-frequency (TF) representations of time series are intrinsically subject to the boundary effects. As a result, the structures of signals that are highlighted by the representations are garbled when approaching the boundaries of the TF domain. In this paper, for the purpose of real-time TF information acquisition of nonstationary oscillatory time series, we propose a numerically efficient approach for the reduction of such boundary effects. The solution relies on an extension of the analyzed signal obtained by a forecasting technique. In the case of the study of a class of locally oscillating signals, we provide a theoretical guarantee of the performance of our approach. Following a numerical verification of the algorithmic performance of our approach, we validate it by implementing it on biomedical signals.

Motivated by analyzing complicated and non-stationary time series, we study a generalization of the scattering transform (ST) that includes broad neural activation functions, which is called neural activation ST (NAST). On the whole, NAST is a transform that comprises a sequence of ``neural processing units'', each of which applies a high pass filter to the input from the previous layer followed by a composition with a nonlinear function as the output to the next neuron. Here, the nonlinear function models how a neuron gets excited by the input signal. In addition to showing properties like non-expansion, horizontal translational invariability and insensitivity to local deformation, the statistical properties of the second order NAST of a Gaussian process with various dependence and (non-)stationarity structure and its interaction with the chosen high pass filters and activation functions are explored and central limit theorem (CLT) and non-CLT results are provided. Numerical simulations are also provided. The results explain how NAST processes complicated and non-stationary time series, and pave a way towards statistical inference based on NAST under the non-null case.

We systematically study the spectrum of kernel-based graph Laplacian (GL) constructed from high-dimensional and noisy random point cloud in the nonnull setup. The problem is motived by studying the model when the clean signal is sampled from a manifold that is embedded in a low-dimensional Euclidean subspace, and corrupted by high-dimensional noise. We quantify how the signal and noise interact over different regions of signal-to-noise ratio (SNR), and report the resulting peculiar spectral behavior of GL. In addition, we explore the impact of chosen kernel bandwidth on the spectrum of GL over different regions of SNR, which lead to an adaptive choice of kernel bandwidth that coincides with the common practice in real data. This result paves the way to a theoretical understanding of how practitioners apply GL when the dataset is noisy.

Kernelized Gram matrix $W$ constructed from data points $\{x_i\}_{i=1}^N$ as $W_{ij}= k_0( \frac{ \| x_i - x_j \|^2} {\sigma^2} )$ is widely used in graph-based geometric data analysis and unsupervised learning. An important question is how to choose the kernel bandwidth $\sigma$, and a common practice called self-tuned kernel adaptively sets a $\sigma_i$ at each point $x_i$ by the $k$-nearest neighbor (kNN) distance. When $x_i$'s are sampled from a $d$-dimensional manifold embedded in a possibly high-dimensional space, unlike with fixed-bandwidth kernels, theoretical results of graph Laplacian convergence with self-tuned kernels have been incomplete. This paper proves the convergence of graph Laplacian operator $L_N$ to manifold (weighted-)Laplacian for a new family of kNN self-tuned kernels $W^{(\alpha)}_{ij} = k_0( \frac{ \| x_i - x_j \|^2}{ \epsilon \hat{\rho}(x_i) \hat{\rho}(x_j)})/\hat{\rho}(x_i)^\alpha \hat{\rho}(x_j)^\alpha$, where $\hat{\rho}$ is the estimated bandwidth function {by kNN}, and the limiting operator is also parametrized by $\alpha$. When $\alpha = 1$, the limiting operator is the weighted manifold Laplacian $\Delta_p$. Specifically, we prove the point-wise convergence of $L_N f $ and convergence of the graph Dirichlet form with rates. Our analysis is based on first establishing a $C^0$ consistency for $\hat{\rho}$ which bounds the relative estimation error $|\hat{\rho} - \bar{\rho}|/\bar{\rho}$ uniformly with high probability, where $\bar{\rho} = p^{-1/d}$, and $p$ is the data density function. Our theoretical results reveal the advantage of self-tuned kernel over fixed-bandwidth kernel via smaller variance error in low-density regions. In the algorithm, no prior knowledge of $d$ or data density is needed. The theoretical results are supported by numerical experiments on simulated data and hand-written digit image data.

In nonparametric regression, it is common for the inputs to fall in a restricted subset of Euclidean space. Typical kernel-based methods that do not take into account the intrinsic geometry of the domain across which observations are collected may produce sub-optimal results. In this article, we focus on solving this problem in the context of Gaussian process (GP) models, proposing a new class of Graph Laplacian based GPs (GL-GPs), which learn a covariance that respects the geometry of the input domain. As the heat kernel is intractable computationally, we approximate the covariance using finitely-many eigenpairs of the Graph Laplacian (GL). The GL is constructed from a kernel which depends only on the Euclidean coordinates of the inputs. Hence, we can benefit from the full knowledge about the kernel to extend the covariance structure to newly arriving samples by a Nystr\"{o}m type extension. We provide substantial theoretical support for the GL-GP methodology, and illustrate performance gains in various applications., Comment: 38 pages, 15 figures

Modern time series are usually composed of multiple oscillatory components, with time-varying frequency and amplitude contaminated by noise. The signal processing mission is further challenged if each component has an oscillatory pattern, or the wave-shape function, far from a sinusoidal function, and the oscillatory pattern is even changing from time to time. In practice, if multiple components exist, it is desirable to robustly decompose the signal into each component for various purposes, and extract desired dynamics information. Such challenges have raised a significant amount of interest in the past decade, but a satisfactory solution is still lacking. We propose a novel {\em nonlinear regression scheme} to robustly decompose a signal into its constituting multiple oscillatory components with time-varying frequency, amplitude and wave-shape function. We coined the algorithm {\em shape-adaptive mode decomposition (SAMD)}. In addition to simulated signals, we apply SAMD to two physiological signals, impedance pneumography and electroencephalography. Comparison with existing solutions, including linear regression, recursive diffeomorphism-based regression and multiresolution mode decomposition, shows that our proposal can provide an accurate and meaningful decomposition with computational efficiency.

Click-evoked otoacoustic emissions (CEOAEs) are clinically used as an objective way to infer whether cochlear functions are normal. However, because the sound pressure level of CEOAEs is typically much lower than the background noise, it usually takes hundreds, if not thousands of repetitions to estimate the signal with sufficient accuracy. In this paper, we propose to improve the signal-to-noise ratio (SNR) of CEOAE signals within limited measurement time by optimal shrinkage (OS) in two different settings: the covariance-based OS (cOS) and the singular value decomposition (SVD)-based OS (sOS). By simulation and analyzing human CEOAE data, the cOS consistently reduced the noise and enhanced the SNR by 1 to 2 dB from a baseline method (BM) that is based on calculating the median. The sOS achieved an SNR enhancement of 2 to 3 dB in simulation, and demonstrated capability to enhance the SNR in real recordings when the SNR achieved by the BM was below 0 dB. An appealing property of OS is that it produces an estimate of every individual column of the signal matrix. This property makes it possible to investigate CEOAE dynamics across a longer period of time when the cochlear conditions are not strictly stationary.

Airflow signal encodes rich information about respiratory system. While the gold standard for measuring airflow is to use a spirometer with an occlusive seal, this is not practical for ambulatory monitoring of patients. Advances in sensor technology have made measurement of motion of the thorax and abdomen feasible with small inexpensive devices, but estimation of airflow from these time series is challenging. We propose to use the nonlinear-type time-frequency analysis tool, synchrosqueezing transform, to properly represent the thoracic and abdominal movement signals as the features, which are used to recover the airflow by the locally stationary Gaussian process. We show that, using a dataset that contains respiratory signals under normal sleep conditions, an accurate prediction can be achieved by fitting the proposed model in the feature space both in the intra- and inter-subject setups. We also apply our method to a more challenging case, where subjects under general anesthesia underwent transitions from pressure support to unassisted ventilation to further demonstrate the utility of the proposed method.

When analyzing modern machine learning algorithms, we may need to handle kernel density estimation (KDE) with intricate kernels that are not designed by the user and might even be irregular and asymmetric. To handle this emerging challenge, we provide a strong uniform consistency result with the $L^\infty$ convergence rate for KDE on Riemannian manifolds with Riemann integrable kernels (in the ambient Euclidean space). We also provide an $L^1$ consistency result for kernel density estimation on Riemannian manifolds with Lebesgue integrable kernels. The isotropic kernels considered in this paper are different from the kernels in the Vapnik-Chervonenkis class that are frequently considered in statistics society. We illustrate the difference when we apply them to estimate the probability density function. Moreover, we elaborate the delicate difference when the kernel is designed on the intrinsic manifold and on the ambient Euclidian space, both might be encountered in practice. At last, we prove the necessary and sufficient condition for an isotropic kernel to be Riemann integrable on a submanifold in the Euclidean space., Comment: 50 pages

Objective: We aimed to fuse the outputs of different electrocardiogram-derived respiration (EDR) algorithms to create one EDR signal that is of higher quality. Methods: We viewed each EDR algorithm as a software sensor that recorded breathing activity from a different vantage point, identified high-quality software sensors based on the respiratory signal quality index, aligned the highest-quality EDRs with a phase synchronization technique based on the graph connection Laplacian, and finally fused those aligned, high-quality EDRs. We refer to the output as the sync-ensembled EDR signal. The proposed algorithm was evaluated on two large-scale databases of whole-night polysomnograms. We evaluated the performance of the proposed algorithm using three respiratory signals recorded from different hardware sensors, and compared it with other existing EDR algorithms. A sensitivity analysis was carried out for a total of five cases: fusion by taking the mean of EDR signals, and the four cases of EDR signal alignment without and with synchronization and without and with signal quality selection. Results: The sync-ensembled EDR algorithm outperforms existing EDR algorithms when evaluated by the synchronized correlation ({\gamma}-score), optimal transport (OT) distance, and estimated average respiratory rate (EARR) score, all with statistical significance. The sensitivity analysis shows that the signal quality selection and EDR signal alignment are both critical for the performance, both with statistical significance. Conclusion: The sync-ensembled EDR provides robust respiratory information from electrocardiogram. Significance: Phase synchronization is not only theoretically rigorous but also practical to design a robust EDR.

Adaptive Local Iterative Filtering (ALIF) is a currently proposed novel time-frequency analysis tool. It has been empirically shown that ALIF is able to separate components and overcome the mode-mixing problem. However, so far its convergence is still an open problem, particularly for highly nonstationary signals, due to the fact that the kernel associated with ALIF is non-translational invariant, non-convolutional and non-symmetric. Our first contribution in this work is providing a convergence analysis of ALIF. From the practical perspective, ALIF depends on a robust frequencies estimator, based on which the decomposition can be achieved. Our second contribution is proposing a robust and adaptive decomposition method for noisy and nonstationary signals, which we coined the Synchrosqueezing Iterative Filtering Technique (SIFT). In SIFT, we apply the synchrosqueezing transform to estimate the instantaneous frequency, and then apply the ALIF to decompose a signal. We show numerically the ability of this new approach in handling highly nonstationary signals.

We consider detecting the evolutionary oscillatory pattern of a signal when it is contaminated by non-stationary noises with complexly time-varying data generating mechanism. A high-dimensional dense progressive periodogram test is proposed to accurately detect all oscillatory frequencies. A further phase-adjusted local change point detection algorithm is applied in the frequency domain to detect the locations at which the oscillatory pattern changes. Our method is shown to be able to detect all oscillatory frequencies and the corresponding change points within an accurate range with a prescribed probability asymptotically. This study is motivated by oscillatory frequency estimation and change point detection problems encountered in physiological time series analysis. An application to spindle detection and estimation in sleep EEG data is used to illustrate the usefulness of the proposed methodology. A Gaussian approximation scheme and an overlapping-block multiplier bootstrap methodology for sums of complex-valued high dimensional non-stationary time series without variance lower bounds are established, which could be of independent interest.

Motivated by analyzing complicated time series, nonlinear-type time-frequency analysis became an active research topic in the past decades. Those developed tools have been applied to various problems. In this article, we review those developed tools and summarize their applications to high-frequency biomedical signals.

To handle time series with complicated oscillatory structure, we propose a novel time-frequency (TF) analysis tool that fuses the short time Fourier transform (STFT) and periodic transform (PT). Since many time series oscillate with time-varying frequency, amplitude and non-sinusoidal oscillatory pattern, a direct application of PT or STFT might not be suitable. However, we show that by combining them in a proper way, we obtain a powerful TF analysis tool. We first combine the Ramanujan sums and $l_1$ penalization to implement the PT. We call the algorithm Ramanujan PT (RPT). The RPT is of its own interest for other applications, like analyzing short signal composed of components with integer periods, but that is not the focus of this paper. Second, the RPT is applied to modify the STFT and generate a novel TF representation of the complicated time series that faithfully reflect the instantaneous frequency information of each oscillatory components. We coin the proposed TF analysis the Ramanujan de-shape (RDS) and vectorized RDS (vRDS). In addition to showing some preliminary analysis results on complicated biomedical signals, we provide theoretical analysis about RPT. Specifically, we show that the RPT is robust to three commonly encountered noises, including envelop fluctuation, jitter and additive noise.

Background: Accurate detection of QRS complexes during mobile, ultra-long-term ECG monitoring is challenged by instances of high heart rate, dramatic and persistent changes in signal amplitude, and intermittent deformations in signal quality that arise due to subject motion, background noise, and misplacement of the ECG electrodes. Purpose: We propose a revised QRS detection algorithm which addresses the above-mentioned challenges. Methods and Results: Our proposed algorithm is based on a state-of-the-art algorithm after applying two key modifications. The first modification is implementing local estimates for the amplitude of the signal. The second modification is a mechanism by which the algorithm becomes adaptive to changes in heart rate. We validated our proposed algorithm against the state-of-the-art algorithm using short-term ECG recordings from eleven annotated databases available at Physionet, as well as four ultra-long-term (14-day) ECG recordings which were visually annotated at a central ECG core laboratory. On the database of ultra-long-term ECG recordings, our proposed algorithm showed a sensitivity of 99.90% and a positive predictive value of 99.73%. Meanwhile, the state-of-the-art QRS detection algorithm achieved a sensitivity of 99.30% and a positive predictive value of 99.68% on the same database. The numerical efficiency of our new algorithm was evident, as a 14-day recording sampled at 200 Hz was analyzed in approximately 157 seconds. Conclusions: We developed a new QRS detection algorithm. The efficiency and accuracy of our algorithm makes it a good fit for mobile health applications, ultra-long-term and pathological ECG recordings, and the batch processing of large ECG databases.

Given a connected set $\Omega_0 \subset \mathbb{R}^2$, define a sequence of sets $(\Omega_n)_{n=0}^{\infty}$ where $\Omega_{n+1}$ is the subset of $\Omega_n$ where the first eigenfunction of the (properly normalized) Neumann $p-$Laplacian $ -\Delta^{(p)} \phi = \lambda_1 |\phi|^{p-2} \phi$ is positive (or negative). For $p=1$, this is also referred to as the Ratio Cut of the domain. We conjecture that, unless $\Omega_0$ is an isosceles right triangle, these sets converge to the set of rectangles with eccentricity bounded by 2 in the Gromov-Hausdorff distance as long as they have a certain distance to the boundary $\partial \Omega_0$. We establish some aspects of this conjecture for $p=1$ where we prove that (1) the 1-Laplacian spectral cut of domains sufficiently close to rectangles of a given aspect ratio is a circular arc that is closer to flat than the original domain (leading eventually to quadrilaterals) and (2) quadrilaterals close to a rectangle of aspect ratio $2$ stay close to quadrilaterals and move closer to rectangles in a suitable metric. We also discuss some numerical aspects and pose many open questions., Comment: 35 pages, 16 figures, comments welcome!