Discrete event data arises in a variety of forms including crimes, health events, neural firings, and social media posts. Frequently each event can be associated with a node in a network, and practitioners aim to learn the relationships between the nodes in the network from the event data. For example, one might observe a sequence of crimes associated with different gangs and seek to learn which crimes are most likely to spark retaliatory violence from rival gangs.

Such problems have attracted widespread interest in recent decades and a variety of point process models have been proposed to model such data. The Hawkes process ([13]) is a popular continuous time framework which has been applied in a number of different contexts (e.g., [34, 33, 8, 31]). In addition, other works works have used a discrete time framework to model time series event data (e.g., [17, 10, 3]).

A central assumption of many of these works is that all the events are observed. However, in many cases we may observe only a subset of the events at random. For example, while point process models have been widely used to model crime incidence ([23, 24, 18]), frequently one only has access to reported crime data. For many crimes the true number of incidents can be substantially higher. The gap between the reported and true crime rates is referred to as “the dark figure of crime” by researchers in Sociology and Criminology who have studied this issue extensively ([6, 15]). Unobserved events also pose a challenge in inference from Electronic Health Record (EHR) data which can be incomplete for a number of different reasons ([30, 29]).

The unobserved events still play a role in the dynamical system governing the time series, making network estimation with incomplete data particularly challenging. In this paper, we contribute to the growing literature on modeling in the presence of unobserved events by proposing a novel method for estimating networks when we only observe a subset of the true events.

1.1 Problem Formulation

Many point process models of time series of discrete events are temporally discretized either because event times were discretized during the data collection process or for computational reasons. In such contexts, the temporal discretization is typically such that either one or zero events are observed in each discrete time block for each network node. With this in mind, we model the true but unobserved observations using a Bernoulli autoregressive process:

Here is a vector indicating whether events occurred in each of the M nodes during time period t. The vector is a constant bias term, and the matrix is the weighted adjacency matrix associated with the network we wish to estimate. We assume that each row a of lies in the ball of radius r, which we denote . We generally consider a case where a is sparse and the magnitude of all its entries are bounded, so that r is a universal constant which is independent of M.

We observe , a corrupted version of (1.1) where only a fraction of

events are observed as follows:

Here denotes the Hadamard product and is a vector where each entry is independently drawn to be one with probability p and zero with probability

Our analysis of (1.1) and (1.2) can be naturally extended to several more complex variants. Instead of assuming each is observed with probability p, we can assume events from each node i are observed with a unique probability . We consider only a first order AR(1) process but our framework can be extended to incorporate more sophisticated types of memory as in [22]. We omit discussion of these extensions in the interest of clarity.

Corrupted or missing data in high-dimensional data sets is a problem which appears in a number of different domains and has attracted widespread interest over the past few decades (see [11] for an application-focused overview). Our focus is on a particular type of corrupted data: partially observed sequences of discrete events. In recent years researchers have started to focus on this problem ([32, 28, 16]). The prior works of which we are aware use a Hawkes process framework and assume knowledge of the time periods when the data is corrupted. In the context of (1.2) this amounts to knowledge of . Our method can operate in a setting where the researcher cannot differentiate between periods when no event actually occurred, and when events potentially occurred but were not recorded. Moreover, because we use a discrete-time framework, we are able to derive sample complexity bounds for the estimation procedure proposed in Section 3. Our theoretical results complement the empirical nature of much of the past work in this area.

This paper is also related a variety of works on regularized estimation in high-dimensional statistics, including [5, 27, 4] and [14]. Many of these works have derived sample complexity guarantees using linear models, and some of these results have been extended to autoregressive generalized linear models ([25, 12]). Another line of research ([20, 1, 21, 19, 25]) has formalized a notion of Restricted Strong Convexity (RSC) which we leverage in Section 4.2. While many loss functions of interest in high-dimensional statistics are not strongly convex, these works have shown that they frequently satisfy an RSC condition which is sufficient to provide statistical and optimization error guarantees. The main technical challenges in our setting lie in establishing results similar to these RSC conditions.

2.1 Missing Data in a High-Dimensional Linear Model

[20] straddles the missing data literature and high-dimensional statistics literature. The

authors consider a missing data linear model

where and one observes pairs and aims to estimate . The authors propose minimizing a loss function of the observed data Z which satisfies the property

for any

denotes the classical Lasso loss function with the unobserved data X and regularization parameter . In other words, the missing data loss function is an unbiased estimator for the full data Lasso loss function we would ideally like to minimize. This idea motivates our construction of a loss function for the observed process (1.2).

Our problem can be viewed as an extension of [20] to autoregressive GLM models without knowledge of W.1 In particular, we cannot distinguish events that were missed () from correctly observed periods with no events (able prove sample complexity bounds as well as optimization bounds which are consistent with the high-dimensional statistical literature in that they scale with rather than the dimension of . We are able to prove analogous bounds for our estimator in Section 4.

2.2 Contributions

This paper makes the following contributions.

• We propose a novel method for estimating a network when only a subset of the true events are observed. In contrast to previous work, our methods do not rely on knowledge of when the data is potentially missing. Our procedure uses Taylor series approximations to an unbiased loss function, and we show that these approximations have controlled bias and lead to accurate and efficient estimators.

• We prove bounds on both the statistical error and optimization error of our proposed estimation method. The results hinge on showing that our loss function satisfies a restricted strong convexity (RSC) condition. Past work on linear inverse problems with corrupted designs also establish RSC conditions, but these conditions do not carry over to the autoregressive GLM setting.

• We demonstrate the effectiveness of our methodology on both simulated data and real crime data.

Given the full data , the negative log-likelihood function is decomposable in the M rows of A. In other words, if

then where is the mth row of A and denote the loss function restricted to a specific row. Throughout the paper we slightly abuse notation and let refer to the entire loss function when A is a matrix, and let refer to the loss function for a specific row when is a row vector. The loss function for the mth row takes the form

where f(x) = log(1 + exp(x)) is the partition function for the Bernoulli GLM.

We do not have access to X and instead we aim to estimate A using the corrupted data . As discussed in Section 2.1, our strategy will be to construct a loss function of Z which is an unbiased estimator for . In other words, we want to find some function such that for any

In contrast to the Gaussian case discussed in Section 2.1, the Bernoulli partition function f(x) = log(1 + exp(x)) is not a polynomial and constructing a function satisfying (3.1) directly is challenging. We adopt a strategy of computing unbiased approximations to truncated Taylor series expansions of and arriving at as a limit of such approximations.

To do this, we first rewrite f using its Taylor series expansion around zero

The constant factor log(2) does not effect estimation in any way so we ignore it for the remainder of our discussion in the interest of simplicity. We let denote the degree q Taylor truncation to . The are binary vectors and we assume each row is sparse, so will not be too far from zero. Thus it is reasonable to hope that for small is a good approximation for whenever . We bound the approximation error in Lemma B.3 in the supplement.

We now consider the problem of constructing a function

Key question we need to address with this approach include (a) can we (approximately) solve this optimization problem efficiently? (b) will the solution to this optimization problem be robust to initialization? (c) will it be a strong estimate of the ground truth?

3.1 Definition of L(2)Z,p

We first derive an unbiased estimator of the degree-two Taylor series expansion

Note that there are straightforward unbiased estimates of the first two terms:

For the third term, , note that

Thus we must estimate the monomials with repeat terms (i = j) differently from the monomials with all distinct terms (). Using Equations (3.3) and (3.4) we can define the degree two unbiased estimator:

3.2 Higher-Order Expansions

The construction of in the previous section suggests a general strategy for constructing satisfying (3.2). Take any monomial

depending on the counts in nodes during time period t. Wherever this monomial appears in , our unbiased loss function will have a term

scaled by where k denotes the number of unique terms in the monomial. For example, in Equation (3.3) each degree two monomial was unique so we scaled everything by . However, in (3.5) some of the degree two monomials had repeated terms and so they were scaled by . In order to formalize these ideas and generalize our estimator to q > 2, we first need to introduce additional notation.

3.2.1 Notation

Let denote the set of all monomials of degree d. We represent an element as a list containing d elements. An element in the list corresponds to the index of a term in the monomial. For an example, the monomial can be represented as the list (1, 1, 3).

For a polynomial function h we let denote the coefficient of the monomial U in h. Finally we define the order of a list to denote the number of unique elements in the list, so |(1, 2)| = 2 whereas |(1, 1)| = 1.

Example Consider the function . We can decompose h as

where with corresponding coefficients and Using this notation we can write

3.2.2 Definition of L(q)Z,p

The degree q likelihood is constructed as follows

Recall that |U| denotes the number of unique terms in the monomial U. In other words, in we adjust by scaling each monomial according the the number of unique terms rather than the number of overall terms. This definition clearly satisfies (3.2). We show in the supplement that if r = 1 and then converges uniformly on to a function we denote as . Extending this loss function on individual rows to an entire matrix, we can define . An additional technical discussion in the supplement guarantees that actually satisfies the desired property in Equation (3.1).

3.3 Proposed Optimization

In practice we can only compute . To estimate we consider the following constrained optimization:

where is an estimate of the missingness parameter p and

In general, is not a convex function. However, we show in Section 4.2 that under certain assumptions all stationary points of (3.6) must lie near . Thus we can approximately solve (3.6) via a simple projected gradient descent algorithm.

In order to apply our algorithm it is necessary to have an estimate of the frequency of missed data. In many cases one may have prior knowledge available. For example, social scientists have attempted to quantify the frequency of unreported crimes ([15, 26]). Moreover, a simulation study in Section 5 suggests our strategy is robust to misspecification of p.

In this section we answer several questions which naturally arise with our proposed estimation procedure in Equation (3.6). In particular, we’d like to know:

• If we can find a solution to Equation (3.6), will it be a good approximation to

• The loss function is not convex. Is it possible to actually find a minimizer to (3.6), or an approximation to it?

In Theorem 4.1 we answer the first question affirmatively. In Theorem 4.2 we show that all stationary points of (3.6) are near one another, so that simple gradient-based methods will give us a good solution.

Throughout this section we assume . All results in the section apply for the loss functions for . In the case we recover the idealized loss function

We use to mean and to mean a = Cb where C is a universal constant. Define

4.1 Statistical Error

The following theorem controls the statistical error of our proposed estimator.

Theorem 4.1 (Accuracy of

where

for with probability at least

The two terms in the upper bound of Theorem 4.2 have a natural interpretation. The first represents the error for the idealized estimator , while the second represents the error due to the Taylor series truncation. Our error scales as which is reasonable in the context of our algorithm because is only well-defined when (see Remark 2 in the supplement). An interesting open question which arises from Theorem 4.1 is whether the process described in (1.1) and (1.2) is unidentifiable for or something specific to our methodology fails for p below this threshold.

The proof of Theorem 4.1 uses ideas from the analysis of high-dimensional GLMs ([12, 22]) as well as ideas from the analysis of missing data in the linear model [20] and Gaussian linear autoregressive processes [14]. The key technical challenge in the proof lies in controlling the gradient of the error term . This is done in Lemmas B.2-B.6 in the supplement.

4.2 Optimization Error

We next focus on the optimization aspects of Equation (3.6). Our loss function is non-convex, so at first glance it may appear to be a poor proxy loss function to optimize. However, a body of research (see [1, 21, 20, 19]) has studied loss functions satisfying a properties known as restricted strong convexity (RSC) and restricted smoothness (RSM). These works have shown that under certain conditions, the optimization of non-convex loss functions may be tractable. The formal definitions of the RSC and RSM conditions we use are as follows.

Definition 1 (Restricted Strong Convexity). Let denote the first order Taylor approximation to a loss function L centered at w. A loss function L satisfies the RSC condition with parameters

for all

Definition 2 (Restricted Smoothness). A loss function L satisfies the RSM condition with parameters

for all

We are able to show these conditions are satisfied for

This in turn gives the following result. As in Theorem 4.1 we assume

Theorem 4.2. Suppose and for at least rows of where C is a universal constant. Let be any stationary point of where

with probability at least

As in Theorem 4.1 the first term in our bound can be interpreted as the error for the idealized estimator while the second term can be thought of as the error due to the Taylor series truncation. The assumption that for at least rows of says that at least a constant fraction of nodes are influenced by other nodes in the network. This assumption allows us to state Theorem 4.2 in terms of s - the support of . In extreme cases where almost all nodes in the network fire independently of the other nodes it is possible for the optimization error to have a slower scaling than s.

The RSC and RSM conditions are closely related to ideas used in our statistical error bounds in Theorem 4.1. Lemma D.2 shows that the conditions are satisfied for on the order of which leads to an overall optimization error bound of the same order. This is a slower convergence rate than in the linear case; whether stronger rates can be obtained in the autoregressive GLM setting is an open question.

The proof of Theorem 4.2 proceeds as follows. We first establish that the RSC/RSM conditions hold for reasonable constants in Lemma D.2. This proofs relies on the technical machinery built up in Lemmas B.2-B.6. We can then combine our RSC/RSM results with Theorem 2 in [1] to conclude that all stationary points of lie in a small neighborhood of with high probability.

In this section we evaluate the proposed method on synthetic data. We generate with s = 50 nonzero entries with locations chosen at random. Each nonzero entry is chosen uniformly in and . We then generate a “true" data set X and an “observed”

Figure 1: MSE vs T for p = .6 (top) and p = .75 (bottom). The blue line uses regularized MLE on full data – i.e. data unavailable in our setup – and represents a kind of oracle estimator. The red line uses incomplete data with (our proposed method). The yellow lines corresponds to minimizing the full data likelihood over the incomplete data – that is, this estimator naively ignores the issue of missing data. Median of 50 trials is shown and error bars denote sample standard deviations.

data set Z according (1.1) and (1.2) with . We perform projected gradient descent with a random initialization and show a median of 50 trials.

Figure 1 shows mean squared error (MSE) vs T for p = .6 (top) and p = .75 (bottom). Our method is shown in red. It uses the loss function on the partially observed data Z. Our method is compared to the loss function using both the full data X (i.e., an oracle estimator with access to the missing data) and the partially observed data Z (i.e. a naive estimator that ignores the possibility of missing data). As expected, with access to the complete data one can get a more accurate estimate of than either method using the partially observed data. However, our method outperforms the full data likelihood when given the partially observed data. In particular, note that the accuracy for the full data likelihood stalls after some time due to the inherent bias in using the corrupted data on the true data likelihood. In contrast our unbiased method continues to converge to the solution, as suggested by the results in Section 4. Finally, observe that for large T there is little variation between trials when using even though each trial was initialized randomly. This agrees Theorem 4.2 which states that all stationary points of lie near one another.

In practical applications one may have strong reason to believe some events are unobserved, but pinning down a precise estimate of the missingness parameter p might be unrealistic. Therefore it is important to see how our algorithm performs as a function of the misspecification . We examine this in Figure 2. We generate data as in the previous section but with p = .7. We then apply our algorithm with the loss function and varying values of

Figure 2 shows that our method is highly robust to small misspecification of the missingness parameter p. Interestingly, underestimating p by more than 10% leads to poor results but our method is still robust to overestimation of over 10%. This suggests there is value in applying our techniques with a conservative estimate of the amount of missed data, even when one has only a rough estimate of the frequency of missed events.

Figure 2: Robustness to missestimation of p using a true value of p = .7. Median of 50 trials is shown and errorbars denote sample standard deviations.

As final experiment we measure how MSE varies as a function of the Taylor series truncation level q. Calculating takes a significant amount of time for high-dimensional problems, so we randomly generate nonzero entries compared to 50 in the previous simulations and run 30 trials. We set

In Figure 3 we show MSE as a function of T for the full data loss function at three different truncation levels: and . Recall that and has no odd degree terms other than 1, so and . We see that the second and fourth degree truncations perform essentially the same as the full data likelihood. We also plot MSE as a function of T for the truncated missing data loss functions . As expected, using the full data gives stronger results than the partially observed data. We again see that the second and fourth order truncations perform nearly the same. The sample standard deviations are also similar - e.g. when T = 2000 the standard deviations of are .184, .178 and .186 respectively while the standard deviations for and are .322 and .311. The similarity between the second and fourth order truncation levels suggests that choosing one of these truncation levels will give us a strong approximation to . Since takes significantly longer to compute, we use the second order approximation in the first two experiments.

This section studies a data set consisting of crimes committed in Chicago between January 2001 and July 2018 ([9]). Point process models have been applied to this data set in the past ([18]). In a missing data setting, in order to validate our model it is important to have a “ground truth” data set. For this reason we limit our study to homicides within the data set. For other crimes recorded data is known to be incomplete ([15, 6]), but we assume that nearly every murder is observed. This allows us to create an incomplete data set by removing murders randomly while still maintaining a ground truth data set for validation.

Figure 3: MSE vs T using different loss functions in (3.6). use the full data X while use the missing data Z. Plots suggest that Taylor series truncations produce nearly identical results to complete loss functions.

The city is divided into 77 community areas and the community area where each murder occurred is recorded. The majority of these community areas experience few murders so we focus on the nine areas with the most murders since 2001. These areas form two clusters: one on the west side of the city and another on the south side. We discretize the events using one week bins, so if a murder occurred in community area i during week t and otherwise. This gives a data matrix which we divide into a train set containing the first 600 weeks in the period, and a test set containing the final 318 weeks. We then create an incomplete data set contains independent realizations of a Bernoulli random variable with mean p = .75.

We learn parameters and using the training set and the full data likelihood . We also learn parameters and using the incomplete train data and the missing data likelihood for various values of

We compare the log-likelihood of these parameters on the test set . The results are shown in Figure 4. The missing data estimates perform nearly as well as the full data estimate when is close to the true value of .75. Note that closely approximates the full data likelihood and the hold out likelihood is substantially worse for compared to for close to .75. In other words, ignoring the missing data entirely gives a weaker estimate than applying the techniques this paper introduces, even when is not a perfect estimate of p. Finally, observe that is more robust to misspecification when compared to when . This is a trend which also appears in Figure 2 and suggests there is value in using conservative estimates of the amount of missed data in practical applications.

Given estimates of A and p we can use density propagation to predict the likelihood of homicides during week n based on observed homicides up to week . We do this for learned from the incomplete data with as well as learned from but with , which corresponds to assuming there is no missing data. We use

Figure 4: Test performance on Chicago crime data. Log-likelihood of events on hold out set using full data with (yellow) and partial data with for various values of For , the proposed estimator performs nearly as well as an oracle estimator with full access to the missing data, and significantly better than a naive method that ignores the possibility of missing data.

particle filtering to construct estimates

These probabilities correspond to the likelihood of homicides during the nth week based on the observations over the first weeks. We construct such estimates for each week in the 318 week test set. As expected assigns higher likelihoods of homicides, with 960 expected homicides in total compared to . As a naive method of correcting for missing data, we divide by a constant scaling factor of 0.75 and report these likelihoods below; by doing this, we ensure that both predictions yield similar average numbers of homicides, so differences in performance between the proposed and naive estimator are not due to a simple difference in averages, but rather because the proposed method is capturing the underlying dynamics of the system.

Figure 5 displays these likelihoods for Community Area 25 (Austin) which has the largest number of homicides recorded during the test period. We use Gaussian smoothing to help visualize the trends. The top panel shows the predicted probability of events using (in red) and the scaled predicted probability of events using (in blue). The bottom panel shows the actual and incompletely observed events during the test period. The true events generally peak at times in which the predicted events for peak. For example, both the predicted event and true event charts have peaks around weeks 60 and 210. In contrast, the predicted events for are nearly constant over time. Since it does not account for the missing data (except via a uniform scaling factor), the network is not able to capture the dynamics of the process and so it cannot predict events with as much precision as

Figure 5: Result of density propagation on Chicago crime data for Community Area 25 (Austin). After a training period used to estimate (proposed estimator) and (naive estimator that doesn’t account for missing data), density propagation is run in subsequent test weeks to predict the likelihood of each community area having a homicide at time n based on observations up to time . The top panel shows the predicted likelihood of a homicide occurring in the Austin community area of Chicago. The network total homicides in the nine community areas during the test period, compared to . The actual number of homicides was 1035. The bottom panel shows the true events as well as the partially observed events (after Gaussian smoothing used for visualization).

We propose a novel estimator for Bernoulli autoregressive models which accounts for partially-observed event data. This model can be used in a variety of contexts in which observed discrete event data exhibits autoregressive structure. We provide mean squared error bounds which suggest that our method can accurately capture a network’s structure in the presence of missing events. Simulations and a real data experiment show that our method yields significant improvement compared with ignoring the missed events and that our method is robust to misspecification of the proportion of missed events. The framework described in this paper suggests a strategy for addressing regression problems with corrupted data in other GLMs, although further work is needed to extend our theoretical analysis beyond binary observations.

[1] A Agarwal, S Negahban, and M. Wainwright. Fast global convergence of gradient methods for high-dimensional statistical recovery. The Annals of Statistics, 40(5):2452– 2482, 2012.

[2] H. Alzer. Sharp bounds for the bernoulli numbers. Archiv der Mathematik, 74(3):207– 211, 2000.

[3] Y. Bao, C. Kwong, P. Peissig, D. Page, and R. Willett. Point process modeling of adverse drug reactions with longitudinal observational data. In Proc. Machine Learning and Healthcare, 2017.

[4] S. Basu and G. Michailidis. Regularized estimation in sparse high-dimensional time series models. Annals of Statistics, 43(4):1535–1567, 2015.

[5] P. Bickel, Y. Ritov, and A. Tsybakov. Simultaneous analysis of Lasso and Dantzig selector. Annals of Statistics, 37(4):1705–1732, 2009.

[6] A. Biderman. Understanding Crime Incidence Statistics. Springer-Verlag, New York, first edition, 1991.

[7] L. Carlitz. Eulerian numbers and polynomials. Mathematics Magazine, 32(5):247–260, 1959.

[8] S. Chen, A. Shojaie, E. Shea-Brown, and D. Witten. The multivariate Hawkes process in high dimensions: Beyond mutual excitation. arXiv preprint arXiv:1707.04928, 2017.

[9] City of Chicago. Crimes - 2001 to present, 2018. https://data.cityofchicago.org.

[10] A. Fletcher and S. Rangan. Scalable inference for neuronal connectivity from calcium imaging. In NIPS, 2014. arXiv:1409.0289.

[11] J. Graham. Missing data analysis: Making it work in the real world. Annual Review of Psychology, 60:549–576, 2009.

[12] E. C. Hall, G. Raskutti, and R. Willett. Inference of high-dimensional autoregressive generalized linear models. arXiv preprint arXiv:1605.02693, 2016.

[13] A. G. Hawkes. Point spectra of some self-exciting and mutually-exciting point processes. Journal of the Royal Statistical Society. Series B (Methodological), 58:83–90, 1971.

[14] A. Jalali and R. Willett. Missing data in sparse transition matrix estimation for subgaussian vector autoregressive processes. In American Control Conference, 2018. arXiv:1802.09511.

[15] L. Langton, M. Berzofsky, C. Kerbs, and H. Smiley-McDonald. Victimizations not reported to the police, 2006-2010. Technical report, Bureau of Justice Statistics, 2012.

[16] T. Le. A multivariate Hawkes process with gaps in observations. IEEE Transactions on Information Theory, 63(3):1800–1811, 2018.

[17] S. Linderman, R. Adams, and J. Pillow. Bayesian latent structure discovery from multi-neuron recordings. In Advances in neural information processing systems, 2016.

[18] S. W. Linderman and R. P. Adams. Discovering latent network structure in point process data. In ICML, 2014. arXiv:1402.0914.

[19] Po-Ling Loh. Statistical consistency and asymptotic normality for high-dimensional robust M-estimators. Annals of Statistics, 42(5):866–896, 2017.

[20] Po-Ling Loh and Martin J Wainwright. High-dimensional regression with noisy and missing data: Provable guarantees with non-convexity. The Annals of Statistics, 40:1637–1664, 2012.

[21] Po-Ling Loh and Martin J Wainwright. Regularized m-estimators with nonconvexity: Statistical and algorithmic theory for local optima. Journal of Machine Learning Research, 16:559–616, 2015.

[22] B. Mark, G. Raskutti, and R. Willett. Network estimation from point process data. arXiv preprint arXiv:1802.09511, 2018.

[23] G. Mohler. Marked point process hotspot maps for homicide and gun crime prediction in chicago. International Journal of Forecasting, 30(3):491–497, 2014.

[24] G Mohler, M Short, Sean Malinowski, Mark Johnson, G Tita, Andrea Bertozzi, and P. Brantingham. Randomized controlled field trials of predictive policing. Journal of the American Statistical Association, 110(512):1399–1411, 2016.

[25] Sahand Negahban, Pradeep Ravikumar, Martin Wainwright, and Bin Yu. A unified framework for high-dimensional analysis of M-estimators with decomposable regularizers. Statistical Science, 27(4):538–557, 2010.

[26] T. Palermo, J. Bleck, and A. Peterman. Tip of the iceberg: Reporting and gender-based violence in developing countries. American Journal of Epidemiology, 179(5):602–612, 2014.

[27] G. Raskutti, M. J. Wainwright, and B. Yu. Restricted eigenvalue conditions for correlated Gaussian designs. Journal of Machine Learning Research, 11:2241–2259, 2010.

[28] C. Shelton, Z. Qin, and C. Shetty. Hawkes process inference with missing data. In Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence, 2018.

[29] NG. Weiskopf and Weng C. Methods and dimensions of electronic health record data quality assessment: Enabling reuse for clinical research. Journal of the American Medical Informatics Association, 20(1):144–151, 2013.

[30] B. Wells, K. Chagin, A. Nowacki, and M. Kattan. Strategies for handling missing data in electronic health record derived data. EGEMS, 1(3):1035, 2013.

[31] H. Xu, D. Luo, X. Chen, and L. Carin. Benefits from superposed Hawkes processes. In Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 2018.

[32] H. Xu, D. Luo, and H. Zha. Learning Hawkes processes from short doubly-censored event sequences. In ICML, 2017. arXiv:1702.07013.

[33] Y. Yang, J. Etesami, N. He, and N. Kiyavash. Online learning for multivariate Hawkes processes. In NIPS, 2017.

[34] K. Zhou, H. Zha, and L. Song. Learning social infectivity in sparse low-rank networks using multi-dimensional Hawkes processes. In Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), 2013.

The supplementary material is organized as follows. In Section A we recall all the notation needed for the proofs. In Section B we restate the main statistical results from the paper and in Section C we provide proofs of these results. In Section D we state the optimization results and give proofs in Section E.

• negative log likelihood of complete data

• Taylor series approximation of

• negative log likelihood of missing data . Loss function is unbiased in the sense that

• Taylor series approximation of

• R

•

• p: fraction of data which is observed

•

Finally, we introduce additional notation which will be helpful in the proofs of Lemmas B.4 and B.5. First let denote the set of all monomials of degree d. We represent an element as a list containing d elements. An element in the list corresponds to the index of a term in the monomial (the list can potentially have repeated elements). For an example, the monomial can be represented as the list (1, 1, 3).

For a polynomial function h we let denote the coefficient of the monomial U in h. Finally we define the order of a list to denote the number of unique elements in the list, so |(1, 2)| = 2 whereas |(1, 1)| = 1.

Example Consider the function . Using all the notation above, we can decompose h as

where with corresponding coefficients and

Remark 1. We next make several observations by applying the notation above to functions which appear in the likelihoods and . First, for a fixed t, m we decompose the following function as a sum of monomials:

and note that

We also have

and we similarly note that

Next consider the function g which appears in the missing data likelihood

where . This observation will be important for our analysis because it allows us to leverage Equation A.1. Similarly we have

where allowing us to use Equation A.2.

Remark 2. Using an identical argument to Lemma B.4 one can show that for . This implies that converges uni- formly on so that is well defined on this ball. Moreover, it implies that converges and so

and satisfies (3.1)

We assume

Theorem B.1 (Accuracy of

for with probability at least

The proof of Theorem B.1 relies on the following supplementary lemmas.

Lemma B.2. Let

Lemma B.3 (Truncation error of

Lemma B.4 (Truncation Error of

Lemma B.5.

with probability at least

Lemma B.6.

with probability at least

Lemma B.7. Let

we have

and

with probability at least

C.1 Proof of Theorem B.1

so by Lemma B.6

Part 2: Setting Up the Standard Equations The next several steps follow standard techniques for regularization in GLMs. Expanding and using the substitution

gives

Rearranging terms yields

Since f(x) = log(1 + exp(x)) is -strongly convex on with we can lower bound Equation (C.1) by . Also note

Using Theorem 2.2 in [12],with probability at least ing these observations to Equation (C.1) gives

and so

Restriction of It remains to lower bound in terms of In order to do this we will rely heavily on the fact that is not an arbitrary vector. Instead we show that must lie in a cone with important properties. In particular, returning to Equation (C.2) and observing that it follows that

Restricted Eigenvalue Condition As mentioned in the previous section, our goal is to lower bound . This is commonly referred to as a restricted eigenvalue condition in the literature, and it is closely related to the restricted strong convexity condition proved in Lemma D.2. In particular, Lemma B.7 guarantees that there exist universal constants

For arbitrary vectors this lower bound can be negative; however, by Equation (C.4)

for a universal constant c. Plugging this in to Equation (C.3) gives that

with probability

C.2 Proof of Lemma B.2

A computation shows that

where the A(q, m) are the Eulerian numbers. The alternating sum of the Eulerian numbers for fixed q can be given as

where th Bernoulli number (see the derivation of Equation 4.8 in [7]). Thus

C.3 Proof of Lemma B.3

For any j we have

where the last two lines use Lemma B.2 along with the fact that

C.4 Proof of Lemma B.4

Differentiating with respect to

We first bound

which is the degree term of . All the terms other than the in are always non-negative and

We conclude that

. This is the sum of the degree d terms of for all . In other words,

as claimed.

C.5 Proof of Lemma B.5

We begin by bounding individual monomials of. We then extend these individual bounds to bounds on the entire expression.

Bounding Individual Monomials Following the notation introduced in Section A a degree monomial of with index U is of the form

Meanwhile, the degree monomial of with index U is of the form

The difference of these monomials is given by

and

We apply Hoeffding’s inequality to conclude that

Thus

with probability at least

Extension to Entire Expression We need to take a union bound so that this holds for all monomials of degree at most log(T). However, since and are binary random variables

whenever

Suppose we have shown that Equation (C.5) holds for all monomials of degree < d. Now to show it holds for all monomials of degree d we only need to show that (C.5) holds for all monomials of degree d that have d distinct terms. The remaining concentrations are already covered by the degree monomials. There aremonomials of degree d with d distinct terms. Hence we need to take a union bound over at most

monomials of degree

We condition on this event and recall that denotes the set of all monomials of degree . Using Equation (C.6), the difference between the degree and can be bounded by

Using Lemma B.2 and Equation A.2 along with the fact that

where the final inequality uses Lemma B.2. Thus we have the bound

with probability at least

C.6 Proof of Lemma B.6

By Lemmas B.3 and B.4 the first two terms can be bounded by

while by Lemma B.5 the final term can be bounded is

with probability at least . We conclude that

with probability at least as claimed.

C.7 Proof of Lemma B.7

We have

By Theorem 2.1 in [12],

Now we define the matrix as follows:

Note that each entry of G is a martingale and

Applying the Azuma-Hoeffding inequality we conclude that for any

Overall we have concluded

with probability at least . Combining Equations (C.7),(C.8) and (C.9) give the desired result.

The main result of this section is the following Theorem.

Theorem D.1. Suppose and for at least rows of where C is a universal constant. Let be any stationary point of where

with probability at least

In order to prove Theorem D.1 we need to introduce notions of Restricted Strong Convexity (RSC) and Restricted Smoothness (RSM) from [1]. To do this we first define the first order Taylor expansion to a function L:

Definition 3 (Restricted Strong Convexity). A loss function L satisfies the RSC condition with parameters

for all

Definition 4 (Restricted Smoothness). A loss function L satisfies the RSM condition with parameters

for all

Lemma D.2. The RSC and RSM conditions are satisfied for with constants

and τ

universal constants.

Combining Lemma D.2 with Theorem 2 in [1] gives the following corollary.

Corollary D.3. Under the conditions of Theorem B.1, let denote the s’th iteration of projected gradient descent using the loss function . There exists some S such that for all s > S we have

E.1 Proof of Lemma D.2

We begin by computing the first order Taylor approximation:

T

We first handle term (1), which is the first order Taylor error for . A computation shows that this is equal to

where again f(x) = log(1 + exp(x)). Since f is -strongly convex on with we can lower bound term (1):

It remains to handle term (2) which is . By the mean value theorem there exists some such that . Thus we can bound term (2) by

By Lemma B.6

with probability . Combining this with our bounds on term 1 in Equations (E.1) and (E.2) gives the final result.

E.2 Proof of Theorem D.1

Since is a stationary point it satisfies

Using this, along with the fact that we get that

Using the RSC condition from Lemma D.2 we conclude

Finally we apply the statistical error bound on from Theorem B.1 along with the triangle inequality to conclude that

Summing over all m and assuming for at least values of m allows us to conclude that

To get the final form of the result, we recall that th iteration of projected gradient descent run with an arbitrary initialization within . In particular, if we initialize at a stationary point then which gives the final form of the Theorem.