text:
\chapter{Robust Estimation of Self-Exciting {Generalized Linear Models}}
\chaptermark{Robust Estimation of GLM's}

\label{chap:hawkes}

In this chapter, we close the gap in theory of compressed sensing for non-i.i.d. data by providing theoretical guarantees on stable estimation of self-exciting generalized linear models. We consider the problem of estimating self-exciting {generalized linear models} from limited binary observations, where the history of the process serves as the covariate. In doing so, we relax the assumptions of i.i.d. covariates and exact sparsity common in CS. Our results indicate that utilizing sparsity recovers important information about the intrinsic frequencies of such processes.  We analyze the performance of two classes of estimators, namely the $\ell_1$-regularized maximum likelihood and greedy estimators, for a canonical self-exciting process and characterize the sampling tradeoffs required for stable recovery in the non-asymptotic regime. Our results extend those of compressed sensing for linear and generalized linear models with i.i.d. covariates to those with highly inter-dependent covariates. We further provide simulation studies as well as application to real spiking data from the mouse's lateral geniculate nucleus and the ferret's retinal ganglion cells under different nonlinear forward models which agree with our theoretical predictions.

\section{Introduction}

The theory of compressed sensing (CS) has provided a novel framework for measuring and estimating statistical models governed by sparse underlying parameters {\cite{donoho2006compressed, candes2006compressive, candes2006stable, candes2008introduction, needell2009cosamp, bruckstein2009sparse}}. In particular, for linear models with random covariates and sparsity of the parameters, the CS theory provides sharp trade-offs between the number of measurement, sparsity, and estimation accuracy. Typical theoretical guarantees imply that when the number of random measurements are roughly proportional to sparsity, then stable recovery of these sparse models is possible.

Beyond those described by linear models, observations from binary phenomena form a large class of data in natural and social sciences. Their ubiquity in disciplines such as neuroscience, physiology, seismology, criminology, and finance has urged researchers to develop formal frameworks to model and analyze these data. In particular, the theory of point processes provides a statistical machinery for modeling and prediction of such phenomena. Traditionally, these models have been employed to predict the likelihood of self-exciting processes such as earthquake occurrences \cite{ogata1988statistical, vere1970stochastic}, but have recently found applications in several other areas. For instance, these models have been used to characterize heart-beat dynamics \cite{barbieri2005point, valenza2013point} and violence among gangs \cite{egesdal}. Self-exciting point process models have also found significant applications in analysis of neuronal data {\cite{brown2001analysis,smith2003estimating, brown2004multiple,  paninski2004maximum, truccolo2005point, paninski2007statistical, pillow2011model}.}

In particular, point process models provide a principled way to regress binary spiking data with respect to extrinsic stimuli and neural covariates, and thereby forming predictive statistical models for neural spiking activity. {Examples include place cells in the hippocampus \cite{brown2001analysis}, spectro-temporally tuned cells in the primary auditory cortex \cite{calabrese2011generalized}, and spontaneous retinal or thalamic neurons spiking under tuned intrinsic frequencies \cite{Borowska,liets2003spontaneous}. Self-exciting point processes have also been utilized in assessing the functional connectivity of neuronal ensembles \cite{kim2011granger,brown_func_conn}.} When fitted to neuronal data, these models exhibit three main features: first, the underlying parameters are nearly sparse or compressible  \cite{Brown_pp, brown_func_conn}; second, the covariates are often highly structured and correlated; and third, the input-output relation is highly nonlinear. Therefore, the theoretical guarantees of compressed sensing do not readily translate to prescriptions for point process estimation.

Estimation of these models is typically carried out by Maximum Likelihood (ML) or regularized ML estimation {in discrete time, where the process is viewed as a Generalized Linear Model (GLM). In order to adjust the regularization level, empirical methods such as cross-validation are typically employed \cite{brown_func_conn}.} In the signal processing and information theory literature, sparse signal recovery under Poisson statistics has been considered in \cite{poi_sp_del} with application to the analysis of ranking data. In \cite{cs_poi}, a similar setting has been studied, with motivation from imaging by photon-counting devices. Finally, in theoretical statistics, high-dimensional $M$-estimators with decomposable regularizers, such as the $\ell_1$-norm, have been studied for GLMs \cite{Negahban}.

A key underlying assumption in the existing theoretical analysis of estimating {GLMs} is the independence and identical distribution (i.i.d.) of covariates. This assumption does not hold for self-exciting processes, since the history of the process takes the role of the covariates. Nevertheless, regularized ML estimators show remarkable performance in fitting {GLMs} to neuronal data with history dependence and highly non-i.i.d. covariates.  In this chapter, we close this gap by presenting new results on robust estimation of compressible GLMs, relaxing the assumptions of i.i.d. covariates and exact sparsity common in CS.

In particular, we will consider a {canonical GLM} and will analyze two classes of estimators for its underlying parameters: the $\ell_1$-regularized maximum likelihood and greedy estimators. We will present theoretical guarantees that extend those of CS theory and characterize fundamental trade-offs between the number of measurements, model compressibility, and estimation error of GLMs in the non-asymptotic regime. Our results reveal that when the number of measurements scale sub-linearly with the product of the ambient dimension and a generalized measure of sparsity (modulo logarithmic factors), then stable recovery of the underlying models is possible, even though the covariates solely depend on the history of the process. We will further discuss the extensions of these results to more general classes of {GLMs}. Finally, we will present applications to simulated as well as real data from {two classes of neurons exhibiting spontaneous activity}, namely the mouse's lateral geniculate nucleus and the ferret's retinal ganglion cells, which agree with our theoretical predictions. Aside from their theoretical significance, our results are particularly important in light of the technological advances in neural prostheses, which require robust neuronal system identification based on compressed data acquisition.

The rest of the chapter is organized as follows: In Section \ref{prelim}, we present our notational conventions, preliminaries and problem formulation. In Section \ref{theoretical}, we discuss the estimation procedures and state the main theoretical results of this chapter. Section \ref{simulations} provides numerical simulations as well as application to real data. In Section \ref{discussions}, we discuss the implications of our results and outline future research directions. Finally, we present the proofs of the main theoretical results and give a brief background on relevant statistical tests in Appendices \ref{appprf}--\ref{appks}.

\section{Preliminaries and Problem Formulation} \label{prelim}

We {first} give a brief introduction to self-exciting GLMs (see \cite{Daley2007} for a detailed treatment).

We consider a sequence of observations in the form of binary spike trains obtained by discretizing continuous-time observations (e.g. electrophysiology recordings), using bins of length $\Delta$. We assume that  not more than one event fall into any given bin. In practice, this can always be achieved by choosing $\Delta$ small enough. The binary observation at bin $i$ is denoted by $x_i$. The observation sequence can be modeled as the outcome of conditionally independent Poisson or Bernoulli trials, with a spiking probability given by $\mathbb{P}(x_i = 1) =: \lambda_{i|H_i}$, where $\lambda_{i|H_i}$ is the spiking probability at bin $i$ given the history of the process $H_i$ up to bin $i$.

These models are widely-used in neural data analysis and are motivated by the continuous time point processes with history dependent conditional intensity functions \cite{Daley2007}. For instance, given the history of a continuous-time point process $H_t$ up to time $t$, a conditional intensity of $\lambda(t | H_t) = \lambda$ corresponds to the {homogeneous Poisson} process.  As another example, a conditional intensity of $\lambda(t | H_t) = \mu + \int_{-\infty}^t \theta(t - \tau) dN(\tau)$ corresponds to a process known as the Hawkes process \cite{Hawkes_orig} with base-line rate $\mu$ and history dependence kernel $\theta(\cdot)$. Under the assumption of the orderliness of a continuous-time point process, a discretized approximation to these processes can be obtained by binning the process by bins of length $\Delta$, and defining the spiking probability by $\lambda_i := \lambda(i \Delta | H_{i\Delta}) \Delta + o(\Delta)$. In this chapter, we consider discrete {random processes} characterized by the spiking probability $\lambda_{i|H_i}$, which are either inherently discrete or employed as an approximation to continuous-time point process models.

Throughout the rest of the chapter, we drop the dependence of $\lambda_{i|H_i}$ on $H_i$ to simplify notation, denote it by $\lambda_i$ and refer to it {as} spiking probability. Given the sequence of {binary} observed data {$\mathbf{x}_1^n$}, the negative log-likelihood function under the {Bernoulli} statistics can be expressed as:

\begin{equation}
\label{L_def}
\mathfrak{L}(\boldsymbol{\theta}) = -\frac{1}{n} \sum\limits_{i=1}^n \left \{ x_i\log\lambda_{i}+(1-x_i)\log(1-\lambda_{i}) \right \}.
\end{equation}
Another common likelihood model used in the analysis of neuronal spiking data corresponds to Poisson statistics \cite{Brown_pp}, for which the negative log-likelihood takes the following form:
\begin{equation}
\label{poi_ll}
\mathfrak{L}(\boldsymbol{\theta}) :=  -\frac{1}{n} \sum\limits_{i=1}^n \left \{ x_i\log\lambda_{i}-\lambda_{i} \right \}.
\end{equation}
Throughout the chapter, we will focus on binary observations governed by Bernoulli statistics, whose negative log-likelihood is given in Eq. (\ref{L_def}). In applications such as electrophysiology in which neuronal spiking activities are recorded at a high sampling rate, the binning size $\Delta$ is very small and the Bernoulli and Poisson statistics coincide.

When the discrete process is viewed as an {approximation} to a continuous-time process, these log-likelihood functions are known as the Jacod log-likelihood approximations \cite{Daley2007}. We will present our analysis for the negative log-likelihood given by (\ref{L_def}), but our results can be extended to other statistics including (\ref{poi_ll}) (See {the remarks of Section \ref{theoretical}}).

Throughout this chapter {${\mathbf{x}_{-p+1}^{n}}$} will be considered as the observed spiking sequence which will be used for estimation purposes. A popular class of models for $\lambda_i$ is given by GLMs. In its general form, a GLM consists of two main components: an observation model and an equation expressing some (possibly nonlinear) function of the observation mean as a \textit{linear} combination of the covariates. In neural systems, the covariates consist of external stimuli as well as the history of the process. Inspired by spontaneous neuronal activity, we consider \textit{fully} self-exciting processes, in which the covariates are only functions of the process history. As for a canonical {GLM} inspired by the Hawkes process, we consider a process for which the spiking probability is a \textit{linear} function of the process history:
\begin{equation}
\label{CIF_models}
\lambda_i := \mu + \boldsymbol{\theta}' {\mathbf{x}_{i-p}^{i-1}},
\end{equation}
where $\mu$ is a positive constant representing the base-line rate, and $\boldsymbol{\theta}=[\theta_1,\theta_2,\cdots,\theta_p]'$ is a parameter vector denoting the history dependence of the process. {We further assume that the process is non-degenerate, i.e., it will not terminate in an infinite sequence of zeros.} We refer to this {GLM, viewed as a random process,} as the \textit{canonical self-exciting process}. Other popular models in the computational neuroscience literature include the log-link model where $\lambda_i = \exp ( \mu + \boldsymbol{\theta}'{\mathbf{x}_{i-p}^{i-1}} )$ and the logistic-link model where $\lambda_i = \frac{\exp(\mu + \boldsymbol{\theta}'{\mathbf{x}_{i-p}^{i-1}})}{1+\exp (\mu + \boldsymbol{\theta}'{\mathbf{x}_{i-p}^{i-1}})}$. The parameter vector $\boldsymbol{\theta}$ can be thought of as the binary equivalent of autoregressive (AR) parameters in linear AR models.

When fitted to neuronal spiking data, the parameter vector $\boldsymbol{\theta}$ exhibits a degree of sparsity \cite{Brown_pp, brown_func_conn}, that is, only certain lags in the history have a significant contribution in determining the statistics of the process. These lags can be thought of as the preferred or intrinsic {delays} in the spontaneous response of a neuron. {To be more precise, for a sparsity level $s < p$, we denote by $\boldsymbol{\theta}_s$ the best $s$-term approximation to $\boldsymbol{\theta}$.

Finally, in this chapter, we are concerned with the compressed sensing regime where $n \ll p$, i.e., the observed data has a much smaller length than the ambient dimension of the parameter vector. The main estimation problem of this chapter is the following: \emph{given observations ${\mathbf{x}_{-p+1}^n}$ from the \textcolor{black}{canonical self-exciting process}, the goal is to estimate the unknown baseline rate $\mu$ and the $p$-dimensional $(s,\xi)$-compressible history dependence parameter vector $\boldsymbol{\theta}$ in a stable fashion (where the estimation error is controlled) when $n \ll p$.}

\section{Theoretical Results}\label{theoretical}
In this section, we consider two estimators for $\boldsymbol{\theta}$, namely, the $\ell_1$-regularized ML estimator and a greedy estimator, and present the main theoretical results of this chapter on the estimation error of these estimators. Note that when $\mu$ is not known, the following results can be applied to the augmented parameter vector $[\mu,\boldsymbol{\theta}']'$. We analyze the case of known $\mu$ for simplicity of presentation.

\subsection{$\ell_1$-Regularized ML Estimation}

The natural estimator for the parameter vector is the ML estimator, which is widely used in neuroscience \cite{Brown_pp}, which by virtue of (\ref{L_def}) is given by:
\begin{equation}
\label{ML_est_pp_L}
\widehat{\boldsymbol{\theta}}_{{\sf ML}}=\argmin\limits_{\boldsymbol{\theta}\in \boldsymbol{\Theta} } \mathfrak{L}(\boldsymbol{\theta}),
\end{equation}
where $\boldsymbol{\Theta}$ is the relaxed closed convex feasible region for which $0 \le \lambda_i \le 1$ given by the conditions:

\vspace{-1mm}
\begin{align}\label{eq:star}
\begin{tabular}{l} ${0 < \pi_{\min} \le \mu -\mathbf{1}' \boldsymbol{\theta}^-}$,\\
 ${\mu +\mathbf{1}' \boldsymbol{\theta}^+ \leq \pi_\max < 1/2}$,
 \end{tabular} {\tag{$\star$}}
\end{align}
for some constants $\pi_{\min}$ and $\pi_{\max}$. {This first inequality incurs minimal loss of generality, as $\pi_{\min}$ can be chosen to be arbitrarily small. The restriction of $\pi_{\max} < 1/2$ ensures that the process is fast mixing and has mainly been adopted for technical convenience. This assumption incurs some loss of generality, as it excludes processes for which the maximum spiking probability exceeds $1/2$. However, due to the low spiking probability of typical neuronal activity, this loss is tolerable for the applications of interest in this chapter (see Section \ref{simulations})}.

In the regime of interest when $n \ll p$, the ML estimator is ill-posed and is typically regularized with a smooth norm. In order to capture the compressibility of the parameters, we consider  the $\ell_1$-regularized ML estimator:
\begin{equation}
\label{sp_est_pp_L}
\widehat{\boldsymbol{\theta}}_{{\sf sp}}:=\argmin\limits_{\boldsymbol{\theta}\in \boldsymbol{\Theta}} \quad \mathfrak{L}(\boldsymbol{\theta})+ \gamma_n\|\boldsymbol{\theta}\|_1.
\end{equation}
where $\gamma_n > 0$ is a regularization parameter. It is easy to verify that {the objective function and constraints in Eq.} (\ref{sp_est_pp_L}) are convex in $\boldsymbol{\theta}$ and hence $\widehat{\boldsymbol{\theta}}_{\sf sp}$ can be obtained using standard numerical solvers. Note that the solution to (\ref{sp_est_pp_L}) might not be unique. However, we will provide error bounds that hold for all possible solutions of (\ref{sp_est_pp_L}), with high probability.

It is known that ML estimates are asymptotically unbiased under mild conditions, and with $p$ fixed, the solution converges to the true parameter vector as $n \rightarrow \infty$. However, it is not clear how fast the convergence rate is for finite $n$ or when $p$ is not fixed and is allowed to scale with $n$. This makes the analysis of ML estimators, and in general regularized M-estimators, very challenging \cite{Negahban}. Nevertheless, such an analysis has significant practical implications, as it will reveal sufficient conditions on $n$ with respect to $p$ as well as a criterion to choose $\gamma_n$, which result in a stable estimation of $\boldsymbol{\theta}$. Finally, note that we are fixing the ambient dimension $p$ throughout the analysis. In practice, the history dependence is typically negligible beyond a certain lag and hence for a large enough $p$, {GLMs} fit the data very well.

\subsection{Greedy Estimation}
Although there exist fast solvers to convex problems of the type given by Eq. (\ref{sp_est_pp_L}), these algorithms are polynomial time in $n$ and $p$, and may not scale well with high-dimensional data.  This motivates us to consider greedy solutions for the estimation of $\boldsymbol{\theta}$. In particular, we will consider a generalization of the Orthogonal Matching Pursuit (OMP) \cite{zhang_omp,OMP} for general convex cost functions. A flowchart of this algorithm is given in Table \ref{tab:1}, which we denote by the Point Process Orthogonal Matching Pursuit (POMP) algorithm. At each iteration, the component in which the objective function has the largest deviation is chosen and added to the current support. The algorithm proceeds for a total of $s^{\star}$ steps, resulting in an estimate with $s^\star$ components.

The main idea behind the generalized OMP is in the greedy selection stage, where the absolute value of the gradient of the cost function at the current solution is considered as the selection metric. Consider an estimate $\widehat{\boldsymbol{\theta}}^{(k-1)}$ at the $(k-1)$-st stage of the generalized OMP for a quadratic cost function of the form $\| \mathbf{b} - \mathbf{A} \boldsymbol{\theta}\|_2^2$, with $\mathbf{b}$ and $\mathbf{A}$ denoting the observation vector and covariates matrix, respectively. Then, the gradient takes the form $\mathbf{A}' (\mathbf{b} - \mathbf{A} \widehat{\boldsymbol{\theta}}^{(k-1)})$ which is exactly the correlation vector between the residual error and the columns of $\mathbf{A}$ as in the original OMP algorithm.

\begin{table}
\centering
\framebox{$\begin{array}{l}
\text{Input: } \mathfrak{L}(\boldsymbol{\theta}) , s^\star\\
\text{Output: } \widehat{\boldsymbol{\theta}}_{\sf POMP}^{(s^\star)}\\
\text{Initialization:}\Big\{\begin{array}{l}
\text{Start with the index set } S^{(0)}=\emptyset\\
\text{and the initial estimate }\widehat{\boldsymbol{\theta}}^{(0)}_{{\sf POMP}} = 0
\end{array}\\
\textbf{for } k=1,2,\cdots,s^\star\\
\text{  }\begin{array}{l}
j = \argmax \limits_i \left| \left( \nabla \mathfrak{L} \; \left(\widehat{\boldsymbol{\theta}}_{{\sf POMP}}^{(k-1)}\right) \right)_i\right|\\
S^{(k)}=S^{(k-1)}\cup \{j\}\\
\widehat{\boldsymbol{\theta}}_{{\sf POMP}}^{(k)} = \argmin \limits_{\support ({\boldsymbol{\theta}}) \subset S^{(k)}} \mathfrak{L}(\boldsymbol{\theta})
\end{array}\\
\textbf{end }\\
\end{array}$
}
\caption{\small{Point Process Orthogonal Matching Pursuit (POMP)}}
\label{tab:1}
\end{table}

\subsection{Theoretical Guarantees}

Recall that the parameter vector $\boldsymbol{\theta} \in \mathbb{R}^p$ is assumed to be $(s,\xi)$-compressible, so that $\sigma_s(\boldsymbol{\theta}) = \|\boldsymbol{\theta}-\boldsymbol{\theta}_S\|_1 = \mathcal{O} (s^{1-\frac{1}{\xi}})$, and the observed data are given by the vector ${\mathbf{x}_{-p+1}^n} \in \{0,1\}^{n+p-1}$, all in the regime of $s, n \ll p$. In the remainder of this chapter, we assume that $\boldsymbol{\theta} \in \boldsymbol{\Theta}$. The main theoretical result regarding the performance of the $\ell_1$-regularized ML estimator is given by the following theorem:
\begin{thm}
\label{negahban_lambda}
If $\sigma_s(\boldsymbol{\theta}) = \mathcal{O}(\sqrt{s})$, there exist constants $d_1,d_2,d_3$ and $d_4$ such that for $n>d_1s^{2/3}p^{2/3} \log p$ and a choice of $\gamma_n=d_2 \sqrt{\frac{\log p}{n}}$, any solution $\widehat{\boldsymbol{\theta}}_{{\sf sp}}$ to (\ref{sp_est_pp_L}) satisfies the bound
\begin{equation}
\label{thm2}
\left \|\widehat{\boldsymbol{\theta}}_{{\sf sp}}-\boldsymbol{\theta}\right \|_2 \leq d_3 \sqrt{\frac{s \log p}{n}}+ \sqrt{d_3\sigma_s(\boldsymbol{\theta})}\sqrt[4]{\frac{\log p}{n}},
\end{equation}
with probability greater than $1-\mathcal{O}\left(\frac{1}{n^{d_4}}\right)$.
\end{thm}

Similarly, the following theorem characterizes the performance bounds for the POMP estimate:

\begin{thm}
\label{thm_OMP}
If $\boldsymbol{\theta}$ is $(s,\xi)$-compressible for some $\xi < 1/2$, there exist constants $d_1',d_2',d_3'$ and $d_4'$ such that for $n>d_1's^{2/3} p^{2/3} \left(\log s \right)^{2/3} \log p$, the POMP estimate satisfies the bound
\begin{equation}
\label{thm2_OMP}
\left \|\widehat{\boldsymbol{\theta}}_{{\sf POMP}}-\boldsymbol{\theta}\right\|_2 \leq d_2' \sqrt{\frac{s \log s \log p}{n}} + d_3' \frac{\log s}{s^{{\frac{1}{\xi}-2}}}
\end{equation}
after $s^\star =\mathcal{O}(s \log s)$ iterations with probability greater than $1-\mathcal{O}\left(\frac{1}{n^{d'_4}}\right)$.
\end{thm}

{Full proofs of Theorems \ref{negahban_lambda} and \ref{thm_OMP} are given in Appendix \ref{appprf}.}

\noindent  \textit{\textbf{Remarks.}} An immediate comparison of the sufficient condition $n = \mathcal{O}(s^{2/3} p^{2/3} \log p)$ of Theorem \ref{negahban_lambda} with those of \cite{Negahban} for GLM models with i.i.d. covariates given by $n = \mathcal{O}(s \log p)$ reveals that a loss of order $\mathcal{O}(p^{2/3} s^{-1/3})$ is incurred due to the inter-dependence of the covariates. However, the sample space of $n$ i.i.d. covariates is $np$-dimensional, whereas in our problem the sample space is only $(n+p)$-dimensional. Hence, the aforementioned loss can be viewed as the price of self-averaging of the process accounting for the low-dimensional nature of the covariate sample space. \textcolor{black}{To the best of our knowledge, the dominant loss of $\mathcal{O}(p^{2/3})$ in both theorems does not seem to be significantly improvable, as self-exciting processes are known to converge quite slowly to their ergodic state \cite{reynaud2007some}. {On a related note, the analysis of the sampling requirements of linear AR models reveals a loss of {$\mathcal{O}(p^{1/2})$} in the number of measurements \cite{arpaper}.}}

The sufficient condition of Theorem \ref{thm_OMP} given by $n = \mathcal{O}(s^{2/3} p^{2/3} \left(\log s \right)^{2/3} \log p)$ implies an extra loss of $\left(\log s \right)^{2/3}$ due to the greedy nature of the solution. Theorem \ref{thm_OMP} also requires a high compressibility level of the parameter vector $\boldsymbol{\theta}$  ($\xi < 1/2$), whereas Theorem \ref{negahban_lambda} does not impose any extra restrictions on $\xi \in (0,1)$. Intuitively speaking, this comparison reveals the trade-off between computational complexity and compressibility requirements for convex optimization vs. greedy techniques, which is well-known for linear models \cite{bruckstein2009sparse}.

The constants $d_i, d'_i$, $i=1,\cdots,4$, $\alpha$ and $\beta$ are explicitly given in the proof of the theorems in Appendix \ref{appprf}. As for a typical numerical example, for $\pi_\min = 0.01$ and $\pi_\max=0.49$, the constants of Theorem \ref{negahban_lambda} can be chosen as $d_1 \approx 10^3, d_2 = 50, d_3 \approx 10^4$ and $d_4 = 4$.

The main ingredient in the proofs of Theorems \ref{negahban_lambda} and \ref{thm_OMP} is inspired by the beautiful treatment of Negahban et al. in \cite{Negahban} in establishing the notion of Restricted Strong Convexity (RSC).  The major technical challenge for the \textcolor{black}{canonical self-exciting process}, as opposed to the GLM models with i.i.d. covariates in \cite{Negahban}, lies in the fact that the covariates are highly inter-dependent as they are formed by the history of the process. Hence, it is not straightforward to establish RSC with high probability, as the large deviation techniques used for i.i.d. random vectors {do} not hold. We establish the RSC for the \textcolor{black}{canonical self-exciting process} in two steps (see Lemma \ref{lemma1} in Appendix \ref{appprf}). First, we show that RSC holds for the expected value of the negative log-likelihood $\mathbb{E} [ \mathfrak{L}(\boldsymbol{\theta})]$, and then by invoking results on concentration of dependent random variables show that the negative log-likelihood $\mathfrak{L}(\boldsymbol{\theta})$ resides in a sufficiently small neighborhood of $\mathbb{E} [ \mathfrak{L}(\boldsymbol{\theta})]$ with high probability, and hence satisfies the RSC.

The remainder of the proof of Theorem \ref{negahban_lambda} establishes that upon satisfying the RSC, the estimation error can be suitably bounded (Proposition \ref{negahban_thm}, Appendix \ref{appprf}). Similarly, Theorem 2 is proven using the RSC of the {canonical self-exciting process} together with the results adopted from \cite{zhang_omp} on the performance of OMP for convex cost functions (Proposition \ref{prop_omp}, Appendix \ref{appprf}).

\medskip

\noindent {\textbf{Extensions.}} For simplicity and clarity of presentation, we have opted to present the proofs for the case of known $\mu$ and for the {canonical self-exciting process} as a canonical {GLM.} The following corollary extends our results to the case of unknown $\mu$.

\begin{corollary}\label{cor:1} The claims of Theorems \ref{negahban_lambda} and \ref{thm_OMP} hold when $\mu$ is not known, except for possibly slightly different constants.
\end{corollary}
\begin{proof}
The proof is given in Appendix \ref{appprf}.
\end{proof}

The {canonical self-exciting process} can be generalized to a larger class of {GLMs} by generalization of its spiking probability function. In a more general form we can consider a spiking probability function given by
\begin{equation*}
\label{nhp}
\lambda_i = \phi\left(\mu+\boldsymbol{\theta}'{\mathbf{x}_{i-p}^{i-1}}\right),
\end{equation*}
where $\phi(\cdot)$ is a possibly nonlinear function for which $0 < \lambda_i < 1$. In their continuous form, such processes are referred to as the \textit{nonlinear Hawkes process} \cite{nlhp}. Two of the commonly-used models in neural data analysis are the log-link and logistic-link models. Our prior numerical studies in \cite{kazemipour} revealed a similar performance improvement of  the $\ell_1$-regularized ML and the greedy solution over the ML estimate for the log-link model. Stationarity of these discrete processes can be proved similar to the canonical self-exciting process (see Appendix \ref{app:hawkes_psd}). The latter fact is key to extending our proofs to other models and is summarized by the following corollary:

\begin{corollary}\label{cor:2} Theorems \ref{negahban_lambda} and \ref{thm_OMP} hold when the spiking probability is given by $\lambda_i = \phi\left(\mu+\boldsymbol{\theta}'{\mathbf{x}_{i-p}^{i-1}}\right)$ for some continuous, bounded, convex and twice-differentiable function $\phi(\cdot)$ (e.g., $\phi(x) = \exp(x)$ or $\phi(x) = {\sf logit}^{-1}(x)$)  for which $0 < \lambda_i < 1/2$, except for  different constants.
\end{corollary}
\begin{proof}
The proof is given in Appendix \ref{appprf}.
\end{proof}

\section{Application to Simulated and Real Data}\label{simulations}

In this section, we study the performance of the conventional ML estimator, the $\ell_1$-regularized ML estimator, and the POMP estimator on simulated data as well as real spiking data recorded from the mouse's lateral geniculate nucleus (LGN) neurons and the ferret's retinal ganglion cells (RGC). {We have archived a MATLAB implementation of the estimators used in this chapter using the CVX package \cite{cvx} on the open source repository GitHub and made it publicly available \cite{hawkes_code}.}

\vspace{-3mm}
\subsection{Simulation Studies}

In order to simulate spiking data governed by the canonical self-exciting process, we sequentially generate spikes using (\ref{CIF_models}).   We have used $\mu = 0.1$, $\pi_\min = 0.01$, $\pi_\max = 0.49$, $p=1000$, $s=3$ and $n = 950$ for simulation purposes. Figure \ref{hawkes_fig} shows $500$ samples of the {canonical self-exciting process} generated using a history dependence parameter vector shown in Figure \ref{estimate_plot_synthetic}(a). The parameter vector $\boldsymbol{\theta}$ is compressible with a sparsity level of $s = 3$ and $\sigma_3(\boldsymbol{\theta}) = 0.05$. A value of $\gamma_n = 0.1$ is used to obtain the $\ell_1$-regularized ML estimate, which is slightly tuned around the theoretical estimate given by Theorem \ref{negahban_lambda}. Figures \ref{estimate_plot_synthetic}(b), \ref{estimate_plot_synthetic}(c), and \ref{estimate_plot_synthetic}(d) show the estimated history dependence parameter vectors using ML, $\ell_1$-regularized ML, and POMP, respectively. It can be readily visually observed that regularized ML and POMP significantly outperform the ML estimate in finding the correct values of $\boldsymbol{\theta}$. More specifically, the components at lags $405$ and $800$ (indicated by the gray arrows) are underestimated by the ML estimator, and their contribution is distributed among several falsely identified smaller lag components.

\begin{figure}[H]
\vspace{-2mm}
\centering
\includegraphics[width=.8\columnwidth]{hawkes}
\caption{{\small{A sample of the simulated {canonical self-exciting process}.}}}\label{hawkes_fig}
\vspace{-3mm}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=.9\columnwidth]{synthetic_gt}
\subcaption*{(a) True}
\vspace{3mm}
\includegraphics[width=.9\columnwidth]{synthetic_ml}
\subcaption*{(b) ML}
\vspace{3mm}
\includegraphics[width=.9\columnwidth]{synthetic_sp}
\subcaption*{(c) $\ell_1$-regularized ML}
\vspace{3mm}
\includegraphics[width=.9\columnwidth]{synthetic_OMP}
\subcaption*{(d) POMP}
\end{center}
\vspace{-3mm}
\caption{{\small{(a) True parameters vs. (b) ML, (c) $\ell_1$-regularized ML, and (d) POMP estimates.}}}\label{estimate_plot_synthetic}
\vspace{-2mm}
\end{figure}

\begin{figure}[H]
\centering
\includegraphics[width=.9\columnwidth]{MSE_p500}
\caption{{\small{MSE performance of the ML, $\ell_1$-regularized ML and POMP estimators.}}}\label{MSE_high_dim}
\vspace{-5mm}
\end{figure}

In order to quantify this performance gain, we repeated the same experiment {by generating  realizations corresponding to randomly chosen supports of size $s=3$ for $\boldsymbol{\theta}$ and spike trains of length $10^2 \le n \le 10^6$. In each case, the magnitudes of the components of $\boldsymbol{\theta}$ were chosen to satisfy the assumptions (\ref{eq:star}). For a given $\boldsymbol{\theta}$, the mean-square-error (MSE) of the estimate $\widehat{\boldsymbol{\theta}}$ is defined as $\widehat{\mathbb{E}}\{ \| \widehat{\boldsymbol{\theta}} - \boldsymbol{\theta} \|_2^2 \}$, where $\widehat{\mathbb{E}} \{ \cdot \}$ is the sample average over the realizations of the process.} Figure \ref{MSE_high_dim} shows the results of this simulation, where a similar systematic performance gain is observed. The left segment of the plot (shaded in yellow) and the right segment correspond to the compressive ($n < p$) and denoising ($n > p$) regimes, respectively. Error bars on the plot indicate 90\% quantiles of the MSE for this simulation obtained by multiple realizations. As it can be inferred from Figure \ref{MSE_high_dim}, the $\ell_1$-regularized ML and POMP have a systematic performance gain over the ML estimate {in the compressive regime, where $n \ll p$}, with the former outperforming the rest. In the denoising regime, the performance of the $\ell_1$-regularized and ML become closer, while the POMP saturates to a higher MSE floor. The latter observation can be explained by the fact that the POMP can only estimate $s^{\star}$ components (including those of $\boldsymbol{\theta}_s$), and fails to capture the $(p - s^{\star})$ compressible components. This results in an MSE floor above that obtained by ML, for large values of $n$.

\begin{figure}[H]
\begin{center}
\includegraphics[width=.8\columnwidth]{ML_test}
\subcaption*{(a) ML}
\vspace{3mm}
\includegraphics[width=.8\columnwidth]{sp_test}
\subcaption*{(b) $\ell_1$-regularized ML}
\vspace{3mm}
\includegraphics[width=.8\columnwidth]{OMP_test}
\subcaption*{(c) POMP}
\end{center}
\vspace{-3mm}
\caption{{\small{KS and ACF tests at {$95\%$ confidence level}, for the ML, $\ell_1$-regularized ML and POMP estimates.}}}
\label{hawkes_ks_acf_synthetic}
\vspace{-5mm}
\end{figure}

{The MSE comparison in Figure \ref{MSE_high_dim}} requires one to know the true parameters. In practice, the true parameters are unknown, and statistical tests are typically used to assess the goodness-of-fit of the estimates to the observed data. We use the Kolmogorov-Smirnov (KS) test and the autocorrelation function (ACF) test to assess the goodness-of-fit. These tests are based on the time-rescaling theorem for point processes \cite{time_rescaling}, which states that if the time axis is rescaled using the {estimated} conditional intensity function of {the} inhomogeneous Poisson process, the resulting point process is a homogeneous Poisson process with unit rate. Thereby, one can test for the validity of the time-rescaling theorem via two statistical tests: the KS test reveals how close the empirical quantiles of the time-rescaled point process to the true quantiles of a unit rate Poisson process, and the ACF test reveals how close the ISI {distribution} of the time-rescaled process is to the true ISI distribution of a unit rate Poisson process. Details of these tests are given in Appendix \ref{appks}.

Figure \ref{hawkes_ks_acf_synthetic} shows the KS and ACF tests {at a $95\%$ confidence level} for the ML $\ell_1$-regularized ML, and the POMP estimates from Figure \ref{estimate_plot_synthetic}. The yellow shades mark the regions below the specified confidence levels. The ML estimate {fails to pass the KS test, while the regularized and POMP estimates pass both tests.}

\subsection{Application to Spontaneous Neuronal Spiking Activity}

\subsubsection{{\textbf{Background and motivation}}}
Early studies of spontaneous neuronal activity from the cat's cochlear nucleus \cite{gerstein1960approach} marked a significant breakthrough in computational neuroscience by going beyond the so-called Poisson hypothesis, by which single neurons were assumed to be firing according to homogeneous Poisson statistics. The diversity of the ISIs deduced from the spontaneous activity of the cochlear neurons led to the development of more sophisticated statistical models based on renewal process theory, resulting in the Gamma and inverse Gaussian ISI descriptions of spontaneous neuronal activity \cite{gerstein1964random,tuckwell2005introduction}. Due to the analytical difficulties involved in working with these models, their generalization to a broader range of spiking statistics is not straightforward.

In light of the more recent discoveries on the role of spontaneous neuronal activity in brain development \cite{xu2011instructive,blankenship2010mechanisms}, its relation to functional architecture \cite{tsodyks1999linking}, and its functional significance in a variety of modalities including retinal \cite{xu2011instructive}, visual \cite{luck1997neural}, auditory \cite{tritsch2007origin}, hippocampal \cite{sombati1995recurrent}, cerebellar \cite{aizenman1999regulation}, and thalamic \cite{pinault1998intracellular} function, the modeling and analysis of this phenomenon has sparked a renewed interest among researchers in recent years. In particular, models based on GLMs have shown to overcome the analytical difficulties of the abovementioned models based on renewal theory, and have been successfuly used in relating the spontaneous neuronal activity to instrinsic and extrinsic neural covariates \cite{paninski2004maximum,paninski2007statistical,time_rescaling,barbieri2001construction} as well as inferring the functional connectivity of neuronal ensembles \cite{kim2011granger,brown_func_conn}. The above-mentioned results rely on the accuracy of the ML estimation of these models. In addition, the estimated parameters are typically sparse. Therefore, the $\ell_1$-regularized ML and POMP estimators are expected to offer a more robust alternative than the ML, especially under the limited observation setting.

{In order to evaluate the performance of these estimators on real data, in the remainder of this section we will compare the performance of the ML, $\ell_1$-regularized ML, and POMP estimators in modeling the spontaneous spiking activity recorded from two different types of neurons, namely the mouse's lateral geniculate nucleus and the ferret's retinal ganglion cells.}

In the following analysis, the regularization parameter $\gamma_n$ was chosen using a two-fold cross-validation refinement around the value obtained from our theoretical results. {The length of the history components $p$ was chosen by first selecting a large enough $p$ as an upper bound for the expected correlation length of neuronal spontaneous activity (estimated as $\sim 1.5~\text{s}$), followed by reducing $p$ to the point where an increase in the history length  does not result in significantly detected history components.}

\subsubsection{{\textbf{Application to LGN spiking activity}}}
\label{nsp}
 We first compare the performance of the estimators on the LGN neurons.  The LGN is part of the thalamus in the brain, which acts as a relay from the retina to the primary visual cortex \cite{thalamus}. The data {were} recorded at $1ms$ resolution from the mouse LGN neurons using single-unit recording \cite{scholl2013emergence}. We used about $5$ seconds of data from one neuron for the analysis. In order to capture the history dependence governing the spontaneous spiking activity of the LGN neuron, we model the spiking probability using the \textcolor{black}{canonical self-exciting process} model with $p=100$ ($\Delta = 1ms$). {Figure \ref{real_fig} shows the spiking data used in the analysis.}

\begin{figure}[H]
\centering
\includegraphics[width=.8\columnwidth]{real-01}
\caption{\small{The LGN spiking data used in the analysis.}}\label{real_fig}
\vspace{-1mm}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=.8\columnwidth]{real_estimate_ml}
\subcaption*{(a) ML}
\vspace{-1mm}
\includegraphics[width=.8\columnwidth]{real_estimate_sp}
\subcaption*{(b) $\ell_1$-regularized ML}
\vspace{3mm}
\includegraphics[width=.8\columnwidth]{real_estimate_pomp}
\subcaption*{(c) POMP}
\end{center}
\vspace{-3mm}
\caption{\small{(a) ML, (b) $\ell_1$-regularized ML, and (c) POMP estimates of the LGN spiking parameters.}}
\label{lgn_mlvssp}
\vspace{-5mm}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=.9\columnwidth]{real_ks_ml}
\subcaption*{(a) ML}
\vspace{3mm}
\includegraphics[width=.9\columnwidth]{real_ks_sp}
\subcaption*{(b) $\ell_1$-regularized ML}
\vspace{3mm}
\includegraphics[width=.9\columnwidth]{real_ks_pomp}
\subcaption*{(c) POMP}
\end{center}
\vspace{-3mm}
\caption{\small{KS and ACF tests at $99\%$ confidence level, for the ML, $\ell_1$-regularized ML and POMP estimates.}}\label{lgn_ks_sp}
\vspace{-5mm}
\end{figure}

Figure \ref{lgn_mlvssp} shows the estimated history dependence parameter vectors using the three methods. Both the regularized ML (Figure \ref{lgn_mlvssp}(b)) and POMP (Figure \ref{lgn_mlvssp}(c)) estimates capture significant history dependence components around a lag of $90$--$95~\text{ms}$ (marked by the upward arrows). In \cite{Borowska}, an intrinsic neuronal oscillation frequency of around $10 \hertz$ has been reported in around $30\%$ of all classes of mouse retinal cells under experiment, using combined two-photon imaging and patch-clamp recording. Our results are indeed consistent with the above mentioned findings about the intrinsic spiking frequency of retinal neurons. \textcolor{black}{To see this, we consider the power spectral density of the \textcolor{black}{canonical self-exciting process} given by:
\begin{equation}
\label{bart_spec}
S(\omega) = \frac{1}{2\pi} \left( \pi_\star^2 \delta(\omega) + \frac{\pi_\star-\pi_\star^2}{\left(1-\mathbf{1}'\boldsymbol{\theta}\right)^2 \left|1 - \Theta(\omega)\right|^2} \right),
\end{equation}
where $\Theta(\omega)$ is the discrete-time Fourier transform of $\boldsymbol{\theta}$ and $\pi_\star = \mu/(1-\mathbf{1}'\boldsymbol{\theta})$ denotes the stationary distribution probability of spiking. The derivation of the power spectral density is given in Appendix \ref{app:hawkes_psd}. The power spectral density of the \textcolor{black}{canonical self-exciting process} resembles the Bartlett spectrum of the Hawkes process \cite{Bartlett1, Bartlett2, Hawkes_orig}, whose peaks correspond to the significant oscillatory components of the underlying process.} Our estimated parameter vectors $\boldsymbol{\theta}$ using the regularized ML and POMP have significant nonzero components around lags of $ 90 \le k \le 95$. As a result, $S(\omega)$ peaks at $\omega = \frac{2\pi}{k\Delta}$. Hence, $f=\frac{1}{k\Delta}$ is an estimate of the significant intrinsic frequency of the underlying self-exciting process. Using the estimated numerical values, the intrinsic frequency is around $10.5$--$11~\text{Hz}$, which is consistent with experimental findings of \cite{Borowska}. Compared to the method in \cite{Borowska}, our estimates are obtained using much shorter recordings of spiking activity and provide a principled framework to study the oscillatory behavior of LGN neurons using {sparse GLM estimation.}

Note that there is a difference in the orders of magnitudes of the POMP estimate compared to the ML and regularized ML estimates. This is due to the fact that the POMP estimate is exactly $s$-sparse, whereas the ML and regularized ML estimates consist of $p = 100$ non-zero values. In order to assess the goodness-of-fit of these estimates, we invoke the KS and ACF tests. Figure \ref{lgn_ks_sp} shows the corresponding KS and ACF test plots. As it is implied from Figure \ref{lgn_ks_sp}(a), the ML estimate fails both tests due to overfitting, whereas the regularized ML (Figure \ref{lgn_ks_sp}(b)) passes both tests at the specified confidence levels. The POMP estimate (Figure \ref{lgn_ks_sp}(c)), however, passes the KS test while marginally failing the ACF test. The latter observation implies that the seemingly negligible components of the parameter vector captured by the regularized ML estimate seem to be important in explaining the statistics of the observed data.

\subsubsection{{\textbf{Application to RGC spiking activity}}} \label{sec:rgc}
{We will next study the performance of the estimators on spiking} data recorded from the RGCs of neonatal and adult ferrets \cite{wong1993transient}. The retinal ganglion cells are located in the innermost layer of the retina. They integrate information from photoreceptors and project them into the brain \cite{bear2007neuroscience}. The data {were} recorded using a multi-electrode array from the ferret retina at $50~\mu s$ \cite{wong1993transient}. We used $5$ seconds of data from one neuron for the analysis (neuron 2, session 1, adult data set, CARMEN data base \cite{eglen2014data}). {Figure \ref{real_2_fig} shows a segment of the spiking data used in our analysis. The RGC activity in the adult ferret is characterized by bursts of activity with a mean firing rate of $9 \pm 7$~Hz, which are separated by $0.5$--$1~\text{s}$ intervals \cite{wong1993transient}.}

 {
\begin{figure}[H]
\centering
\includegraphics[width=.95\columnwidth]{real_2-01}
\caption{\small{Segment of the RGC spiking data used in the analysis.}}\label{real_2_fig}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=.8\columnwidth]{RGC_ML_linear_estimates}
\subcaption*{(a) ML}
\includegraphics[width=.8\columnwidth]{RGC_sp_linear_estimates}
\subcaption*{(b) $\ell_1$-regularized ML}
\includegraphics[width=.8\columnwidth]{RGC_OMP_linear_estimates}
\subcaption*{(c) POMP}
\end{center}
\caption{\small{(a) ML, (b) $\ell_1$-regularized ML, and (c) POMP estimates of the RGC spiking parameters {using the canonical self-exciting process model.}}}
\label{rgc_mlvssp_linear}
\end{figure}

In order to capture the history dependence governing the spontaneous spiking activity of the RGC neuron, we model the spiking probability using two different link models to further corroborate the generalization of our results to models beyond the canonical self-exciting process studied in this chapter. First, we consider the canonical self-exciting process model. We have chosen  $\pi_\max = 0.49$,  $p=50$ ($\Delta = 25~\text{ms}$) and $s_\star = 3$. The baseline parameter $\mu$ is estimated from the data and is set to be equal to half of empirical mean firing rate of the neuron.

\begin{figure}[H]
\begin{center}
\includegraphics[width=.8\columnwidth]{RGC_ML_test}
\subcaption*{(a) ML}
\includegraphics[width=.8\columnwidth]{RGC_sp_test}
\subcaption*{(b) $\ell_1$-regularized ML}
\includegraphics[width=.8\columnwidth]{RGC_OMP_test}
\subcaption*{(c) POMP}
\end{center}
\caption{\small{KS and ACF tests at $95\%$ confidence level, for the ML, $\ell_1$-regularized ML and POMP estimates {using the canonical self-exciting process model.} }}\label{rgc_ks_sp_linear}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=.8\columnwidth]{rgc_exp_ml}
\subcaption*{(a) ML}
\includegraphics[width=.8\columnwidth]{rgc_exp_sp}
\subcaption*{(b) $\ell_1$-regularized ML}
\includegraphics[width=.8\columnwidth]{rgc_exp_pomp}
\subcaption*{(c) POMP}
\end{center}
\caption{\small{(a) ML, (b) $\ell_1$-regularized ML, and (c) POMP estimates of the RGC spiking parameters  {using the logistic link model}.}}
\label{rgc_mlvssp}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[width=.8\columnwidth]{rgc_exp_rates_ml}
\subcaption*{(a) ML}
\includegraphics[width=.8\columnwidth]{rgc_exp_rates_sp}
\subcaption*{(b) $\ell_1$-regularized ML}
\includegraphics[width=.8\columnwidth]{rgc_exp_rates_OMP}
\subcaption*{(c) POMP}
\end{center}
\caption{\small{(a) ML, (b) $\ell_1$-regularized ML, and (c) POMP estimates of the RGC spiking {probability using the unconstrained logistic link model}. Blue vertical lines show the locations of the spikes, and red traces show the estimated {probabilities}.}}
\label{rgc_rates}
\end{figure}

Figure \ref{rgc_mlvssp_linear} shows the estimated history components using the three estimators. {All three estimates} capture significant self-exciting history dependence components around the lags of $150~\text{ms}$ and $0.65$--$0.75~\text{s}$ (marked by the upward arrows). Invoking the foregoing argument for the LGN neuron regarding the power spectral density of the process (\ref{bart_spec}), these estimated lag components are consistent with the empirical estimates of \cite{wong1993transient}, as they indicate that the data can be characterized by a combination of $\frac{1}{150~\text{ms}} = 6.66$~Hz bursts separated by gaps of length $0.65$--$0.75~\text{s}$.
{The ML estimator predicts an extra self-inhibitory (negative) component, which results in over-fitting the data.}  This phenomenon can be observed by noting that the ML estimate fails the KS test shown in Figure \ref{rgc_ks_sp_linear}.

We will next consider a logistic link model of the form $\lambda_i = \frac{\exp(\mu + \boldsymbol{\theta}'x_{i-p}^{i-1})}{C+\exp (\mu + \boldsymbol{\theta}' {\mathbf{x}_{i-p}^{i-1}})}$, with $C = 100$. This model is widely used in neuronal modeling literature (e.g., \cite{truccolo2005point, Brown_pp}), where the assumptions given by (\ref{eq:star}) are dropped and the optimization is performed in an unconstrained fashion. We adopt this approach and obtain all the estimates by dropping the assumptions of $(\star)$. Figure \ref{rgc_mlvssp} shows the estimated history components {using the unconstrained estimators. } {Compared to the canonical self-exciting process model with a linear link, both the regularized ML (Figure \ref{rgc_mlvssp}(b)) and POMP (Figure \ref{rgc_mlvssp}(c)) estimates capture similar significant self-exciting history dependence components, which are consistent across the two sets of estimates.}

{The KS and ACF test results for this case are very similar to Figure \ref{rgc_ks_sp_linear} are are thus omitted for brevity.} In order to further inspect the goodness-of-fit of these methods, we plot the estimated spiking probabilities in Figure \ref{rgc_rates}. The ML estimate shown in Figure \ref{rgc_rates}(a) overfits the spiking events by rapidly saturating the rate to either 0 and 1, which results in undesired high rate estimates where there are no spikes. On the contrary, the regularized ML (Figure \ref{rgc_rates}(b)) and POMP (Figure \ref{rgc_rates}(c)) provide a more reliable estimate of the rates consistent with the spiking events. {This analysis suggests that the sufficient assumptions of $(\star)$ are not necessary for the superior performance of the regularized and POMP estimators over that of ML.}

\section{Concluding Remarks}\label{discussions}
In this chapter, we studied the sampling properties of $\ell_1$-regularized ML and greedy estimators for a \textcolor{black}{canonical self-exciting process}. The main theorems provide non-asymptotic sampling bounds on the number of measurements, which lead to stable recovery of the parameters of the process. To the best of our knowledge, our results are the first of this kind, and can be readily generalized to various other classes of self-exciting {GLMs}, such as processes with logarithmic or logistic links.

Compared to the existing literature, our results bring about two major contributions. First, we provide a theoretical underpinning for the advantage of $\ell_1$-regularization in ML estimation as well as greedy estimation in problems involving {binary} observations. These methods have been used in neuroscience in an ad-hoc fashion. Our results establish the utility of these techniques by characterizing the underlying sampling trade-offs. Second, our analysis relaxes the widely-assumed hypotheses of i.i.d. covariates. This assumption is often violated when working with history-dependent data such as neural spiking data.

We also verified the validity of our theoretical results through simulation studies as well as application to real neuronal spiking data {from mouse's LGN and ferret's RGC neurons.} These results show that both the regularized ML and the greedy estimates significantly outperform the widely-used ML estimate. In particular, through making a connection with the spectrum of discrete point processes, we were able to quantify the estimation of the intrinsic firing frequency of LGN neurons. {In the spirit of easing reproducibility, we have archived a MATLAB implementation of the estimators studied in this work using the CVX package \cite{cvx} on the open source repository GitHub and made it publicly available \cite{hawkes_code}.}

One of the limitations of our analysis is the assumption that the spiking probabilities are bounded by $1/2$, which results in loss of generality. This assumption is made for the sake of theoretical analysis in bounding the mixing rate of the canonical self-exciting process. Our numerical experiments suggest that it is not necessary for the operation of the $\ell_1$-regularized and POMP estimators. We consider further inspection of the mixing properties of this process and thus relaxing this assumption as future work. Our future work also includes generalization of our analysis to multivariate GLMs, which will allow to infer network properties from multi-unit recordings of neuronal ensembles.


TITLE: Is the image of this strip equal to the modular curve $Y(2)$
QUESTION [1 upvotes]: Let $B$ be the image of the strip $$\{x+iy : -1\leq x< 1, y>\frac{1}{2} \} \subset \mathbf{H}$$ in the modular curve $Y(2)$ under the quotient map $\mathbf{H} \to Y(2)$.
Is $B=Y(2)$?
Is $B=Y(2)$ iff the strip $\{x+iy : -1\leq x< 1, y>\frac{1}{2}\} $ contains a fundamental domain for $\Gamma(2)$?

REPLY [1 votes]: No, your strip does not surject onto $Y(2)$.
Domains $A \subseteq \mathbf{H}$ have a $\operatorname{PSL}(2,\mathbb{R})$-invariant "area" given by
$$ \int_{A} \frac{\mathrm{d}x\, \mathrm{d}y}{y^2} $$
Your strip has area 4, but the standard fundamental domain for $\operatorname{PSL}(2, \mathbb{Z})$ has area $\frac{\pi}{3}$. As $\Gamma(2)$ has index 6 in $\operatorname{PSL}(2,\mathbb{Z})$ and $6 \times \frac{\pi}{3} > 4$, there's no way it can surject onto $Y(2)$.
(PS: Here's another reason: $\Gamma(2)$ has 3 cusps, and your set contains a neighbourhood of $\infty$ but not of any other cusp. So there's no way that any set contained in $\{ z: \operatorname{Im}(z) > r\}$, for any $r > 0$, can surject onto $Y(2)$.)


TITLE: U-substitution step to solve Integral
QUESTION [0 upvotes]: I'm looking at how an integral was computed in "Paul's online math notes" in the section on U-substitution and am a bit confused.  In the original expression there are $2$ $x$'s, one in front of $\cos$ and the other as the numerator in the second term.  When he substitutes $xdx=\dfrac{1}{2}du$ why does $\dfrac{1}{2}du$ replace both $x$ terms?
thanks

REPLY [1 votes]: They are factoring the $x$ out of both terms so there will only be one $x$ in front of parenthesis.
It is just a matter of convenience. You do not have to factor if you do not like it. Then of course you have  two $ xdx$ to deal with.


TITLE: Calculating the radius of convergence for Taylor series
QUESTION [2 upvotes]: $$g(z)=\dfrac {\sin z-z+\dfrac{z^{3}}{6}}{\cos z-1}$$
What is the radius of convergence of the Taylor series of $g$ centered at 0.
My thought was to use the Cauchy-Hadamard formula to calculate it. But I don't know how to expand the $g(x)$ into power series.

REPLY [2 votes]: In this Wikipedia article it is stated that

The radius of convergence of a power series $f$ centered on a point $a$ is equal to the distance from $a$ to the nearest point where $f$ cannot be defined in a way that makes it holomorphic.

This can be proven using Cauchy's Integral Formula.
Since the closest non-removable singularity of $g$ to $z=0$ is $z=2\pi$, the radius of convergence should be $2\pi$.


TITLE: Show that there are no even natural numbers between $\frac{\pi}{10^{-k}}$ and $\frac{\pi}{\arctan (10^{-k})}$ for integer $k \ge 2$
QUESTION [2 upvotes]: Show that there are no even natural numbers between $\frac{\pi}{10^{-k}}$ and $\frac{\pi}{\arctan (10^{-k})}$ for integer $k \ge 2$. I encountered this while solving a physics problem and I am not even sure if it is true or not.
I tried to bound the difference of the two using the fact that $x > \arctan (x)$ for $x > 0$, but still there can be cases like $\frac{\pi}{10^{-k}} = n + 0.999999999$ for some integer $n$, regardless of how small the bound is.
Can anyone help?
EDIT : I have found Galperin's paper and it seems it is still left as a conjecture. All sources are in the comments of the first answer.

REPLY [2 votes]: We have $$ \frac{\pi}{\arctan (x)}-\frac{\pi}{x}=\frac{\pi}{x}\left(\frac{x}{\arctan(x)}-1\right)=\frac{\pi}{x}\left(\frac{x^2}{3}+O(x^4)\right)=\frac{\pi}{3}x + O(x^3) $$
Then the difference is around $\frac{\pi}{3}x \approx 10^{-k}$
To find an integer inbetween, then, roughly amounts to find, after $k$ digits in the decimal expansion of $\pi$, a sequence of at least $k$ consecutive nines, for some $k$.
This seems very unlikely. But also very difficult to prove. I don't think probabilistic approaches, assuming $\pi$ is normal, would lead anywhere.
Edit: as pointed out in the comments, the paper "Playing pool with π (the number π from a billiard point of view)" studies this (sec 10, question 1). It basically agrees with the previous paragraphs (the property "holds if and only if the string of $2N$ first decimal
digits of $\pi$ contains $N − 1$ nines in its right half"). It concedes that "modern mathematics is powerless to answer it" but it also conjectures -with almost certainity- that it's not true. Some probabilistic analysis is provided.


TITLE: How do i prove $(1-\frac{t}{n})^n < e^{-t}$ when $0<t<n$?
QUESTION [1 upvotes]: How do i prove : $(1-\frac{t}{n})^n < e^{-t}$ when $0<t<n$?
It seems obvious but i am absurdly struggling this for an hour.. Please help!
It's trivial that $\lim_{x\to\infty} (1-\frac{t}{x})^x = e^{-t}$. So it could be proven if it is monotonically increasing. But even the derivative is hard one know its positivity and negativity..

REPLY [1 votes]: Extending and detailing kmitov's answer:
Getting $\ln$ from both sides gives:
$$n\ln(1-\frac tn)<-t$$
Dividing by $n$ and putting $x=\frac t n$ gives:
$$\ln(1-x)<-x,\space 0<x<1$$
Or $f(x)=\ln(1-x)+x<0$.
The derivative of $f(x)$ is $1-\frac{1}{1-x}$ which is negative for $0<x<1$, so the function is decreasing. $f(0)=\ln(1)+0=0$. So $f(x)$ is negative for $0<x<1$ and the inequality holds as needed.


text:
\begin{document}

	\title[Analogues of cyclic insertion type identities for multiple zeta star values]{Analogues of cyclic insertion type identities \\ for multiple zeta star values}
	\date{June 27, 2018}
	\author{Steven Charlton}
	\address{Max-Planck-Institut f\"ur Mathematik \\
		Vivatsgasse 7 \\
		53111 Bonn \\
		Germany}

	\email{charlton@mpim-bonn.mpg.de}

	\keywords{Multiple zeta star values, cyclic insertion, block decomposition, generalised 2-1 formula}
	\subjclass[2010]{Primary 11M32}

	\begin{abstract}
		We prove an identity for multiple zeta star values, which generalises some identities due to Imatomi, Tanaka, Tasaka and Wakabayashi.  This identity gives an analogue of cyclic insertion type identities, for multiple zeta star values, and connects the block decomposition with Zhao's generalised 2-1 formula.
	\end{abstract}

	\maketitle

	\section{Introduction}

	Multiple zeta values (MZV's), and multiple zeta star values (MZSV's) are defined by the following series respectively
	\begin{align*}
		\zeta(s_1, \ldots, s_r) &\coloneqq \sum_{0 < n_1 < n_2 < \cdots < n_r} \frac{1}{n_1^{s_1} \cdots n_r^{s_r}} \, , \\
		\zeta^\star(s_1, \ldots, s_r) &\coloneqq \sum_{0 < n_1 \leq n_2 \leq \cdots \leq n_r} \frac{1}{n_1^{s_1} \cdots n_r^{s_r}} \, .
	\end{align*}
	In each case \( r \) is called the \emph{depth}, and \( s_1 + \cdots + s_r \) is called the \emph{weight}.  We make use of the shorthand `\( \wt \)' for the weight of the MZV's appearing in an identity, and write \( \{s\}^n \) to mean \( s \) repeated \( n \) times.  \medskip

	In \cite{imatomi2009some}, the following identities for multiple zeta star values are conjectured
	\begin{Conj}[Conjectures 4.1 and 4.3, \cite{imatomi2009some}]\label{conj:imatomi}
		For any integers \( a_0, \ldots, a_{2n} \geq 0 \) , we have
		\begin{gather*}
			\sum_{\text{permute \( a_0, \ldots, a_{2n} \)}} \zeta^\star(\{2\}^{a_0 + 1}, 1, \{2\}^{a_1}, 3, \{2\}^{a_2}, \ldots, 1, \{2\}^{a_{2n-1}}, 3, \{2\}^{a_{2n}}) \overset{?}{\in} \Q \pi^{\wt} \\
			\sum_{\text{permute \( a_1, \ldots, a_{2n} \)}} \zeta^\star(1, \{2\}^{a_1}, 3, \{2\}^{a_2}, \ldots, 1, \{2\}^{a_{2n-1}}, 3, \{2\}^{a_{2n}}) \overset{?}{\in} \Q \pi^{\wt} \, ,
		\end{gather*}
		 Notice the blocks of \( 2 \) all have lengths \( a_i \), except for the initial one; it has length \( a_0 + 1 \) in the first identity and length 0 in the second identity.
	\end{Conj}

	These identities are similar in structure to the cyclic insertion conjecture of \cite{borwein1998combinatorial}, on classical MZV's, and should perhaps be regarded as an analogue.  The cyclic insertion conjecture states
	\begin{Conj}[Conjecture 1, \cite{borwein1998combinatorial}]\label{conj:bbbl}
		For any integers \( a_0, \ldots, a_{2n} \geq 0 \) , we have
		\[
		\sum_{\text{cycle \( a_0, \ldots, a_{2n} \)}} \zeta(\{2\}^{a_0}, 1, \{2\}^{a_1}, 3, \{2\}^{a_2}, \ldots, 1, \{2\}^{a_{2n-1}}, 3, \{2\}^{a_{2n}}) \overset{?}{=} \frac{\pi^{\wt}}{(\wt + 1)!} \, .
		\]
	\end{Conj}

	In \cite{charlton2015zeta}, a symmetrised version of \autoref{conj:bbbl} was proven by the author, up to a rational, using the motivic MZV framework of Brown \cite{brown2012mixed,brown2012decomposition}.  This symmetrised result was generalised to a wider class of MZV's using the \emph{alternating block decomposition} of iterated integrals \cite{charlton2017alternating}, along with conjectural cyclic version, since proven exactly by Hirose and Sato. \medskip

	Zhao's generalised 2-1 formula, Theorem 1.4 in \cite{zhao2016identity}, gives an expression for an arbitrary MZSV as a sum of alternating MZV's, with arguments from a certain indexing set \( \Pi(\vec{s}^{(1)}) \).  As a consequence, Zhao gives a concise proof of \autoref{conj:imatomi}.  The goal of this paper is generalise Zhao's proof of \autoref{conj:bbbl}, by connecting \( \vec{s}^{(1)} \) with the `block decomposition' of the multiple zeta value \( \zeta(\vec{s}) \).  This allows us to give analogues of other MZV cyclic insertion identities in the MZSV case. \medskip

	Before stating the main result, we must first recall the construction of the block decomposition from \cite{charlton2017alternating}.  Any word in \( \{0,1\}^\times \) can be written as a concatenation of some number of `alternating words' \( 0, 1, 01, 10, 010, 101, 0101, 1010, \ldots \).  By deconcatenating \( w \) at a repeated letter, one obtains the (unique) decomposition of \( w \) into the minimal possible number of such words.  Moreover, by assuming \( w \) starts with 0, the lengths of the alternating words uniquely determine \( w \) since the concatenation occurs at a repeated letter.

	\begin{Def}[Block decomposition, \cite{charlton2017alternating}]
		For \( w \in \{0, 1\}^\times \), starting with a 0, write \( w \) as a concatenation of the fewest alternating words \( w_1, \ldots, w_n \), with lengths \( \ell_1, \ldots, \ell_n \) respectively.  The block decomposition of \( w \) is
		\[
			\bl(w) \coloneqq (\ell_1, \ldots, \ell_n) \, .
		\]
	\end{Def}

	\begin{Eg}
		For \( w = 01100101010010101 \), we have
		\[
			\bl(w) = \bl(\underbrace{01}_{2} \mid \underbrace{10}_{2} \mid \underbrace{0101010}_{7} \mid \underbrace{010101}_{6}) = (2, 2, 7, 6) \, .
		\]
		From the lengths \( (2, 2, 7, 6) \) we recover a unique word starting with 0, by writing
		\[
			\bl^{-1}(2, 2, 7, 6) = 0\underbracket{1 \mid 1}_{\mathclap{\text{repeat} \hspace{1em}}}\underbracket{0 \mid 0}_{\mathclap{\hspace{1em} \text{repeat}}}10101\underbracket{0 \mid 0}_{\text{repeat}}10101 = w
		\]
	\end{Eg}

	We can now state the main result of this paper.

	\begin{Thm}\label{thm:mzsvcyc}
		For integers \( \ell_i > 1 \), the following identity on MZSV's holds
		\[
		\sum_{\sigma \in S_n} \zeta^\star(\bl^{-1}( 2 \circ \ell_{\sigma(1)}, \ldots, \ell_{\sigma(n)}) ) = \sum_{\vec{r} \in \Part_{\text{odd}}(n)} 2^{\# r} \prod_i (\# r_i - 1)! \prod_i \widetilde{\zeta}\bigg(\sum\limits_{j \in r_i} \ell_j \bigg) \, .
		\]
		Here \( \Part_{\text{odd}}(n) \) is the set of partitions of \( \Set{1, \ldots, n} \) into odd-sized parts.  Moreover
		\[
		\zeta^\star(0 \, 10^{k_1-1} \cdots 10^{k_d-1} \, 1) \coloneqq \zeta^\star(k_1,\ldots,k_d) \, ,
		\]
		as in the iterated integral representation of an MZV, although this also includes the \emph{bounds} of integration.
		Finally, \( \widetilde{\zeta} \) and \( \circ \) are defined by
		\[
		\widetilde{\zeta}(n) = \begin{cases}
		\zeta(n) & \text{\( n \) odd} \\
		\frac{1}{2} \zeta^\star(\{2\}^{n/2}) & \text{\( n \) even,}
		\end{cases}
		\quad \text{ and } \quad
		\circ = \begin{cases}
		, & n \equiv \sum \ell_i \pmod*{2} \\
		+ & n \not\equiv \sum \ell_i \pmod*{2} \, .
		\end{cases}
		\]

		In particular, the sum is always a polynomial in Riemann zeta values, since \( \zeta^\star(\{2\}^n) \in \Q \pi^{2n} \).
	\end{Thm}

	In the case all \( \ell_i \) even, we recover the identities in \autoref{conj:imatomi}, and can give explicit terms in the right hand side in various cases.

	\begin{Eg}
		If \( (\ell_1, \ell_2, \ell_3) = (2a + 2, 2a + 2, 2a + 2) \), we are in the case \( n = 3 \), and \( \circ = \text{``\,+\,''} \).  Then
		\[
			\zeta^\star(\bl^{-1}(2 + 2a + 2, 2b + 2, 2c + 2)) = \zeta^\star(\{2\}^{a + 1}, 1, \{2\}^{b}, 3, \{2\}^{c}) \, .
		\]
		The theorem tell us that we need to sum over \( \vec{r} \in \Part_{\text{odd}}(3) = \Set{ 1\mid2\mid3 \, , 123 } \), and so we obtain
		\begin{align*}
		& {}{} = 2^3(1-1)!^3 \cdot \tfrac{1}{2} \zeta^\star(\{2\}^{a + 1}) \cdot \tfrac{1}{2} \zeta^\star(\{2\}^{b+1}) \cdot \tfrac{1}{2} \zeta^\star(\{2\}^{c+1}) \tag*{\( \leftsquigarrow \vec{r} = \{ 1 \mid 2 \mid 3 \} \)} + {}  \\[1ex]
		& \quad\quad {} + 2^1 (3-1)! \cdot \tfrac{1}{2} \zeta^\star(\{2\}^{a + b + c + 3}) \tag*{\( \leftsquigarrow \vec{r} = \{123\} \)}
		\end{align*}
		which simplifies to
		\[
		= \zeta^\star(\{2\}^{a+1})\zeta^\star(\{2\}^{b+1})\zeta^\star(\{2\}^{c+1}) + 2\zeta^\star(\{2\}^{a+b+c+3}) \, .
		\]

		This gives the identity
		\begin{align*}
			&\sum_{\text{permute \( a, b, c \)}} \zeta^\star(\{2\}^{a + 1}, 1, \{2\}^{b}, 3, \{2\}^{c}) \\
			& \quad = \zeta^\star(\{2\}^{a+1})\zeta^\star(\{2\}^{b+1})\zeta^\star(\{2\}^{c+1}) + 2\zeta^\star(\{2\}^{a+b+c+3}) \in \Q \pi^\wt \, ,
		\end{align*}
		as in case 1 of \autoref{conj:imatomi}.  \medskip

		If \( (\ell_1,\ldots,\ell_4) = (2a + 2, 2b + 2, 2c + 2, 2d+2) \), we are in the case \( n = 4 \), and \( \circ = \text{``\,,\,''} \).  Then
		\[
			\zeta^\star(\bl^{-1}(2, 2a + 2, 2b + 2, 2c + 2, 2d+2)) = \zeta^\star(1,\{2\}^{a}, 3, \{2\}^{b}, 1, \{2\}^{c}, 3, \{2\}^d) \, .
		\]
		We sum over \( \vec{r} \in \Part_{\text{odd}}(4) = \Set{1\mid234, 2 \mid 134, 3 \mid 124, 4 \mid 234, 1 \mid 2 \mid 3 \mid 4} \), and obtain
		\begin{align*}
			&\sum_{\text{permute \( a, b, c, d \)}} \zeta^\star(1,\{2\}^{a}, 3, \{2\}^{b}, 1, \{2\}^{c}, 3, \{2\}^d) \\
			&\quad = 2 \big( \zeta(\{2\}^{a + b + c + 3}) \zeta^\star(\{2\}^{d+1}) + \zeta(\{2\}^{a + b + d + 3}) \zeta^\star(\{2\}^{c+1}) + {} \\
			& \quad\quad {} + \zeta(\{2\}^{a + c + d + 3}) \zeta^\star(\{2\}^{b+1}) + \zeta(\{2\}^{b + c + d + 3}) \zeta^\star(\{2\}^{a+1}) \big) + {}  \\
			& \quad\quad {} + \zeta^\star(\{2\}^{a+1}) \zeta^\star(\{2\}^{b+1}) \zeta^\star(\{2\}^{c+1}) \zeta^\star(\{2\}^{d+1}) \in \Q \pi^\wt \, ,
			\end{align*}
			as in case 2 of \autoref{conj:imatomi}.  \medskip
	\end{Eg}

	\begin{Eg}[Hoffman's identity]
		For integers \( a, b, c \geq 0 \), Hoffman's identity on MZV's (generalised and proven up to \( \mathbb{Q} \) in \cite{charlton2017alternating}, and proven exactly in \cite{hirose2017hoffman}) states
		\begin{align*}
			&\zeta(\{2\}^a, 3, \{2\}^b, 3, \{2\}^c) - \zeta(\{2\}^b, 3, \{2\}^c, 1, 2, \{2\}^a) + \zeta(\{2\}^c, 1, 2, \{2\}^a, 1, 2, \{2\}^b) \\
			& \quad = -\zeta(\{2\}^{a+b+c+3}) \, .
		\end{align*}
		It arises from the block decomposition \( (\ell_1,\ell_2,\ell_3) = (2a + 3, 2b + 3, 2c + 2) \) of the first MZV above. \medskip

		We can apply \autoref{thm:mzsvcyc} to \( (\ell_1,\ell_2,\ell_3) = (2a + 3, 2b + 3, 2c + 2) \) to obtain an analogue on MZSV's.  We are in the case \( \circ = \text{``\,+\,''} \) , and we obtain the following combination of MZSV's.
	\begin{align*}
		& \zeta^\star(\{2\}^{a+1}, 3, \{2\}^b, 3, \{2\}^c) && {} + \zeta^\star(\{2\}^{b+1}, 3, \{2\}^a, 3, \{2\}^c) + {} \\
		&{} + \zeta^\star(\{2\}^{b + 1},3,\{2\}^c,1,2, \{2\}^{a}) && {} + \zeta^\star(\{2\}^{a + 1},3,\{2\}^c,1,2, \{2\}^{b}) + {} \\
		&{} + \zeta^\star(\{2\}^{c + 1}, 1,2, \{2\}^a, 1, 2, \{2\}^b) \hspace{-1em}  && {} + \zeta^\star(\{2\}^{c + 1}, 1,2, \{2\}^b, 1, 2, \{2\}^a)
	\end{align*}
		The theorem tell us that we need to sum over \( \vec{r} \in \Part_{\text{odd}}(3) = \Set{ 1\mid2\mid3 \, , 123 } \), and so we obtain
		\begin{align*}
			& {}{} = 2^3(1-1)!^3 \zeta(2a + 3)\zeta(2b + 3) \cdot \tfrac{1}{2}\zeta^\star(\{2\}^{c+1}) \tag*{\( \leftsquigarrow \vec{r} = \{ 1 \mid 2 \mid 3 \} \)} + {}  \\[1ex]
			& \quad\quad {} + 2^1 (3-1)! \cdot \tfrac{1}{2} \zeta^\star(\{2\}^{a + b + c + 4}) \tag*{\( \leftsquigarrow \vec{r} = \{123\} \)}
		\end{align*}
		which simplifies to
		\[
		= 4 \zeta(2a + 3)\zeta(2b + 3)\zeta^\star(\{2\}^{c+1}) + 2 \zeta^\star(\{2\}^{a + b + c + 4})
		\]
	\end{Eg}

	Similar identities can be given for a wide range of initial block lengths, allowing one to produce identities for many MZSV's with indices 1, 2, 3.  For example
	\begin{Eg}
		Starting with \( \zeta^\star(1, 3, 3, \{2\}^m) \), one reads off the block decomposition
		\[
		\bl^{-1}(0 ; 1 \, 100 \, 100 \, (10)^m ; 1) = (2, 2, 3, 2m+2) \, .
		\]
		By taking \( (\ell_1, \ell_2, \ell_3) = (2, 3, 2m+2) \), we are in the case \( \circ = \text{``\,,\,''} \) as \( \sum \ell_i \equiv 3 \pmod*{2} \).  The theorem gives us the following identity containing \( \zeta^\star(1,3,3,\{2\}^m) \):
		\begin{align*}
			&\zeta^\star(1, 3, 3, \{2\}^m) + \zeta^\star(1, 2, 1, \{2\}^m, 3) + \zeta^\star(1, \{2\}^m,3, 3) + \\
			&\zeta^\star(1, 2, 1, 3, \{2\}^m) + \zeta^\star(1, 3, \{2\}^m, 1, 2) + \zeta^\star(1, \{2\}^m, 3, 1, 2) \\
			& \quad = 2 \zeta(2) \zeta(3) \zeta^\star(\{2\}^{m+1}) + 4 \zeta(2m + 7) \, .
		\end{align*}
	\end{Eg}

	\paragraph{Acknowledgements}  This work was completed during the trimester program ``Periods in Number Theory, Algebraic Geometry and Physics'' at the Hausdorff Institute for Mathematics, concurrent with the author's stay at the Max Planck Institute for Mathematics.  It was motivated by observations made during the author's stay at the MZV research centre, Kyushu university.  I am grateful all three institutes for their hospitality and excellent working conditions.

	I am also grateful to Nobuo Sato for helpful discussions during the trimester, and for directing me to Zhao's generalised 2-1 identity which plays a key role in the proof.

	\section{Proof of the theorem}

	Much of the proof relies on Zhao's generalisation of the 2-1 formula, proven in \cite{zhao2016identity}.  The key step is to relate the block decomposition to Zhao's \( \vec{s}^{(1)} \), with the following lemma.

	\paragraph{Warning:} Since I use the opposite contention for MZV's, the version of \( \vec{s}^{(1)} \) used here is the reverse of the one obtained from Zhao's definition.  These changes are incorporated into the proofs.

	\begin{Lem}
		For \( \circ = \text{``\,,\,'' or ``\,+\,''} \), let \( \vec{L} = (2 \circ \ell_1, \ldots, \ell_n) \) be a block decomposition with all \( \ell_i > 1 \), and corresponding MZV arguments \( \vec{s} = (s_1, \ldots, s_{n'}) \).  Then Zhao's \( \vec{s}^{(1)} \) associated to \( \vec{s} \) is given by
		\[
			\vec{s}^{(1)} = (\widetilde{\ell_1}, \ldots, \widetilde{\ell_i}) \,
		\]
		where
		\[
			\widetilde{\ell_i} = \begin{cases}
				\ell_i & \text{if \( \ell_i \) odd} \\
				\overline{\ell_i} & \text{if \( \ell_i \) even}
			\end{cases}
		\]

		\begin{proof}
			Since block notation, and arguments strings are in \( 1 \leftrightarrow 1 \) correspondence (and block decompositions with \( \ell_i > 1 \) correspond to \( s_i \in \Set{1,2,3} \), with no \( 1, 1 \)), we can build up the block decomposition, and \( s^{(i)} \) term-by-term, and see they match as claimed above.

			Suppose that both \( s^{(2)} = (\widetilde{t_1},\ldots,\widetilde{t_k}) \) and \( \mathbf{L} = (t_1,\ldots,t_k) \) have been constructed for \( \zeta(s_1,\ldots,s_n) \), and that they agree as claimed.  Then for \( \zeta(s_1,\ldots,s_n,s_{n+1}) \) we obtain by Zhao's inductive definition
			\[
				s^{(1)} = \begin{cases}
					(\widetilde{t_1},\ldots,\widetilde{t_k}) \cdot (1) = (\widetilde{t_1},\ldots,\widetilde{t_k},1) & s_{n+1} = 1 \\
					(\widetilde{t_1},\ldots,\widetilde{t_k}) \oplus (2) = (\widetilde{t_1},\ldots,\widetilde{t_{k-1}},\widetilde{t_k} \oplus 2) & s_{n+1} = 2 \\
					(\widetilde{t_1},\ldots,\widetilde{t_k}) \oplus (\overline{1},\overline{2}) = (\widetilde{t_1},\ldots,\widetilde{t_{k-1}}, \widetilde{t_k} \oplus \overline{1}, \overline{2}) & s_{n+1} = 3
				\end{cases}
			\]
			If \( t_k \) is even, then \( \widetilde{t_k} = \overline{t_k} \), so we have \( \widetilde{t_k} \oplus 2 = \overline{t_k} \oplus 2 = \overline{t_k + 2} = \widetilde{t_k + 2} \), and \( \widetilde{t_k} \oplus 1 = \overline{t_k} \oplus \overline{1} = t_k + 1 = \widetilde{t_k + 1} \).  If \( t_k \) is odd, then \( \widetilde{t_k} = t_k \), so we have \( \widetilde{t_k} \oplus 2 = t_k \oplus 2 = t_k + 2 = \widetilde{t_k + 2} \), and \( \widetilde{t_k} \oplus 1 = t_k \oplus \overline{1} = \overline{t_k + 1} = \widetilde{t_k + 1} \).

			So
			\[
			s^{(1)} = \begin{cases}
			 (\widetilde{t_1},\ldots,\widetilde{t_k},1) & s_{n+1} = 1 \\
			(\widetilde{t_1},\ldots,\widetilde{t_{k-1}},\widetilde{t_k+2} & s_{n+1} = 2 \\
			(\widetilde{t_1},\ldots,\widetilde{t_{k-1}}, \widetilde{t_k + 1}, \widetilde{2}) & s_{n+1} = 3
			\end{cases}
			\]

			Now consider the corresponding construction on \( \mathbf{L} \).  If \( s_n = 1 \), then the integral word is
			\[
				\cdots 0 \, 1 \, ; \, 1 = \cdots01 \mid 1 \quad \leftrightarrow \quad \mathbf{L} = (\ldots, \ell_{n-1}, 1) \, .
			\]
			We cannot have \( s_{n+1} = 1 \), so we have either
			\begin{align*}
				s_{n+1} = 2 \quad &\rightsquigarrow\quad  \cdots01 \, \underline{10} \, ; \, 1 = \cdots 01 \mid 101 \\
				& \Rightarrow \quad \mathbf{L} = (\ldots, \ell_{n-1}, 3) = (\ldots, \ell_{n-1}, \ell_n + 2) \, .
			\end{align*}
			Or we have
			\begin{align*}
				s_{n+1} = 3 \quad &\rightsquigarrow \quad \cdots 01 \, \underline{100} \, ; \, 1 = \cdots 01 \mid 10 \mid 01 \\
				&\Rightarrow \quad \mathbf{L} = (\ldots, \ell_{n-1}, 2, 2) = (\ldots, \ell_{n-1}, \ell_n + 1, 2) \, .
			\end{align*}

			If \( s_n \geq 2 \), then the integral word is
			\[
				\ldots 10 ; 1 = \ldots 101 \quad \leftrightarrow \quad \mathbf{L} = (\ldots, \ell_n)
			\]
			We may have \( s_{n+1} = 1, 2, 3 \), so we obtain
			\begin{align*}
				&&& s_{n+1} = 1 && \rightsquigarrow \quad \ldots 10 \, \underline{1} \, ; \, 1 = 101 \mid 1 &&\quad \Rightarrow \quad \mathbf{L} = (\ldots, \ell_n, 1) &\\
				&&& s_{n+1} = 2 && \rightsquigarrow\quad  \ldots 10 \, \underline{10} \, ; \, 1 = 10101 &&\quad \Rightarrow \quad \mathbf{L} = (\ldots, \ell_n + 2) & \\
				&&& s_{n+1} = 3 && \rightsquigarrow\quad  \ldots 10 \, \underline{100} \, ; \, 1 = 1010 \mid 01 &&\quad \Rightarrow \quad \mathbf{L} = (\ldots, \ell_n + 1, 2) \, . &
			\end{align*}

			We see that after inserting the new argument \( s_{n+1} \), we still have matching between \( s^{(1)} \) and \( \mathbf{L} \).

			To complete the proof, we must check the base case holds.  If we start with \( \zeta(1) \) then
			\[
				s^{(1)} = (1) \quad \text{ and } \quad \zeta(1) = I(0;1;1) \quad \Rightarrow \quad \mathbf{L} = (2,1)
			\]

			If we start with \( \zeta(s_1>1) \), then
			\begin{align*}
				s^{(1)} = (\{1\}^{s_1 - 2}, \widetilde{2}) \quad &\text{and} \quad \zeta(s_1) = I(0;10^{s_1-1};1) = I(010 \mid (0 \mid)^{s_1-3} \ldots \mid 01) \\
				&\Rightarrow \quad \mathbf{L} = (2+1, \{1\}^{s_1-3},2) \, .
			\end{align*}
			In both cases, \( \mathbf{L} = (2 \circ \ell_1,\ldots,\ell_n) \) matches with \( s^{(1)} \) as claimed, and the lemma is proven.
		\end{proof}
	\end{Lem}

	\begin{Lem}[Zhao]\label{lem:zhao}
		For any arguments \( \vec{s} \), we have
		\[
			\zeta^\star(\vec{s}) = \epsilon(\vec{s}) \sum_{\vec{p} \in \Pi(\vec{s}^{(1)})} 2^{\#\vec{p}} \zeta(\vec{p}) \, ,
		\]
		where \( \Pi(s_1,\ldots,s_\ell) \) is the set of all indices of the form \( (s_1 \circ \cdots \circ s_\ell) \), where \( \circ \) is either ``\,,\,'' or ``\,$\oplus$\,'', and \( \epsilon(\vec{s}) = 1 \) if \( s_1 = 1 \), \( \epsilon(\vec{s}) = -1 \) if \( s_1 \geq 2 \).

		\begin{proof}
			Apply Theorem 1.4 of Zhao \cite{zhao2016identity}, and pass to the limit \( n \to \infty \) using Lemma 4.5 of Zhao \cite{zhao2016identity}.
		\end{proof}
	\end{Lem}

	\begin{Prop}
		Let \( \vec{s} = (s_1, \ldots, s_{n}) \) be given, and assume \( S_{n} \) acts on the indices \( 1, \ldots, n \) in the standard way.  Then
		\[
			\sum_{\sigma \in S_{n}} \sum_{\vec{p} \in \Pi(\sigma \cdot \vec{s})} 2^{\# \vec{p}} \zeta(\vec{p}) = \sum_{\vec{q} \in \Part(n)} 2^{\# \vec{q}} \prod_i (\# q_i)! \zeta_\sym\bigg(\bigoplus\limits_{j \in q_1} s_j, \ldots, \bigoplus\limits_{j \in q_t} s_j \bigg) \, ,
		\]
0		where \( \zeta_\sym(a_1,\ldots, a_t) \coloneqq \sum_{\tau \in S_td} \zeta(a_{\tau(1)}, \ldots, a_{\tau(t)}) \) symmetrises the arguments of the MZV. \medskip

		That is, one can move the \( S_{n} \) action from \( \mathbf{L} \) to the arguments of the zeta, at the expense of some coefficients.

		\begin{proof}
			Firstly, observe that \( p \in \Pi(\sigma \cdot \vec{s}) \) means \( p = s_{\sigma(1)} \circ \cdots \circ s_{\sigma(n)} \), where each \( \circ \) is \( , \) or \( \oplus \).  But this is equivalent to \( p = \sigma \cdot (s_1 \circ \cdots \circ s_n) \), under the induced \( S_n \) action, and \( (s_1 \circ \cdots \circ s_n) \in \Pi(s_1,\ldots,s_n) \).  So we can write
			\[
				= \sum_{\sigma \in S_n} \sum_{p \in \Pi(\vec{s})} 2^{\# \vec{p}} \zeta(\sigma \cdot  \vec{p}) \, .
			\]
			Warning: \( \sigma \) acts on the elements \( s_i \) inside \( \vec{p} = (s_1 \circ \cdots \circ s_n) \), and not on the comma separated blocks.  So \( \zeta(\sigma \cdot \vec{p}) \) is not simply \( \zeta_\sym(\vec{p}) \).  We need to do further manipulation to obtain the desired form. \medskip

			An element \( \vec{p}\in \Pi(\vec{s}) \) is of the form \( s_1 \circ \cdots \circ s_n \) for some choices \( \circ = \mathop{,} \text{ or } \oplus \).  That is
			\begin{align*}
				\vec{p} &= ( s_1 \oplus \cdots \oplus s_{i_1} , s_{i_1+1} \oplus \cdots \oplus s_{i_2} , \ldots,
				s_{i_{t-1}+1} \oplus \cdots \oplus s_{\underbrace{i_t}_{=n}} ) \text{ , and }\\
				\sigma \cdot \vec{p} &= ( s_{\sigma(1)} \oplus \cdots \oplus s_{\sigma(i_1)} , s_{\sigma(i_1+1)} \oplus \cdots \oplus s_{\sigma(i_2)} , \ldots,
				s_{\sigma(i_{t-1}+1)} \oplus \cdots \oplus s_{\underbrace{\sigma(i_t)}_{i_t=n}} ) \, .
			\end{align*}

			We can therefore define a surjective map
			\begin{align*}
				\phi \colon (\Pi(\vec{p}) , S_n) & \to \Part^\ast(n) \\
				(\vec{p}, \sigma) & \mapsto [\Set{\sigma(1),\ldots,\sigma(i_1)}, \Set{\sigma(i_1 + 1), \ldots, \sigma(i_2)}, \ldots, \Set{\sigma(i_{t-1}+1),\ldots,\sigma(i_t)}]
			\end{align*}
			where \( i_1,\ldots,i_t \) are given by the expression for \( \sigma \cdot \vec{p} \) above.  Here \( \Part^\ast(n) \) is the set-partitions of \( \{ 1, \ldots, n \} \), where the order of the elements of the parts is not important, but the order of the parts themselves is.  That is \( [\Set{a,b}, \Set{c}] = [\Set{b,a},\Set{c}] \), but these are different from \( [\Set{c}, \Set{a,b}] \).

			If \( q = \phi(\vec{p}, \sigma) = [q_1,\ldots,q_t] \), then \( \zeta(\sigma \cdot \vec{p}) = \zeta(\bigoplus_{j \in q_1} s_j, \ldots, \bigoplus_{j \in q_t} s_j) \), and \( \# \vec{q} = \# \vec{p} \).  Notice that \( \# \phi^{-1}(\vec{q}) = \#q_1! \cdots \#q_t! \), since any permutation which respects the parts of \( \vec{p}\) maps to the same \( \vec{q} \).  So we can write that the desired sum is
			\[
				= \sum_{\vec{q} \in \Part^\ast(n)} 2^{\# q} \left( \prod_i \#q_i! \right) \zeta\bigg(\bigoplus_{j \in q_1} s_j, \ldots, \bigoplus_{j \in q_t} s_j\bigg) \, .
			\]

			Finally, we have a surjective map
			\begin{align*}
				\psi \colon \Part^\ast(n) &\to \Part(n) \\
					[q_1,\ldots,q_t] &\mapsto \Set{q_1,\ldots,q_t} \, ,
			\end{align*}
			with \( \psi^{-1}(\Set{q_1,\ldots,q_t})  = \Set{ [q_{\sigma(1)},\ldots,q_{\sigma(t)}] | \sigma \in S_t } \).

			So the sum can be written
			\[
				= \sum_{\vec{q} \in \Part(n)} 2^{\# q} \#q_1! \cdots \#q_t! \underbrace{\sum_{\sigma \in S_t} \zeta\bigg(\bigoplus_{j \in q_{\sigma(1)}} s_j, \ldots, \bigoplus_{j \in q_{\sigma(t)}} s_j\bigg)}_{\eqqcolon \zeta_\sym} \, .
			\]
			This completes the proof.
		\end{proof}
		\end{Prop}

		Using the symmetric sum formula \cite{hoffman1992multiple} (or rather Zhao's generalisation to alternating MZV's in Lemma 5.1 of \cite{zhao2016identity}), one can evaluate the RHS above.  We obtain the following

		\begin{Prop}\label{prop:oddpartsum}
			The following evaluation holds
			\[
				\sum_{\vec{q} \in \Part(n)} 2^{\# q} \bigg( \prod_{i=1}^{\#\vec{q}} \#q_i! \bigg) \; \zeta_\sym\bigg(\bigoplus_{j \in q_{1}} s_j, \ldots, \bigoplus_{j \in q_{t}} s_j\bigg) = \sum_{\vec{r} \in \Part_{\text{odd}}(n)} 2^{\# \vec{r}} \prod_{i=1}^{\#\vec{r}} (\# r_i - 1)! \prod_{i=1}^{\#\vec{r}} \zeta\bigg(\bigoplus_{j \in r_i} s_j\bigg)
			\]
		\end{Prop}

		\begin{proof}
			Firstly, we must apply the symmetric sum theorem to evaluate the left hand side.  It gives
			\[
				= \sum_{\vec{q}\in \Part(n)} 2^{\#\vec{q}} \bigg( \prod_{i=1}^{\#\vec{q}} \# q_i! \bigg) \sum_{\vec{t} \in \Part(\#\vec{q})} (-1)^{\#\vec{q} - \#\vec{t}} \bigg( \prod_{j = 1}^{\#\vec{t}} (\# t_j - 1)! \bigg) \prod_{j=1}^{\#\vec{t}} \zeta\bigg( \bigoplus_{\alpha \in t_j} \bigoplus_{\beta \in q_\alpha} s_\beta \bigg) \, .
			\]
			As the parts \( q_\alpha \) are disjoint, the \( \zeta \) argument
			\[
				\bigoplus_{\alpha \in t_j} \bigoplus_{\beta \in q_\alpha} s_\beta
			\]
			can be written as
			\[
				\bigoplus_{\alpha \in r_j} s_\alpha \, .
			\]
			for some partition \( \vec{j} \in \Part(n) \).  This partition is obtained from \( (\vec{q}, \vec{t}) \) by `flattening' in the following sense
			\begin{align*}
				 f(\vec{q},\vec{t}) = & \vec{r} = \{ r_1,\ldots,r_{\#\vec{t}} \} \text{, where} \\
				& r_i = \bigcup_{j \in t_i} q_j \, .
			\end{align*}
			For example
			\begin{align*}
				& (\vec{q} = \Set{ q_1 = \Set{1, 2, 4} , q_2 = \Set{3, 5}, q_3 = \Set{6,8}, q_4 = \Set{7} }, \vec{t} = \Set{ \Set{1,3,4}, \Set{2} }) \\
				& \mapsto \Set{ q_1 \cup q_3 \cup q_4 , q_2 } = \Set{ \Set{1,2,4,6,7,8}, \Set{3,5} } \, .
			\end{align*}

			So we may formally write the sum as
			\[
				\sum_{\vec{r} \in \Part(n)} \sum_{(\vec{q},\vec{t}) \in f^{-1}(\vec{r})} 2^{\#\vec{q}} (-1)^{\#\vec{q} - \#\vec{t}} \bigg( \prod_{i=1}^{\#\vec{q}} \# q_i! \bigg) \bigg( \prod_{j = 1}^{\#\vec{t}} (\# t_j - 1)! \bigg) \prod_{k=1}^{\#\vec{r}} \zeta\bigg( \bigoplus_{\alpha \in r_k} s_\alpha \bigg) \, .
			\]
			We thus need to evaluate the coefficient
			\[
				c_\vec{r} \coloneqq \sum_{(\vec{q},\vec{t}) \in f^{-1}(\vec{r})} 2^{\#\vec{q}} (-1)^{\#\vec{q} - \#\vec{t}} \bigg( \prod_{i=1}^{\#\vec{q}} \# q_i! \bigg) \bigg( \prod_{j = 1}^{\#\vec{t}} (\# t_j - 1)! \bigg) \, .
			\] We want to show two things: firstly that if the partition \( \vec{r}\) has any even size parts, then the coefficient is 0.  Secondly if the partition only has odd size parts, the coefficient is as indicated in the statement of the proposition.

			Firstly, we can describe \( f^{-1}(\vec{r}) \) more explicitly, as follows.  The elements \( \vec{q} \) which flatten to \( \vec{r} \) are obtained as \( \vec{q} = \coprod_{i=j}^{\#\vec{r}} \vec{T}_j \), where \( \vec{T}_j \) is any partition of \( r_j \).  This choice of partitions \( \vec{T}_1,\ldots,\vec{T}_{\#\vec{r}} \) determines \( \vec{t} \), since \( r_j = \bigcup \vec{T}_j \).  For example
			\begin{align*}
				& \vec{r} = \Set{ \Set{1, 3, 4, 5, 7} , \Set{ 2, 6, 8} } \\
				& \xrightarrow{\text{\( f^{-1} \) contains}} \{ \underbrace{\Set{ 1,4 \mid 3 \mid 5,7 }}_{\vec{T}_1} , \underbrace{\Set{ 2,6 \mid 8 }}_{\vec{T}_2} \} = \Set{\Set{1,4}, \Set{2,6}, \Set{3}, \Set{5,7}, \Set{8}}
			\end{align*}
			and \( \vec{t} = \Set{\Set{1,3,4}, \Set{2,5} } \).

			Under this construction we have
			\begin{gather*}
				\#\vec{t} = \#\vec{r} \, , \quad \quad \# t_i = \# \vec{T_i} \, , \quad \quad \#\vec{q} = \sum_{i=1}^{\#\vec{r}} \# \vec{T}_i \\
				q_j = T_{k,l} \text{ , some \( k, l \), so that} \\
				\prod_{j=1}^{\#\vec{q}} \#q_j! = \prod_{k=1}^{\#\vec{r}} \prod_{l=1}^{\#\vec{T}_k} \#T_{k,l}!
			\end{gather*}
			Thus
			\[
					c_\vec{r} = (-1)^{\#\vec{r}}\sum_{T_1 \in \Part(\#r_1)} \cdots \sum_{T_{\#\vec{r}} \in \Part(\#r_{\#\vec{r}})} (-2)^{\sum_{i} \# T_i} \cdot \prod_{k=1}^{\#\vec{r}} \prod_{l=1}^{\#\vec{T}_k} \#T_{k,l}! \cdot \prod_{j = 1}^{\#\vec{r}} (\# T_j - 1)! \, .
			\]
			This sum can now be factored into a product of the following form
			\[
				c_\vec{r}= (-1)^{\#\vec{r}} \prod_{i=1}^{\#\vec{r}} g(\#r_i) \, ,
			\]
			where
			\[
				g(i) \coloneqq \sum_{\vec{w} \in \Part(i)} (-2)^{\#\vec{w}} \prod_{l=1}^{\#\vec{w}} \#w_l! \cdot (\#\vec{w} - 1)! \, .
			\]

			I claim that \( g \) can be evaluated as follows
			\[
				g(i) = \begin{cases} -2 (i-1)! & \text{ \( i \) odd } \\
					0 & \text{ \( i \) even.}
					\end{cases}
			\]
			If this claim does hold, then
			\begin{align*}
				c_\vec{r} &= \begin{cases}
					0 & \text{some \( \#r_i \) even} \\
					(-1)^{\#\vec{r}} \prod_{i=1}^{\#\vec{r}} (-2) (\#r_i-1)! & \text{all \( \#r_i \) odd}
				\end{cases} \\
				&= \begin{cases}
					0 & \text{some \( \#r_i \) even} \\
					2^{\#\vec{r}} \prod_{i=1}^{\#\vec{r}} (\#r_i-1)! & \text{all \( \#r_i \) odd.}
				\end{cases}
			\end{align*}
			So the proposition will follow.
		\end{proof}

		For the proof to be complete, we need to show the following claim.

		\begin{Clm}
			Let
			\[
				g(n) \coloneqq \sum_{\vec{w} \in \Part(n)} (-2)^{\#\vec{w}} \prod_{l=1}^{\#\vec{w}} \#w_l! \cdot (\#\vec{w} - 1)! \, ,
			\]
			then
			\[
				g(n) = \begin{cases} -2 (n-1)! & \text{ \( n \) odd } \\
				0 & \text{ \( n \) even.}
				\end{cases}
			\]

			\begin{proof}
				We show the generalised identity
				\[
					\frac{1}{n!} g(n, x) = -\frac{1}{n} + \frac{1}{n} (1 + x)^n \, ,
				\]
				where
				\[
					g(n,x) \coloneqq \sum_{\vec{w} \in \Part(n)} x^{\#\vec{w}} \prod_{l=1}^{\#\vec{w}} \#w_l! \cdot (\#\vec{w} - 1)!
				\]
				Hence for \( x = -2 \), we obtain
				\[
					g(n,-2) = -(n-1)! + (n-1)! (-1)^n = \begin{cases}
					-2(n-1)! & \text{\( n \) odd} \\
						0 & \text{\( n \) even,}
					\end{cases}
				\]
				as claimed. \medskip

				To show the generalised identity, we can first show that the derivatives agree.  Then integrating gives
				\[
					\frac{1}{n!} g(n,x) = c + \frac{1}{n} (1 + x)^n \, ,
				\]
				for some constant \( c \).  One sees that \( c = -\frac{1}{n} \) by setting \( x = 0 \); the left hand side is 0 and the right hand side is \( c + \frac{1}{n} \), which proves the claim. \medskip

				To see the derivatives agree, we must show
				\begin{align}
 					& \frac{1}{n!} g'(n,x) = (1 + x)^{n-1} \, \text{ equivalently, } \nonumber \\
					& g'(n,x) = n! (1 + x)^{n-1} \label{eqn:gder} \, .
				\end{align}
				Term-by-term differentiation of \( g(n,x) \) gives
				\[
					g'(n,x) = \sum_{\vec{w} \in \Part(n)} x^{\#\vec{w}-1} \prod_{l=1}^{\#\vec{w}} \#w_l! \cdot \#\vec{w}! \, .
				\]
				The coefficient of \( x^{i-1} \) on the right hand side of \autoref{eqn:gder} is
				\[
					n! \binom{n-1}{i-1} \, .
				\]
				The coefficient of \( x^{i-1} \) on the left hand side is
				\[
					\sum_{\substack{\vec{w} \in \Part(n) \\
						\#\vec{w} = i}} \prod_{l=1}^{\#\vec{w}} w_l! \cdot \underbrace{\#\vec{w}!}_{i!}
				\]
				These two expressions give two different ways to count the number of ordered partitions of \( [1,\ldots,n] \) into \( i \) non-empty ordered parts, and hence are equal.  Here an ordered partition with ordered parts means that \( [[1,2],[3]] \), \( [[2,1],[3]] \), \([[3],[1,2]] \) and \( [[3],[2,1]] \) are all counted as distinct.  For simplicity, refer to such a partition as an ordered/ordered partition. \medskip

				We can form an ordered/ordered partition of \( [1,\ldots,n] \) into \( i \) parts by first taking any permutation of \( [1,\ldots,n] \), then inserting \( i-1 \) bars into any choice of the \( n-1 \) gaps, breaking the \( i \) non-empty parts.
				\begin{align*}
					[1,\ldots,8] & {} \xrightarrow[\phantom{\text{insert bars}}]{\text{permute}} [4\underset{\uparrow}{,} 5, 2, 3\underset{\uparrow}{,} 7, 6\underset{\uparrow}{,} 1, 8] \\
					& {} \xrightarrow{\text{insert bars}} [[4],[5,2,3],[7,6],[1,8]] \, .
				\end{align*}

				There are \( n! \) permutations, and \( \binom{n-1}{i-1} \) ways of choosing \( i - 1 \) positions from the \( n - 1 \) gaps.  This gives the right hand side.

				Alternatively, we can form an ordered/ordered partition of \( [1,\ldots,n] \) by taking a set-partition of \( \Set{1,\ldots,n} \) into \( i \) parts, then reordering the \( i \) parts arbitrarily, as well as arbitrarily reordering the elements of each part.  Every such ordered/ordered partition of \( [1,\ldots,n] \) arises in this way, for some unique \( \vec{w} \), as forgetting about both orderings gives a surjection onto \( \Part(n) \).
				\begin{align*}
					\vec{w} \in \Part(n) & \xrightarrow[\phantom{\text{permute parts}}]{} \Set{\Set{1,8},\Set{2,3,5},\Set{4},\Set{6,7}} \\
					& \xrightarrow[\text{and elements}]{\text{permute parts}} [[4],[5,2,3],[7,6],[1,8]]
				\end{align*}
				Let \( \vec{w} \in \Part(n) \) be a (set-)partition of \( \Set{1,\ldots,n} \) into \( i \) parts, with sizes of each part \( \#w_1, \ldots, \#w_i \) respectively.  Then there are \( i! \prod_{l=1}^{\#\vec{w}} \#w_l! \) such ordered/ordered partitions arising from \( \vec{w} \).  We must sum over all such \( \vec{w} \in \Part(n) \), giving
				\[
					\sum_{\substack{\vec{w} \in \Part(n) \\ \#\vec{w} = i}}  i! \prod_{l=1}^{\#\vec{w}} \#w_1! \, .
				\]
				This is the left hand side. \medskip

				The coefficients of both sides agree, hence we get the required equality of derivatives, and so claim follows.
			\end{proof}
		\end{Clm}

		Finally, we can uses these results to prove the theorem.

		\begin{proof}[Proof of theorem]
			Using \autoref{lem:zhao}, we have that
			\[
				\sum_{\sigma} \zeta^\star(\bl^{-1}(2 \circ \ell_{\sigma(1)}, \ldots, \ell_{\sigma(n)})) = \epsilon(\circ) \sum_{\sigma \in S_n} \sum_{\vec{p} \in \Pi(\sigma \cdot (\widetilde{\ell_1}, \ldots, \widetilde{\ell_n}))} 2^{\# \vec{p}} \zeta(\vec{p}) \, .
			\]
			Notice \( \epsilon(\circ) \) matches with \( \epsilon(s) \), for the following reason.  If \( \circ = \text{``\,,\,''} \), then \( \zeta^\star(\bl^{-1}(2, \ell_1, \ldots)) = I(01 \mid 1 \cdots) = \zeta(1,\ldots) \) and \( \epsilon(\vec{s}) = 1 \) since \( s_1 = 1 \).  Otherwise \( \circ = ``\,+\,'' \), then \( \zeta_\star(\bl^{-1}(2+\ell_1,\ldots)) = I(0101\cdots) = \zeta(2,\ldots) \), and \( \epsilon(\vec{s}) = -1 \) since \( s_1 = 2 \).

			Interchange the summations, write the result as a sum over odd-sized partitions using the \autoref{prop:oddpartsum}
			\[
				 = 	\epsilon(\circ) \sum_{\vec{r} \in \Part_{\text{odd}}(n)} 2^{\# \vec{r}} \prod_i (\# r_i - 1)! \prod_i \zeta\Big(\bigoplus_{j \in r_i} \widetilde{\ell_j}\Big) \, .
			\]

			Since the size of each partition is odd, we can explicitly evaluate \( \bigoplus_{j \in p_i} \widetilde{\ell_j} \), and the resulting \( \zeta \) as follows.

			\paragraph{Case \( \# \{ \ell_i | \text{ even } \} \equiv 0 \pmod*{2} \):}  Then
			\[
				\bigoplus_{j \in p_i} \widetilde{\ell_j} = \sum_{j \in p_i} \ell_j \, ,
			\]
			and this is sum is odd.  This is because the number of bars is additive (i.e. the sign is multiplicative), and there are an even number of bars in total.  So the \( \oplus \) sum agrees with the sum of the undecorated \( \ell_i \).  Moreover, the total is odd since we sum an odd number of odd numbers.

			Overall, this means
			\[
				\zeta\bigg(\bigoplus_{j \in p_i} \widetilde{\ell_j} \bigg) = \zeta\bigg( \sum_{j \in p_i} \ell_j \bigg) \, .
			\]

			\paragraph{Case \( \# \{ \ell_i | \text{ even } \} \equiv 1 \pmod*{2} \):}  Then
			\[
				\bigoplus_{j \in p_i} \widetilde{\ell_j} = \overline{\sum_{j \in p_i} \ell_j} \, ,
			\]
			 and this sum is even.  This is because there are an odd number of bars in total, so one remains after doing the \( \oplus \)-sum.  Consequently the \( \oplus \)-sum agrees with the bar of the undecorated sum.  Moreover, the total is even, since we add an even number of odd numbers.

			 This means
			 \[
			 	\zeta\bigg( \bigoplus_{j \in \pi_i} \widetilde{\ell_j} \bigg) = \zeta\bigg( \overline{\sum_{j \in p_i}} \ell_j \bigg) \, .
			 \]  We can now use Zlobin's evaluation \cite{zlobin2005brief} \( \zeta^\star(\{2\}^n) = -2\zeta(\overline{2n}) \),  (which is also contained in Zhao's \(2\)-\(1\) theorem) to write
			 \[
			 	\zeta\bigg( \overline{\sum_{j \in p_i}} \ell_j \bigg) = -\frac{1}{2} \zeta^\star\bigg(\{2\}^{\text{\normalsize $ \frac{1}{2} \sum\nolimits_{j \in p_i} \ell_j $}} \bigg) \, .
			 \]
			 This is almost our definition of \( \widetilde{\zeta} \).  I claim that the number of \( -1 \) signs between \( \epsilon(\circ) \) and all \( -\zeta^\star(\{2\}^n) \)'s is even. We may discard it to obtain an equivalent formula with our original definition of \( \widetilde{\zeta} \). \medskip

			 Why is the total number of \( -1 \)'s even?

			 \paragraph{Case \( \circ = \text{``\,,\,''} \):} Here \( \epsilon(\mathop{,}) = 1 \), since the MZV's begin \( \zeta^\star(\bl^{-1}(2,\geq2)) = \zeta^\star(01\mid10\cdots) = \zeta^\star(1,\geq2) \).  I claim that the number of `even-sum' parts is even, hence the total number of \( -1 \)'s is even as claimed.  In this case \( n \equiv \sum \ell_i \pmod*{2} \), and we can check \( n \) odd or even separately.

			 Suppose \( n \) odd, then \( \sum \ell_i \) also odd.  By counting the number \( n \) of \( \ell_i \), we see \( \vec{p} \in  \Part_{\text{odd}}(\ell_i) \) has an odd number of parts.  If an odd number of parts have even sum, we would have \( \sum \ell_i \) even, a contradiction.

			 Similarly if \( n \) even, then \( \sum \ell_i \) also even.  By counting the number \( n \) of \( \ell_i \), we see \( \vec{p}\in \Part_{\text{odd}}(\ell_i) \) has an even number of parts.  If an odd number of parts have even sum, we would again have \( \sum \ell_i \) odd, a contradiction. \medskip

			 \paragraph{Case \( \circ = \text{``\,+\,''} \):} Here \( \epsilon(\mathop{+}) = -1 \), since the MZV's begin \( \zeta^\star(\bl^{-1}(2+\geq2, \ldots)) = \zeta^\star(0101 \cdots) = \zeta^\star(2, \ldots) \).  I claim that the number of `even-sum' parts is odd, hence the total number of \( -1 \)'s is even as claimed.  In this case \( n \not\equiv \sum \ell_i \pmod*{2} \), so just check \( n \) odd or even separately.

			 Suppose \( n \) odd, then \( \sum \ell_i \) is even.  By counting the number \( n \) of \( \ell_i \), we see that \( \vec{p}\in \Part_{\text{odd}}(\ell_i) \) has an odd number of parts.  If an even number of parts have even sum, we would obtain \( \sum \ell_i \) odd.

			 Finally \( n \) even, so \( \sum \ell_i \) is odd.  By counting the number \( n \) of \( \ell_i \), we see that \( \vec{p}\in \Part_{\text{odd}}(\ell_i) \) has an even number of parts.  If an even number of parts have even sum, we would obtain \( \sum \ell_i \) even. \medskip

			 In all cases the overall number of \( -1 \)'s is even and we can drop the \( -1 \) from the definition of \( \widetilde{\zeta} \), to obtain the required result.  This completes the proof of the theorem.
		\end{proof}

	\bibliography{cyc_mzsv}
	\bibliographystyle{plainurl}

\end{document}

config: arxiv


TITLE: What does the "interior angle of the loop" mean for a curve that self-intersects?
QUESTION [2 upvotes]: How does one measure that angle? I've looked but couldn't find a definition for it anywhere.

REPLY [1 votes]: The double point $p=X(s_1)=X(s_2)$ divides the curve $X$ into two loops, one of which we call $\Gamma.$ This $\Gamma$ is then a closed curve that is smooth everywhere except at p, where it necessarily must have a sharp corner (since the self-intersection is transversal). The angle $\alpha$ is simply the interior angle at this corner, as illustrated here:

To be more rigorous, $\alpha$ can be defined as the angle between the tangent vectors $X'(s_1)$ and $-X'(s_2).$ Note that if the loop $\Gamma$ was smooth then we would have $\alpha = \pi;$ so the angle defect (or equivalently, the point mass of curvature) at $p$ is $k_p = \pi - \alpha.$ Thus the claim $$\int_\Gamma k = \pi + \alpha$$ can be rewritten as $$k_p + \int_\Gamma k = 2\pi,$$ which is the Theorem of Turning Tangents; i.e. the familiar fact that the total curvature of a simple closed curve is $2\pi.$


TITLE: How to plot frequency vs time graph of an audio signal/recording?
QUESTION [2 upvotes]: I need to plot a frequency variation with time graph for an audio signal that I will record, to observe the variation in frequency with time. I'm completely new to audio signal analysis and have no experience with any related software or matlab. I'd really appreciate some info on the simplest software to use for making such a plot with least the amount of technical knowledge required. Thanks in advance.

REPLY [3 votes]: One of the simplest ways to observe frequency variation with respect to time is called a spectrogram.  This type of plots stacks Fourier spectra together into an image that looks very much like a heat map.
These are built-in to a number of software packages.  You can use the FREE (and reasonably user-friendly!) audio editing program called Audacity to view your signal's spectrogram.  Here is a link to using that functionality:
http://manual.audacityteam.org/man/spectrogram_view.html
When you do start using something like MATLAB for your tasks, you will find spectrogram tools readily available and very easy to use as well.  MATLAB's tool is called just that, 'spectrogram'.


TITLE: Evaluate $\sin^5 x$ using Imaginary part of $e^{\iota x}$
QUESTION [0 upvotes]: A few days ago somebody posted a problem on evaluating $\int sin^5(x) dx$. The answer posted by Jack D'Aurizio Sir used Complex Numbers to represent $\sin(x)$. This lead to me trying to think of a way to evaluate

$$\int \sin^5 x \text{ using } \Im(e^{\iota x}) $$

An example of this method is as follows:
$$\int e^x \cos x \, dx$$
Since $\cos x$ is the real part of $e^{\iota x}$, we know that
$$\int e^x \cos x dx = \Re\int e^x e^{\iota x}\, dx.$$
Is there any way we can evaluate $\int \sin^5 xdx$ using a similar technique?

REPLY [1 votes]: $$2\iota \sin(\theta) = x -\dfrac{1}{x},$$ where $x = e^{\iota \theta}$. Then,
$$\Rightarrow 2^n \iota^n \sin^n(\theta) = \left(x-\dfrac{1}{x}\right)^n $$
and expanding the R.H.S yields
$$2^4\sin^5(θ) = \sin(5θ)-5 \sin(3θ)+10 \sin(θ)$$
For integration this should be easy.


TITLE: Looking for example of function between metric spaces such that...
QUESTION [0 upvotes]: $f: X \to Y$ is continuous and $E \subset Y$
such that
$f^{-1}(\overline{E}) \not\subset \overline{f^{-1}(E)}$
any tips of what to try?

REPLY [2 votes]: HINT: Arrange matters so that $E$ is not in the range of $f$, but $E$ has at least one limit point that is in the range of $f$.