Nonstationary Sparse Spectral Permanental Process

Authors: zicheng  sun
yixuan  zhang
zenan  ling
xuhui  fan
feng  zhou
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 08:34:37

Abstract: Existing permanental processes often impose constraints on kernel types or stationarity, limiting the model's expressiveness. To overcome these limitations, we propose a novel approach utilizing the sparse spectral representation of nonstationary kernels. 
This technique relaxes the constraints on kernel types and stationarity, allowing for more flexible modeling while reducing computational complexity to the linear level. 
Additionally, we introduce a deep kernel variant by hierarchically stacking multiple spectral feature mappings, further enhancing the model's expressiveness to capture complex patterns in data. Experimental results on both synthetic and real-world datasets demonstrate the effectiveness of our approach, particularly in scenarios with pronounced data nonstationarity. Additionally, ablation studies are conducted to provide insights into the impact of various hyperparameters on model performance.

Main paper: 
 introduction: 
 Many application domains involve point process data, which records the times or locations of events occurring within a region. Point process models are used to analyze these event data, aiming to uncover patterns of event occurrences. The Poisson process [15] , as an important model in the field of point processes, plays a significant role in neuroscience [4] , finance [11] , criminology [31] , epidemiology [5] , and seismology [8] . Traditional Poisson processes assume parameterized intensity functions, which severely restricts the flexibility of the model. To address this issue, a viable solution is to employ Bayesian nonparametric methods, by imposing a nonparametric Gaussian process (GP) prior on the intensity function, resulting in the Gaussian Cox process [3] . This greatly enhances the flexibility of the model, while also endowing it with uncertainty quantification capabilities.In Gaussian Cox processes, posterior inference of the intensity is challenging because the GP prior is not conjugate to the Poisson process likelihood that includes an intensity integral. Furthermore, to ensure the non-negativity of the intensity, we need to use a link function to transform the GP prior. Commonly used link functions include exponential [20; 12] , sigmoid [1; 9; 6] , square [18; 7; 29] , ReLU [16] , and softplus [24] . Among these, using the square link function, i.e., the permanental process [19] , allows for the analytical computation of the intensity integral [18; 33; 29] , thus receiving widespread attention in recent years. For this reason, this work focuses on the permanental process.Currently, the permanental process faces three issues: (1) The permanental process inherits the notorious cubic computational complexity of GPs, making it impractical for use with a large amount 
permanental process: 
 For a Poisson process, if x ∈ X ⊂ R D , the intensity function λ(x) = lim δx→0 E[N ([x, x + δ x ])]/δ x can be used to quantify the rate of events occurring at location x. A Cox process can be viewed as a Poisson process whose intensity function itself is a random process. The Gaussian Cox process employs the GP to model the intensity function: λ(x) = l • f (x) where f is a GP function and l : R → R + is a link function to ensure the non-negativity of the intensity function.The nonparametric nature of GP enables the Gaussian Cox process to be flexible in modeling events in space or time. However, the posterior inference of Gaussian Cox process is challenging. According to the likelihood function of a Poisson process:EQUATIONand the Bayesian framework, the posterior of latent function f can be written as:EQUATIONIn Eq. ( 2), there are two integrals. The first one X λ(x)dx cannot be solved due to the stochastic nature of f . Additionally, the integral over f in the denominator is a functional integral in the function space, which is also infeasible. This issue is the well-known doubly-intractable problem.The link function l : R → R + has many choices, such as exponential, sigmoid, square, ReLU, and softplus. This work focuses on the utilization of the square link function, also known as the permanental process. Furthermore, if we set the kernel in the GP prior to be stationary: k(x 1 , x 2 ) = k(x 1 -x 2 ), then the model is referred to as a stationary permanental process.The sparse spectral representation is a common method for the low-rank approximation of kernels. This approach is based on Bochner's theorem. Theorem 3.1. [2] A stationary kernel function k(x 1 , x 2 ) = k(x 1 -x 2 ) : R D → R is bounded, continuous, and positive definite if and only if it can be represented as:k(x 1 -x 2 ) = R D exp(iω ⊤ (x 1 -x 2 ))dµ(ω),where µ(ω) is a bounded non-negative measure associated to the spectral density p(ω) = µ(ω) µ(R D ) .Defining ϕ(x) = exp(iω ⊤ x), we have k(x 1 -x 2 ) = σ 2 E p(ω) [ϕ(x 1 )ϕ(x 2 ) * ] where * denotes the complex conjugate and σ 2 = µ(R D ). If we utilize the Monte Carlo to sample {ω r } R r=1 independently from p(ω), then the kernel can be approximated byEQUATIONwhere ϕ r (x) = exp(iω ⊤ r x),Φ (R) (x) = σ √ R [ϕ 1 (x), • • • , ϕ R (x)] ⊤ , which corresponds to the random Fourier features (RFF) proposed by [25] . It can be proved that Eq. ( 3) is equivalent to a 2R-sized trigonometric representation (proof is provided in Appendix A):EQUATION)When we specify the functional form of the kernel (spectral measure µ), the RFF method can learn the parameters of the kernel (spectral measure µ) from the data. However, when the functional form of the kernel (spectral measure µ) is unknown, RFF can not flexibly learn the kernel (spectral measure µ) from the data. Another approach is to treat the frequencies {ω r } R r=1 as kernel parameters. By optimizing these frequencies, we can directly learn the kernel from the data in a flexible manner, which corresponds to the sparse spectrum kernel method proposed by [17] . 
our model: 
 The permanental process inherits the cubic computational complexity of GPs, rendering it impractical for large-scale datasets. To efficiently conduct posterior inference, many works employed various low-rank approximation methods [18; 7; 33] . These methods require certain standard types of kernels, such as the squared exponential kernel, polynomial kernel, etc., to ensure that the intensity integral has analytical solutions. This limitation severely affects the flexibility of the model. To address this issue, sparse spectral permanental processes are proposed in recent years [14; 29] . The advantage of this method lies in its ability to analytically compute the intensity integral, while not requiring restrictions on the kernel's functional form. However, the drawback is that it is only applicable to stationary kernels. Therefore, a natural question arises: can we further generalize the sparse spectral permanental processes to not only remove restrictions on the kernel's functional form but also make them applicable to nonstationary kernels?Existing sparse spectral permanental processes rely on the prerequisite of Bochner's theorem, which limits the extension of the sparse spectral method to nonstationary kernels. Fortunately, for more general nonstationary kernels, the spectral representation similar to Bochner's theorem also exists. Theorem 4.1. [37] A nonstationary kernel function k(x 1 , x 2 ) : R D × R D → R is bounded, continuous, and positive definite if and only if it can be represented as:k(x 1 , x 2 ) = R D ×R D exp i(ω ⊤ 1 x 1 -ω ⊤ 2 x 2 ) dµ(ω 1 , ω 2 ),where µ(ω 1 , ω 2 ) is a Lebesgue-Stieltjes measure associated to some bounded positive semi-definite spectral density p(ω1 , ω 2 ) = µ(ω1,ω2) µ(R D ,R D ) .Similarly, we can use the Monte Carlo method to approximate it. However, it is worth noting that, to ensure that the spectral density p(ω 1 , ω 2 ) is positive semi-definite, we must first symmetrize it, i.e., p(ω 1 , ω 2 ) = p(ω 2 , ω 1 ); and secondly introduce diagonal components p(ω 1 , ω 1 ) and p(ω 2 , ω 2 ). We can take a general density function g on the product space and then parameterize p(ω 1 , ω 2 ) in the following form to ensure its positive semi-definite property [26; 32] :EQUATIONThen, we can obtain the sparse spectral feature representation of any nonstationary kernel:EQUATIONEQUATIONIf the GP's kernel in the permanental process (Eq. ( 2)) adopts the above sparse spectral nonstationary kernel, then we refer to this model as NSSPP.We can further enhance the expressive power of the nonstationary kernel by stacking multiple spectral feature mappings hierarchically to construct a deep kernel. The spectral feature mapping φ (R) (x) has another equivalent form: The spectral feature mapping in Eq. ( 8) can be viewed as a more complex single-layer neural network: it employs two sets of weights (ω) and biases (b) to linearly transform the input, followed by the application of cosine activation functions, and then adds the results to obtain a single output dimension. By stacking multiple layers on top of each other, we naturally construct a deep nonstationary kernel:φ (R) (x) = σ √ 2R [EQUATIONwhere Ψ (R) (x) is the spectral feature mapping with L layers and φ (R) l denotes the l-th layer § . This actually recovers the deep spectral kernel proposed by [36] . If the GP's kernel in the permanental process (Eq. ( 2)) adopts the above deep sparse spectral nonstationary kernel, then we refer to this model as DNSSPP.DNSSPP exhibits enhanced expressiveness compared to NSSPP due to its deep architecture. In subsequent sections, we consistently use Ψ (R) as the spectral feature mapping. When it has a single layer, the model corresponds to NSSPP; when multiple layers, the model corresponds to DNSSPP. 
inference: 
 For the posterior inference of (D)NSSPP, this work employs a fast Laplace approximation exploiting the sparse spectral representation of nonstationary kernels. Other inference methods, such as variational inference, are also feasible. Specifically, we extend the method proposed in [29] from stationary kernels to nonstationary kernels.When we approximate a kernel using Eq. (6) or Eq. ( 9), the GP prior can also be approximated as:EQUATIONAdditionally, we set the intensity function as λ(x) = f (x) + α 2 . Following [14] , we introduce an offset α to mitigate the nodal line problem caused by the non-injectiveness of λ = f 2 . This issue results in the intensity between positive and negative values of f being forced to zero, despite the underlying intensity being positive.Combing the likelihood in Eq. ( 1) and the equivalent GP prior in Eq. (10) , we obtain the joint density:EQUATION)where C is a constant and Θ denotes the model hyperparameters. The intensity integral term can be further expressed as:EQUATION§ For simplicity, here we assume all layers have the same width R, but it can also vary.where |X | is the window size, M is a R × R matrix, and m is a R-sized vector with each entry:EQUATIONIt is worth noting that the intensity integral is potentially analytically solvable. For a single-layer spectral feature mapping (NSSPP), analytical solutions for M and m exist (see Appendix C). However, for multiple-layer mapping (DNSSPP), analytical solutions are not possible. In the following, we analytically solve M and m for NSSPP, while using numerical integration for DNSSPP.We use the Laplace's method to approximate the posterior p(β| {x i } N i=1 , Θ). The Laplace's method approximates the posterior with a Gaussian distribution by performing a second-order Taylor expansion around the maximum of the log posterior. Following the common practice:EQUATION)where β is the mode of the true posterior and Q is the negative inverse Hessian of the true posterior evaluated at the mode:EQUATIONThe mode β is typically not directly solvable based on Eq. ( 15), so we employ optimization methods to find it. Then, the precision matrix Q -1 can be computed analytically (see Appendix D).The model hyperparameter Θ consists of the kernel parameters σ, ω, b (if there are multiple layers, these hyperparameters include parameters from all layers) and the bias of the intensity function α. We can learn these hyperparameters from the data by maximizing the marginal likelihood: 16) Substituting Eqs. ( 11) and ( 14) into Eq. ( 16), we can obtain a marginal likelihood approximation. Optimizing this yields the optimal hyperparameter Θ.log p({x i } N i=1 |Θ) = log p({x i } N i=1 , β|Θ) -log p(β| {x i } N i=1 , Θ) ≈ log p({x i } N i=1 , β|Θ) -log N (β| β, Q). (Similar to common low-rank approximation methods for kernels, the computational complexity of the proposed inference method is reduced from O(N 3 ) to O(N R 2 ), where R ≪ N , i.e., linearly with N . A complete algorithm pseudocode is provided in Appendix E.After obtaining the posterior β| {x i } N i=1 , Θ ∼ N (β| β, Q), we can compute the posterior of the intensity function. For any x * ∈ X , the predictive distribution of f(x * ) + α is f (x * ) + α| {x i } N i=1 , Θ ∼ N (µ(x * ), σ 2 (x * )), where µ(x * ) = β⊤ Ψ (R) (x * ) + α, σ 2 (x * ) = Ψ (R) (x * ) ⊤ QΨ (R) (x * ).Then, the intensity λ(x * ) = f (x * ) + α 2 can be proven to have the following mean and variance:EQUATION) 
experiments: 
 In this section, we mainly analyze the superiority in performance of (D)NSSPP over baseline models on both synthetic and real-world datasets, as well as the impact of various hyperparameters. We perform all experiments using the server with two GPUs (NVIDIA TITAN V with 12GB memory), two CPUs (each with 8 cores, Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10GHz), and 251GB memory. We employ several previous works as baselines for comparison. These baselines include VBPP which utilizes the inducing points method [18] , LBPP employing the Nyström approximation [33] , as well as SSPP and GSSPP utilizing the Fourier representation [29] . We use GSSPP based on Gaussian kernel, Matérn kernel with parameter 1/2 and 5/2 respectively, and we denote them as GSSPP, GSSPP-M12 and GSSPP-M52. Besides, we offer a baseline that employs stacked mappings of stationary kernels which we name it Deep Sparse Spectral Permanental Process (DSSPP). All of these methods utilize stationary kernels. As far as we know, there is currently no work employing nonstationary kernels in permanental processes. To compare with nonstationary baselines, we implement a nonstationary permanental process using the nonstationary spectral mixture kernel [26] , referred to as NSMPP.We employ two metrics to evaluate the performance of various baselines. One is the expected log-likelihood on the test data, denoted as L test , with details provided in Appendix F. Additionally, when the ground-truth intensity is available, e.g., for the synthetic data, the root mean square error (RMSE) between the expected posterior intensity and the ground truth serves as another measure.Datasets To better compare the performance of (D)NSSPP with baseline models on point process data in both stationary and nonstationary scenarios, we simulated a stationary and a nonstationary permanental process on the interval [0, 10], respectively. Specifically, we assume two kernels:k 1 (x 1 , x 2 ) = exp(- 1 2 ∥x 1 -x 2 ∥ 2 ), k 2 (x 1 , x 2 ) = ( x ⊤ 1 x 2 100 + 1) 3 exp(- 1 2 ∥x 1 -x 2 ∥ 2 ),where k 1 is a stationary Gaussian kernel and k 2 is a nonstationary kernel obtained by multiplying a Gaussian kernel with a polynomial kernel. We use these two kernels to construct a stationary GP and a nonstationary GP, respectively, and randomly sample two corresponding latent functions. Based on λ(x) = (f (x) + 2) 2 , we construct two corresponding intensity functions. Finally, we use the thinning algorithm [22] to simulate two sets of synthetic data. For each intensity, we simulate ten datasets and use each dataset alternately as the training set and the remaining ones as the test sets.Setup For NSSPP, the number of frequencies (network width) is set to 50. For DNSSPP, we experimented with three different configurations: DNSSPP-[50,30], DNSSPP-[100,50], and DNSSPP- [30, 50, 30] . Each number represents the width of the corresponding network layer. Therefore, DNSSPP-[50,30] implies a model with two layers, where the first layer has a width of 50 and the second layer has a width of 30, and so on. For baseline models, we need to set the number of rank for kernel low-rank approximation. We set the number of inducing points to 50 for VBPP, the number of eigenvalues to 50 for LBPP, and the number of frequencies to 50 for SSPP, GSSPP, and NSMPP. For DSSPP, to facilitate comparison with DNSSPP, we adopt the same layer configurations.The fitting results of intensity functions for all models are depicted in Fig. 1 . To quantitatively compare the fitting performance of different methods in both stationary and nonstationary scenarios, we employ L test and RMSE as metrics. The comparison results are presented in Table 1 . DNSSPP with three different configurations exhibits significant advantages on both stationary and nonstationary datasets. NSSPP and NSMPP show better performance in the nonstationary scenario compared to most stationary baselines. This is reasonable, as stationary baselines face challenges in accurately modeling nonstationary data. Because stationary kernels are a subset of nonstationary kernels, nonstationary models typically perform well on stationary data as well. Compared to relatively shallow models, DNSSPP achieves better performance but also incurs a longer running time.Datasets In this section, we use three sets of common real-world datasets to evaluate the performance of our method and the baselines. Coal Mining Disaster [13] : This is a 1-dimensional dataset containing 191 incidents that occurred between March 15, 1875, and March 22, 1962 . Each event represents a mining accident that killed more than 10 people. Redwoods [27] : This is a 2-dimensional dataset that describes the distribution of redwoods in an area, consisting of 195 events. Porto Taxi [21] : This is a large 2-dimensional dataset containing the tracks of 7,000 taxis in Porto during 2013/2014. We focus on the area with coordinates between (41.147, -8.58) and (41.18, -8.65), and extract 3,000 pick-up points as our dataset. For each dataset, we randomly partition the data into training set and test set of approximately equal size.Setup For real datasets, because we do not know the ground-truth intensity, we use L test as the evaluation metric. Since NSSPP generally performs worse than DNSSPP, we report only the results of DNSSPP in this section. For DNSSPP, we consistently use the same three configurations as in the synthetic data. For baseline models, we choose different numbers of rank for kernel low-rank approximation for 1-dimensional and 2-dimensional datasets. Specifically, for the 1-dimensional dataset (Coal), we set the number of rank, i.e., the number of inducing points for VBPP, the number of eigenvalues for LBPP, and the number of frequencies for SSPP, GSSPP, and NSMPP, to 10. For the 2-dimensional datasets (Redwoods and Taxi), we set the number of rank to 50. We discuss the reasonableness of the above setup in Appendix H.Results Fig. 2 shows the fitting results of the intensity functions from LBPP and DNSSPP on the Redwoods and Taxi datasets. The quantitative comparison results are shown in Table 2 . Since real-world data is more or less nonstationary to varying degrees, nonstationary models generally perform better on real-world data. DNSSPP secures both first and second places on the Coal and Taxi datasets and achieves second place on the Redwoods dataset. For simpler datasets, such as Coal and Redwoods, which have relatively simple data patterns, the performance of DNSSPP is comparable to that of the simpler baselines. However, for the more complex Taxi dataset, DNSSPP demonstrates significant advantages over the simpler baselines. This indicates that DNSSPP has a clear advantage in handling data with more complex patterns, such as large-scale or high-dimensional datasets.In this section, we analyze the impact of four configuration hyperparameters on the performance of DNSSPP: the network width and depth, as well as the number of epochs and learning rate during the update of hyperparameter Θ.We investigate DNSSPP's performance with varying network sizes by adjusting width and depth on nonstationary synthetic data. In our experiments, we maintain consistent width across layers. Results are illustrated in Fig. 1c . As the number of layers and the width of the network increase, DNSSPP's performance initially improves and then declines, reflecting the phenomenon of overfitting as the network size increases. Regarding the number of layers, DNSSPP generally exhibits underfitting when there is only a single layer, while overfitting tends to occur with four layers. In terms of network width, due to the relatively small size of the synthetic data, overfitting becomes apparent when the width is large. However, for larger datasets, a more complex network structure should be chosen.During the inference process, we need to perform numerical optimization on the hyperparameter Θ, implying the need to determine the epoch and learning rate for optimization. We conduct experiments on the nonstationary synthetic data to investigate the coordination of epoch and learning rate. Results are illustrated in Fig. 1d . When the learning rate is too low or the number of epochs is insufficient, the update of Θ is sluggish, leading to inferior performance. However, when the learning rate is too high or the number of epochs is excessive, model performance generally declines as well. This is because our model training is a bi-level optimization: the inner level is the Laplace approximation, and the outer level is to update the hyperparameter Θ. An excessively high learning rate or too many epochs can cause Θ to converge too quickly, making the overall model training more prone to getting stuck in local optima, thus resulting in suboptimal performance.Question: Does the paper specify all the training and test details (e.g., data splits, hyperparameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?Answer: [Yes] Justification: We specify all the training and test details necessary to understand the results. See the experimental setup in the Experiments section.Guidelines:• The answer NA means that the paper does not include experiments.• The experimental setting should be presented in the core of the paper to a level of detail that is necessary to appreciate the results and make sense of them. • The full details can be provided either with the code, in appendix, or as supplemental material. 
limitations: 
 NSSPP not only removes restrictions on the kernel's functional form but also makes it applicable to nonstationary kernels. DNSSPP further enhances its expressive power through a deep architecture. However, the deep architecture of DNSSPP prevents analytical computation of the intensity integral, necessitating numerical integration. This reliance on numerical integration imposes a limitation on the model's computational speed.Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?Answer: [Yes] Justification: We report the statistical significance of the experiments. In the tables presenting the experimental results, we provide the standard deviations across multiple runs, while in some figures, we also include confidence intervals.Guidelines:• The answer NA means that the paper does not include experiments.• The authors should answer "Yes" if the results are accompanied by error bars, confidence intervals, or statistical significance tests, at least for the experiments that support the main claims of the paper.• The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions). • The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.) • The assumptions made should be given (e.g., Normally distributed errors). • It should be clear whether the error bar is the standard deviation or the standard error of the mean. Justification: The research conducted in the paper conform with the NeurIPS Code of Ethics.Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
conclusions: 
 In this study, we introduced NSSPP and its deep kernel variant, DNSSPP, to address limitations in modeling the permanental process. This approach overcomes the stationary assumption, allowing for flexible learning of nonstationary kernels directly from data without restricting their form. By leveraging the sparse spectral representation of nonstationary kernels, we achieved a low-rank approximation, effectively reducing computational complexity. Our experiments on synthetic and real-world datasets demonstrated that (D)NSSPP achieved competitive performance on stationary data while excelling on nonstationary datasets. This study underscores the importance of considering nonstationarity in the permanental process and demonstrates the efficacy of (D)NSSPP in this context. 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [No] Justification: This paper presents work whose goal is to advance the field of machine learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: The paper poses no such risks.Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: All code, models, and datasets mentioned in the text are appropriately cited with their original papers.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 
new assets: 
 Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?Answer: [Yes]Justification: New assets introduced in the paper, such as code, are well documented. The documentation is provided alongside the assets in the supplementary material.Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file. 
crowdsourcing and research with human subjects: 
 Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects. • Including this information in the supplemental material is fine, but if the main contribution of the paper involves human subjects, then as much detail as possible should be included in the main paper. • According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector. 
institutional review board (irb) approvals or equivalent for research with human subjects: 
 Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?Answer: [NA]Justification: The paper does not involve crowdsourcing nor research with human subjects.Guidelines:• The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.• Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: The derivation of Eq. ( 4) is provided below:EQUATIONIn the derivation, due to the symmetry of the spectrum distribution, it can be assumed that 2R points are symmetrically sampled, which allows the elimination of the imaginary part in the equation without altering the expectation.The derivation of Eq. ( 8) is provided below:k(x 1 , x 2 ) ≈ σ 2 4R R r=1 exp i(ω ⊤ 1r x 1 -ω ⊤ 2r x 2 ) + . . . + exp i(ω ⊤ 2r x 1 -ω ⊤ 2r x 2 ) = σ 2 4R R r=1 cos(ω ⊤ 1r x 1 -ω ⊤ 2r x 2 ) + . . . + cos(ω ⊤ 2r x 1 -ω ⊤ 2r x 2 ) (19) = σ 2 4R R r=1 i,j=1,2 [cos(ω ⊤ ir x 1 -ω ⊤ jr x 2 )],where the second line is because the kernel is real-valued, so we can safely eliminate the imaginary part. Additionally, considering the law of total expectation, we haveEQUATIONwhere b ∼ Uniform(0, 2π). And consider the periodicity of cosine function, we haveEQUATIONso Eq. ( 20) is finally equal to 0. Therefore, we can further obtain another equivalent form of the spectral feature mapping, through the following procedure:k(x 1 , x 2 ) ≈ σ 2 4R R r=1 (i,j)={1,2} 2 [cos(ω ⊤ ir x 1 -ω ⊤ jr x 2 ) + cos(ω ⊤ ir x 1 + ω ⊤ jr x 2 + b ir + b jr )] = σ 2 4R R r=1 (i,j)={1,2} 2 2 cos(ω ⊤ ir x 1 + b ir ) cos(ω ⊤ jr x 1 + b jr ) (22) = φ (R) (x 1 ) ⊤ φ (R) (x 2 ),where {b 1r , b 2r } R r=1 are uniformly sampled from [0, 2π].The intensity function of the permanental process has the following expression λ(x) = (β ⊤ Ψ (R) (x)+ α) 2 . For ease of derivation, we adopt the form of Eq. ( 7) for the spectral feature mapping Ψ (R) (x) in subsequent steps, i.e., a 2R-sized vector. When Ψ (R) (x) takes the form of Eq. ( 8), i.e., an R-sized vector, the calculation result of the intensity integral remains the same. The intensity integral on X is:X λ(x)dx = X (β ⊤ Ψ (R) (x) + α) 2 dx = β ⊤ X Ψ (R) (x)Ψ (R) (x) ⊤ dxβ + 2αβ ⊤ X Ψ (R) (x)dx + X α 2 dx = β ⊤ Mβ + 2αβ ⊤ m + α 2 |X |,whereEQUATIONIt is worth noting that M and m are analytically solvable. The analytical solution for the stationary case is provided in Appendix C.2 in [29] , and we provide the analytical expression for the nonstationary case here.For M, the product of two spectral feature mappings is:EQUATIONWe can further transform the product term of trigonometric functions:EQUATIONWithout loss of generality, we consider the integral domain X as [-d, d] D . Then we can compute the integral of the above expression analytically. Since we only consider the D = 1 and D = 2 point process data in this paper, we specifically demonstrate the analytical integral expressions in the one-dimensional and two-dimensional cases here.One-dimension:EQUATIONEQUATION)Two-dimension:EQUATIONEQUATIONwhere η 1 denotes the first ordinate of η, and η 2 denotes the second ordinate of η. It should be noted that theoretically, the integral in Eq. ( 28) only has two possible results as explained in the right side. The analytical solution of the integral of Eq. ( 25) can be obtained by replacing η with ω ai -ω bj or ω ai + ω bj . Note that the case a = b and i = j corresponds to the second case in Eq. ( 26) and Eq. ( 28), which means the computation of the diagonal entries of M is different from that elsewhere.Next, we can compute m:EQUATIONIt is easy to see that if D = 1, m can be analytically computed using Eq. ( 26) and Eq. ( 27), while if D = 2, m can be analytically computed using Eq. ( 28) and Eq. ( 29). It is worth noting that the analytical expression of the intensity integral above applies only to NSSPP (single-layer network).For DNSSPP (multi-layer network), the nested structure of trigonometric functions leads to M and m lacking analytical solutions, and we need resort to numerical integration.The gradient and the Hessian matrix of the posterior of β are:EQUATIONEQUATION)In theory, setting Eq. ( 31) to 0 can solve for the mode of the posterior, denoted as β. However, in practice, we cannot analytically solve it, so we use numerical optimization to obtain β. Then according to Eq. ( 32), we can obtain the precision matrixQ -1 = -∇ 2 β log p({x i } N i=1 , β|Θ)| β= β .The pseudocode of our inference algorithm is provided in Algorithm 1.Algorithm 1: Inference for (D)NSSPP Define the network depth L and width R, and initialize the hyperparameter Θ; for Iteration do Calculate M and m by Eq. ( 13); Update the mode β by maximizing Eq. ( 14); Update the covariance matrix Q by Eq. ( 15); Update the hyperparameters Θ by maximizing Eq. ( 16); end For any x * , output the predicted mean and variance of λ(x * ) by Eq. ( 17).The computation of expected test log-likelihood is provided in Appendix E.2 in [29] . For the completeness of the article, we restate the specific derivation process here. Considering the posterior of the intensity function, we have the approximation of the expected test log-likelihood:EQUATIONderived from Eq. ( 1).First, we compute the first term in Eq. ( 33):EQUATIONwhere M and m are defined as Eq. ( 13).Then, we consider the second term in Eq. (33) . Applying the computational method used by [18] , we can obtain the following expression:EQUATION= -G(-EQUATIONwhere Const ≈ 0.57721566 is the Euler-Mascheroni constant. G is defined asG(z) = 2z ∞ j=0 j!z j (2) j (1/2) j ,where (•) j denotes the rising Pochhammer series:(a) 0 = 1, (a) k = a(a + 1)(a + 2)...(a + k + 1).In practical computation, we first prepare a lookup table for the G function and then obtain numerical approximations through linear interpolation.In this section, we provide additional experimental results for the real data. Specifically, we present the fitting results of the intensity functions from all models on the Coal dataset in Fig. 3 ; the fitting results of the intensity functions from VBPP, SSPP, GSSPP, and NSMPP on the Redwoods dataset in Fig. 4 ; and the fitting results of the intensity functions from VBPP, SSPP, GSSPP, and NSMPP on the Taxi dataset in Fig. 5 . H Some Discussion about Experimental Setup for Real DataFor the real data, we chose different numbers of ranks for the 1-dimensional dataset (10 for Coal) and the 2-dimensional datasets (50 for Redwoods and Taxi) for all baseline methods, except for DNSSPP. The reason is that the Coal dataset is very small (only 191 events) and has a simple data pattern. A higher number of ranks would not further improve performance but would significantly increase the algorithm's running time. Therefore, we set the number of ranks to 10 for Coal, but 50 for Redwoods and Taxi because the latter two datasets are larger and more complex, requiring a higher number of ranks to improve performance.To demonstrate this, we further increased the number of ranks for all baseline models to 50 on the Coal dataset and compared the result to rank=10. The result is shown in Table 3 . As can be seen from the results, when we increased the number of ranks from 10 to 50, the performance of all baseline models did not improve significantly, but the algorithm's running time increased substantially.In summary, for the same dataset, the number of ranks for all baseline methods is kept consistent, except for DNSSPP because it is not a single-layer architecture. The choice of the number of ranks for the baseline methods was made considering the trade-off between performance and running time. Further increasing the number of ranks would not significantly improve model performance but would result in a substantial increase in computation time.Table 3 : The performance of L test and runtime for all baselines on the Coal dataset, with the number of rank set to 10 and 50. When we increase the number of ranks from 10 to 50, the performance of all baseline models does not improve significantly, but the algorithm's running time increases substantially.Coal (rank = 10)Coal (rank = 50) Justification: We provide the data and code in the supplemental material to reproduce the main experimental results.LGuidelines:• The answer NA means that paper does not include experiments requiring code. • Providing as much information as possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted. 


Decision: This paper proposes an approach to using sparse spectral representations of non-stationary kernels to improve the expressiveness of permanental processes, providing more flexible modeling via less constrained kernels. The majority of the reviewers initially assessed this paper as having clear motivation, strong writing, strong novelty, effective experiment designs, promising results and considerable significance. Following the author response and discussion, all of the reviewers were in favor of accepting the paper and all concerns were addressed. The scores support an accept decision. The authors are encouraged to incorporate clarifications from the discussion and reviewer suggestions into the final manuscript.

Reviews: review: paper_summary: paper_summary: The paper introduces an approach to modeling permanental processes by utilizing a sparse spectral representation of nonstationary kernels, termed as Nonstationary Sparse Spectral Permanental Process (NSSPP) and its deep kernel variant (DNSSPP). This method addresses the limitations of traditional permanental processes which often require specific kernel types and assume stationarity, thus restricting the model's flexibility and computational efficiency. The deep kernel variant (DNSSPP) uses hierarchically stacked spectral feature mappings, significantly enhancing the model's ability to capture complex patterns in data. The paper presents experimental results on both synthetic and real-world datasets, showing that NSSPP and DNSSPP perform competitively on stationary data and show improvements on nonstationary datasets.
strength_weakness: strength_weakness: - The experimental results (Table 1) indicate that performance enhancements are observed mainly with the deep kernel variant (DNSSPP) when compared to stationary kernels. It would be beneficial for the authors to include a comparison with a baseline that employs stacked mappings of stationary kernels to more clearly ascertain the source of these improvements.
questions: questions: - Performance comparison with a baseline that employs stacked mappings of stationary kernels
limitations: limitations: Yes.
paper_summary: paper_summary: Point process models are a commonly used technique for analysis of event-based data. Gaussian Cox processes are an example which use GPs to model the intensity function in a Cox process, which itself is a specific case of a Poisson process where the intensity function is a stochastic process. Generally, these types of processes are used to model rates of events occurring at given input locations. The link function within the Gaussian Cox process is a design choice; if we choose a square link function and a stationary kernel for the GP prior, this yields a stationary permanental process. The authors aim to address three key problems with these existing approaches; the cubic complexity of exact GP inference, the requirement to use kernels which allow for analytical solutions to the intensity integral, and the reliance on shallow kernels which limit the expressiveness and flexibility of the model. Introduced is a new form of permanental process model which uses a sparse spectral representation to improve computational efficiency, allows for nonstationary kernels to be used without restriction on the form of said kernel, and allows deep kernels to be formed by composing spectral feature mappings.
strength_weakness: strength_weakness: -	This is just a reoccurring minor typo but I’m fairly certain it’s the Porto taxi dataset (as in the city of Porto, Portugal), not Proto, so just edit all mentions of this.
-	I wasn’t entirely sure of the rationale behind some of the baseline selections, for example, for some of the real-world experiments, you use as low as 10 inducing points/frequencies for the baselines, but use more than this for the DNSSPP (‘same configurations as for the synthetic data’). It would be good to have some discussion of this, and/or how the performance of the proposed methods compares to that of the baselines as we increase the number of parameters. Obviously what constitutes a “fair” comparison between models is highly subjective and depends on various factors such as computational complexity as well as parameter count, but I just wanted a little more context on this.
-	I’m aware that there has been some work by  on understanding and exploring overfitting in the context of deep kernel learning (Promises and Pitfalls of DKL, Ober et al, 2021), it would be useful to touch on some work from this area during your discussion on overfitting just to clarify to the readers that this phenomenon is not arising solely due to the size of the datasets you are modelling.
questions: questions: -	As mentioned earlier, I think the rationale for using of a very small number of inducing points etc. for some of the baselines needs to be clarified in the text; is it the case that if you increase this number to 50+ that the baselines begin to outperform the models proposed by the authors?
-	You mention other variational methods are feasible, was there a certain specific reason that you decided upon a Laplace approximation-based approach?
-	See earlier point about discussing the impact of the deep kernel learning framework in general on overfitting.
-	See typos in weaknesses section, please fix.
limitations: limitations: The Limitations section of the paper briefly discusses the fact that the deep formulation of the DNSSPP does not allow for analytical computation of the intensity integral, which necessitates computationally inefficient numerical integration. There are also some practical considerations/limitations discussed elsewhere such as in Section 6.5.2.
paper_summary: paper_summary: Point process are generative models for datasets of points in (typically 1D, 2D, 3D Euclidean) space, e.g. datasets of rain drops, taxi pickup locations, and a Poisson process consists of an intensity function $ \lambda: \mathcal{X} \to \mathbb{R}^+ $ that roughly shows how likely a point is to appear at $x$. The number of points $N_S$ in a region $S\subset \mathcal {X}$ is Poisson distributed $N_S\sim \text{Poi}[\lambda_S]$ with mean parameter $\lambda_S = \int_S \lambda (x) dx$.

A Poisson Process with a (squared for positive output) Gaussian Process prior over $\lambda(x)$ is a permanental process, from point clouds in 1D, 2D or 3D space, one can learn smooth underlying intensity functions. There has been much study on fitting GPs to learn intensity functions using well known GP tools such as inducing point methods and this work builds closely upon the work of [27] which uses the popular method of Random Fourier features for stationary kernels.

This paper makes two novel contributions, firstly using Random Fourier features for non-stationary kernels and  deriving the corresponding approximate marginal likelihood. Then secondly using deep kernel by nesting Fourier features that increases model capacity at the cost of requiring using numerical integrations for the approximate marginal likelihood. Experiments are performed and positive a range of results and model hyperparameters showing results both with and without improvements over baselines.

[27] [Sellier 2023](https://proceedings.mlr.press/v206/sellier23a/sellier23a.pdf)
strength_weakness: strength_weakness: My only concern is the lack of substantial novel contribution, particularly when compared to the work of [27] that introduced GSSPP.

- Replacing the approximation based on Bochner's Theorem with the approximation based on Yanglon 1987 [37] significantly increases model capacity and flexibility but, unless I am mistaken, appears to be a rather small change from GSSPP, we sample over pairs of frequencies instead.
- consequently, deriving the Laplace approximate marginal Likelihood for NSSPP also appears to be not too difficult given the result for GSSPP, use trig identities to combine frequencies and I believe the results follow from the stationary case.
- using deep features understandably helps model fitting however the analytic result of NSSPP is lost.

Overall, this felt like GSSPP with different kernels, which introduced some difficulties that seemingly had fairly simple solutions.


I feel that it is an elegant although not-too-difficult extension of GSSPP and I do not feel the contribution is significant enough for publication at NeurIPS. I enjoyed the paper and I feel the contribution is certainly very respectable and well designed and demonstrates the efficacy of the method.


[27] [Sellier 2023](https://proceedings.mlr.press/v206/sellier23a/sellier23a.pdf)

[33] [Fast Bayesian Intensity Estimation for the Permanental Process, Walder and Bishop, ICML 2017](https://proceedings.mlr.press/v70/walder17a.html)

[37] [Correlation Theory of Stationary and Related Random Functions, Akiv aYaglom](https://link.springer.com/book/9781461290865)
questions: questions: - how is the marginal likelihood for DNSSPP computed? (The paper just says "numerically" and I couldn't see details in the appendix)
- do the error bars represent a one/two standard errors or standard deviation over test set likelihoods?
limitations: limitations: The authors have described the method in full.

Comments: title: Thanks
Thank you for the response
Thanks
Rebuttal response
comment: Thank you very much for your recognition and support. We will revise the paper according to your suggestions. Thank you once again for your constructive feedback and for increasing your rating.
Thank you for the response.

It seems all reviewers agree on technical aspects and only I differ by my subjective judgement on the impact for which it seems I may have been too harsh!

I am have to raised my score.
Thank you very much for your thoughtful review and for taking the time to carefully consider our responses. We greatly appreciate your kind words and are pleased that you find our work to be well-presented, thorough, and novel. 

We will certainly include the additional results in the updated manuscript, as you suggested. Your feedback is valuable in improving the quality of our work, and we are grateful for your support.

Thank you once again for your constructive feedback and for increasing your rating to accept.
I appreciate the time and effort taken by the authors on the response to my review, and for addressing each of my concerns in turn. This is a well presented, thorough and novel piece of work, and I am satisfied based on the additional results provided (which should be included in the updated manuscript) that the approach presents an improvement over existing methods. As such, I increase my rating to accept.