\section{Introduction}
In neuroscience, statistical modeling of the influence of external covariates and self-history to neuronal activity provides an avenue to study the neural code used by the neural system for processing sensory stimuli, cognitive computation, and motor control.
For neurons where the observations are sequence of all-or-none action potentials, or spike trains, the field has widely adapted the use of an autoregressive point process model, often referred to as the generalized linear model (GLM) or the nonlinear Hawkes process~\cite{Truccolo2005,pillow_spatio-temporal_2008,Truccolo2016,Hart2020}.

To incorporate the time-varying modulation controlled by the covariates, GLMs and GLM-like autoregressive point process models are parametrized with (time) convolutional filters~\cite{tripathy_intermediate_2013,latimer_inferring_2019}.
Given time-varying covariates $x_i(t)$, the neural spiking activity $y(t)$ is fully specified as a point process with the conditional intensity function~\cite{Truccolo2005}:
\begin{equation}
    \lambda(t;\mathscr{H}_t) = G
        \left(
            \sum_i (f_i \ast x_i)(t) + d
        \right)  %
    \label{eqn:conditional_intensity}
\end{equation}
where $f \ast x = \int f(\tau) x(t-\tau) \mathrm{d}\tau$ denotes convolution, $f_i(\tau): \field{R} \to \field{R}$ are temporal filters, $d \in \field{R}$ is the bias term that defines a baseline firing rate, $G(\cdot): \field{R} \to [0, \infty)$ is a pointwise nonlinearity, and $\mathscr{H}_t$ denotes the filtration on the past~\cite{Daley1988}.
The conditional intensity function $\lambda(t;\mathscr{H}_t)$ describes the instantaneous firing rate at time $t$ given its history.
In practice, the filters are further parameterized by a finite set of parameters:
\begin{align}\label{eq:GLMwBasis}
    f_i(\tau) = \sum_{j=1}^{N_b} \beta_{i,j} b_j(\tau)
\end{align}
where $N_b$ number of fixed basis functions $b_j: \field{R} \to \field{R}$ are weighted by coefficients $\beta_{i,j} \in \field{R}$.
This allows (1) efficient model fitting, (2) smoothness of $f_i$ controlled by the smoothness of $\{b_j\}$, and (3) temporal support of $f_i$ to be prespecified.
For example, in neural data analysis, the raised cosine basis set is used for covariates and its temporally log-scaled version in combination with a boxcar basis, to capture fast transition to refractory period, is used for spiking history.~\cite{pillow_spatio-temporal_2008}.
In practice, the span of the filter can last for a few hundred milliseconds in the early sensory neurons to a few seconds in higher order cortical areas.
However, the burden for the neurostatistician is to tune the \textit{hyperparameters}, e.g., the smoothness, the temporal span, and the number of basis functions.
This creates a large landscape over the space of possible models and requires extensive model selection to choose the best performing one.
These hyperparameter choices are among the top frequently asked question by the new users of GLM framework for neural data analysis, and improper choice can lead to misleading scientific conclusions or poor inference as seen in Fig.~\ref{fig:filter_space}

\begin{figure}[t!b]
    \centering
    \includegraphics[width=1.0\textwidth]{fig1_supp_arxiv.pdf}
    \caption{
        Mishaps that can occur using a fixed set of bases. \textbf{A)} For a fixed set of bases only a subset of possible filters are actually realizable.  \textbf{B)} Too many bases used resulting in difficulty capturing smooth properties seen in the ground truth.  \textbf{C)} The temporal span was under-specified and so the resulting inference is not able to capture the full extent of the filter.  \textbf{D)} A linear combination of an inappropriate set of bases has trouble inferring characteristics outside its realm. Shown below each example is the set of bases used.}
    \label{fig:filter_space}
\end{figure}

Here, we propose an automated method that alleviates the GLM user from manual tuning of parameters. We make use of the Gaussian Process (GP) framework to impose flexible and expressive priors over the set of filters specified in the model. Being nonparametric in nature this approach requires little prior parameterization and the ability to update hyperparameters during optimization endows us with a rich space of filters that can be captured.

\section{Method}
\subsection{Gaussian processes and sparse approximation}
A Gaussian process is a time continuous stochastic process of which every finite collection of random variables follows a multivariate Gaussian distribution~\cite{Rasmussen2005}.
GPs are widely-used for modeling functions with a probabilistic view that quantifies uncertainty and suits Bayesian inference.
We model each temporal filter $f(\tau)$ over time lag $\tau \in \mathfrak{T} \subset \mathbb{R}$ with a GP:
\begin{equation}
    f(\tau) \sim \mathcal{GP}(m(\tau), k(\tau, \tau'))
\end{equation}
which is fully specified by the mean function $m(\tau)$ and covariance function $k(\tau, \tau')$ (a.k.a. kernel).

For a finite set of points (lag times) $\boldsymbol{\tau} = \{\tau_n\}_1^N$, the corresponding values $\vf = \{f(\tau_n)\}_1^N$ are Gaussian distributed s.t.
\begin{equation}
    \vf \sim \mathcal{N}(\vm_{\boldsymbol{\tau}}, \vK_{\boldsymbol{\tau}})
\end{equation}
where $\vm_{\boldsymbol{\tau}} = [m(\tau_n)]_{n=1,\ldots,N}^\top$ and $\vK_{\boldsymbol{\tau}} = [k(\tau, \tau')]_{\tau, \tau' \in \boldsymbol{\tau}}$.
One can interpolate or extrapolate the value of the filter at any new time $\tau^*$ by
\begin{equation}
    p(f^* \mid \vf) = \mathcal{N}(m(\tau^*) + \vk_{*\boldsymbol{\tau}}^\top \vK_{\boldsymbol{\tau}}^{-1}\vm_{\boldsymbol{\tau}}, k_{*} - \vk_{*\boldsymbol{\tau}}^\top \vK_{\boldsymbol{\tau}}^{-1}  \vk_{*\boldsymbol{\tau}})
\end{equation}
where $f^* \triangleq f(\tau^*)$, $\vk_{*\boldsymbol{\tau}} = [k(\tau^*, \tau)]_{\tau \in \boldsymbol{\tau}}^\top$ and $k_{*} = k(\tau^*, \tau^*)$.

Though Gaussian processes have nice properties, the $\mathcal{O}(N^3)$ computational complexity of training obstructs scalable and efficient application to even moderately sized problems.
To alleviate this obstacle, we use the sparse approximation to Gaussian processes~\cite{snelson2006sparse, titsias2009variational}.
The sparse GP approach introduces a set of inducing variables $\vu = \{f(z_m)\}_{m=1}^M$ that are the filter evaluated at so-called \textit{inducing points} $\{z_m\}_{m=1}^M \subset \mathfrak{T}$.
Moreover these inducing variables are assumed to be a sufficient statistic for $f^*$ s.t.
\begin{equation}
    p(f^* \mid \vf, \vu) =  p(f^* \mid \vu).
\end{equation}
The sparse approximation provides a finite small set of support and reduces the complexity of posterior inference to $\mathcal{O}(N M^2)$ where $M \ll N$.

\subsection{Discretization}

With the conditional intensity as defined in \eqref{eqn:conditional_intensity}, given the spike times $\{t_n\}_1^N$, the point process log-likelihood of the interval $[0,T]$ is
\begin{equation}
    \log p(\{t_n\} \mid \{x_i(t)\}, \{f_i(\tau)\}, \theta) = \sum_n\log\lambda(t_n) - \int_0^T \! \lambda(t) \,\mathrm{d} t    \label{eqn:pp_ll}
\end{equation}
where $\theta$ represents the parameters~\cite{Daley1988}.
In order to evaluate~\eqref{eqn:pp_ll}, we have to calculate the integral $\int_0^T \! \lambda(t) \,\mathrm{d} t$. This is however often computationally costly. Moreover, many covariates such as stimuli are given or can be discretized piecewise constant along with the filter.
Dividing the interval $[0,T]$ into $K$ evenly spaced bins of width $\Delta$, and defining $\lambda_k = \lambda(k\Delta)$ for $k=0,\ldots, K-1$ the log-likelihood is commonly approximated as~\cite{Truccolo2005},
\begin{equation}
    \log p(\{t_n\} \mid \{x_i(t)\}, \{f_i(\tau)\}, \theta) = \sum_n\log\lambda(t_n) - \sum_k \Delta\lambda_k
\end{equation}

Meanwhile the discretization also involves discrete convolution of which requires the value of filter to be evaluated for each time bin.
Limited by finite amount of data and computational resources, the filters need to be truncated into a finite window $[\tau^{(L)}_i, \tau^{(U)}_i]$ at the tails.
Supposing the evaluation of the $i$-th filter to be $\tilde{\vf}_i = \{f_i(\tau)\}_{\tau \in \{\tau^{(L)}_i, \tau^{(L)}_i + \Delta, \ldots, \tau^{(U)}_i\}}$. Then the discrete convolutions can be written as matrix vector multiplications s.t. $\vX_i \tilde{\vf}_i$ where $\vX_i$ is the design matrix $\vX_{i,nj} = x_i(n\Delta - (\tau_i^{(L)} + j\Delta))$. The log-likelihood can then be written compactly as
\begin{equation}
    \log p(\vy \mid \{\vX_i\}, \{\tilde{\vf}_i\}, \theta) = \vy^T (\sum_i \vX_i \tilde{\vf}_i + \mathbbm{1} d) - \Delta \mathbbm{1}^T \exp(\sum_i \vX_i \tilde{\vf}_i + \mathbbm{1}d)
\end{equation}
where $\mathbbm{1}$ is the all-ones vector and $\vy \in \mathbb{R^K}$ is the binary vector of discretized spikes.

\subsection{GP mean and covariance kernel selection}
The kernel of the GP prior has a strong influence on the characteristics of the posterior distribution.
Biophysical characteristics of the neural system is such that the filters for external covariates are (1) local in time lag, (2) temporally smooth, and (3) have fading memory and typically decay to zero exponentially fast.
The decaying squared exponential (DSE) kernel~\cite{tobar_learning_2015} fulfills these conditions. Moreover, to capture filters that are not centered around zero we introduce an offset $\beta$ to the DSE kernel as follows:
\begin{equation}
    k_{\text{DSE}}(\tau, \tau') = \sigma^2 e^{-\alpha(\tau - \beta)^2 -\nu(\tau - \tau')^2 -\alpha(\tau' - \beta)^2}
\end{equation}
where $\alpha$ controls the exponential decay rate, $\nu$ controls the spectral content, and $\sigma^2$ controls the finite power.

In many cases a null mean is used to specify the GP prior as this does not constrain the posterior mean.  However, it can be beneficial to incorporate prior knowledge and bias the model toward an asymptotic regime by introducing a nontrivial parameterization of the mean.  For example, the refractory period exhibited by neurons makes it very unlikely that spikes should occur within close proximity of each other.  This manifests itself in a history filter as an exponential climb from extremely negative values to zero over a short time succeeding a spike. ~\cite{10.1162/neco_a_01021}. Since this is well known we can incorporate our neuroscientific knowledge and parameterize spike history filters as such

\begin{equation}
    f(\tau) \sim \mathcal{GP}(a e^{-\frac{\tau}{b}}, k_{\text{DSE}}).
\end{equation}
which allows for the resulting inference to more easily capture the characteristic refractory period.  Indeed, this serves a similar function to the use of a 'boxcar' basis under the alternative approach.  Furthermore, without this specification the resulting inference may gravitate towards a small length scale parameter to adequately capture any transient sharp rise and thus impede inference about smooth long timescale varying characteristics.

\subsection{Variational inference}
In Bayesian inference all of the uncertainty over a filter $f(\tau)$ once the data $\vy$ has been observed is embedded into its' posterior distribution
\begin{equation}
    p(f(\tau) \mid \vy) = \frac{p(\vy \mid f(\tau)) p(f(\tau))}{p(\vy)}
\end{equation}

In GLMs, unfortunately, the evidence, $p(\vy)$, and thus the posterior $p(f(\tau) \mid \vy)$ are intractable for GP priors. Sampling methods are able to tackle the intractability. However they require intensive computation and do not produce analytical solutions.
Therefore we seek approximate solutions.
More recent works have used sparse variational approaches which take care of both the intractable posterior as well as the computational complexity of full GP~\cite{lloyd_variational_2015-1, hensman_mcmc_2015}.  Variational inference~\cite{Bishop_PRML} usually assumes a parametric form of an approximation to the posterior distribution and cast the inference problem to optimization.

Continuing with our formulation where $\vu$ are the inducing variables of the filter $f(\tau)$, a lower bound on the log marginal likelihood $\log p(\vy)$ can be derived~\cite{lloyd_variational_2015-1, dezfouli_scalable_2015} yielding the variational objective
\begin{equation}\label{eq:elbo_u}
    \mathcal{L} = \mathbb{E}_{p(f(\tau) \mid \vu) q(\vu)} \log p(\mathbf{y} \mid f(\tau)) - \KL(q(\vu) \,\|\, p(\vu))
\end{equation}
which is often referred to as the evidence lower bound (ELBO).
The $q(\vu)$ closest in KL divergence to the true posterior can then be obtained by maximizing the ELBO.

This lower bound can be evaluated analytically if $q(\vu)$ is selected to be a parameterized Gaussian of the form $\mathcal{N}(\vu \mid \vm, \vS)$ and the joint approximate distribution $q(f(\tau), \vu)$ follows
\begin{equation}
    q(f(\tau), \vu) = p(f(\tau) \mid \vu) q(\vu),
\end{equation}

The resulting marginal $q(f(\tau)) = \int q(f(\tau), \vu) \,\mathrm{d} \vu$ is now an approximation to the true posterior and is also a Gaussian process with mean and covariance given by
\begin{align}
    \mu(\tau) &= K_u (\tau)^T K_{uu}^{-1} \vm \\
    \Sigma(\tau, \tau') &= K(\tau, \tau') - K_u(\tau)^T K_{uu}^{-1 }K_u(\tau) + K_u(\tau)^T K_{uu}^{-1}\vS K_{uu}^{-1} K_u(\tau).
\end{align}

Letting $\vu_i = \{f_i(z_m^i)\}_{m=1}^M, \, z_m^i \in \mathfrak{T}_i$ be the inducing variables of filter $f_i(\tau)$. Eq.~\eqref{eq:elbo_u} can be extended to
\begin{equation}
\begin{split}
    \mathcal{L} = \vy^T \Bigg(\sum_i \vX_i \mu_i & + \mathbbm{1}d \Bigg)
    - \Delta \mathbbm{1}^T \exp\left(d + \sum_i \vX_i \mu_i
    + \frac{1}{2}\mathrm{diag}(\vX_i \Sigma_i \vX_i^T) \right) \\
    &+ \sum_i \KL\left(q(\vu_i) \| p(\vu_i)\right)
\end{split}
\end{equation}
where \[\mu_i = K_{X_i, X_i} K_{u_i, u_i}^{-1} m_i\] and
\[\Sigma_i = K_{X_i,X_i} - K_{X_i,u_i}K_{u_i,u_i}^{-1}K_{u_i,X_i} + K_{X_i,u_i}K_{u_i,u_i}^{-1} S_i K_{u_i,u_i}^{-1}K_{u_i, X_i}\]
are respectively the mean and covariance of the variational approximation $q(\vu_i)$. Now, full inference can be performed as the  gradients of the ELBO w.r.t. the variational parameters $\{\vz^i, \vm_i, \vS_i\}$ and hyperparameters $\{\boldsymbol{\theta}_i\}$ can be calculated in closed form

\subsection{Optimization}
As the variational objective is differentiable w.r.t. hyperparameters, variational parameters, and model parameters, we can maximize the ELBO
without resorting to sampling schemes or approximations. We elect to use a full rank parameterization of the covariance of each variational Gaussian, $\vS_i$, by optimizing its Cholesky factor $\vL_i$, such that $\vS_i = \vL_i \vL_i^T$ remains in the positive semi-definite cone.
While the optimization can be performed jointly we found a coordinate ascent approach more practical and consistent.  This involves partitioning the parameters of the model into two groups -- the first group $\vV_i$, containing the variational parameters, and a second group $\vH_i$, containing the kernel hyperparameters.

In the optimization process we also incorporate updates of the convolutional design matrix for each of the covariates.  Each design matrix, $\vX_i$, and the associated set of lags over the span  $[\tau^{(L)}_i, \tau^{(U)}_i]$ are byproducts of having to choose a finite window over the lag space to realize the convolutional interaction under our model.  When the decision to restrict this window to a certain temporal region is made it is possible that the temporal extent will be underspecified.  Depending on the extent of the underspecification the model will lose the capability to explain phenomena of the firing rate that occur on long time scales or large latencies with respect to a particular covariate.  While the GP prior on the filters provides us a posterior that is well defined over the entire lag space, optimization of the model parameters i.e. hyperparameters and variational parameters, require evaluation of the expected log-likelihood which in turn is directly affected by the temporal extent of the filter used to perform the convolutions.

\begin{algorithm}[H]
\newcommand\mycommfont[1]{\footnotesize\ttfamily\textcolor{black}{#1}}
\SetCommentSty{mycommfont}
\SetKwInput{kwReturn}{Return}
\SetKwBlock{kwRepeatUntil}{Repeat Until Convergence}{end}
\SetAlgoLined
initialize parameters: \hspace{1pt} $\vV^{(0)}_k = \{\vm_k^{(0)}, \vS_k^{(0)},
\vz_k^{(0)}\}$,\hspace{2pt} $\vH^{(0)}_k = \{\boldsymbol{\theta}_k^{(0)}, d_k^{(0)}\}$, \hspace{2pt} $k=1,\ldots,K$ \\
initialize design matrices: $\{X_i\}$ \\
Set $j=1$ \\

\kwRepeatUntil{
\tcp{for each set of variational/inducing variables and hyperparameters}
    \For{$\vP^{(j-1)} \in \{\vV_k^{(j-1)}, \vH_k^{(j-1)}\}$} {
    \tcp{for each covariate}
        \For{$i=1,\ldots,K$}{
            Set $\vP^{(j)} = \underset{\vP}{\argmax} \, \mathcal{L}$
        }
    }
    \tcp{update the design matrix for each covariate}
    \For{$i=1,\ldots,K$}{
    Set $\tau_i^{(L)} = \floor{\underset{m}\min \, {z^{(j)}_{i,m}} / \Delta} \Delta$ \\
    Set $\tau_i^{(U)} = \ceil{\underset{m}\max \, {z^{(j)}_{i,m}} / \Delta} \Delta$ \\
    Update $\mathfrak{T}_i$, $X_i$
    }
    Set $j = j + 1$
}
\caption{NPGLM Sequential Training Algorithm}\label{algo}
\end{algorithm}

The shortcomings of having to choose this finite duration time window can be amended during the optimization procedure thanks to the inferential power granted by having used the GP framework.  In a binned time description of the model this only amounts to evaluating the variational posterior at additional time bins.  As is natural though, we would like to be greedy and use the smallest amount of evaluations we can in order to realize the convolution operation.  In that light, we propose viewing the optimized inducing points as ``guides'' to regions of high posterior density that contribute significantly to the resulting convolution.  More specifically, if the optimized inducing points for a particular covariate's filter carve out a window $[\tau^{(L)}_i, \tau^{(U)}_i]$,  then we choose to evaluate our convolutions with those as the evaluation points of the filter's posterior.  This affords us to be greedy because under-specification of this evaluation window can be amended in an iterative procedure, and more than that we avoid evaluations at points in time that have no relevance and only serve to harm the optimization procedure. Algorithm~\ref{algo} summarizes the optimization procedure.  In practice we had the most success with this style of optimization where one set of parameters ($\vV_i$ or $\vH_i$) is optimized for one covariate at a time, although we note that joint optimization of all covariates and all parameters at once worked but was less robust.

\subsection{Adaptive Design Matrix Intuition}
Quickly, we would like to talk through the process of updating the design matrix, why it is necessary, and clarify further how it is done.  Take for example, inference about the underlying filter in Fig.~\ref{fig:toy_evolution}A.  Under the alternative approach of picking temporally fixed basis functions if you were to select a window of 50 $ms$ then inference about the filter would turn out poor.  Consider the approach described and an initial set of inducing points spaced uniformly from 0 $ms$ to 50 $ms$.  Now, to properly evaluate the ELBO we must convolve the discretized stimuli $x_i(n)$ with the posterior mean of the filter $\bar{\vf}$ for all $K$ time bins of observations $\vy$.  Since $\vf$ has support over all of $\mathbb{R}$ the only way to realize the convolution is to window $\vf$ which will make it 'appear' to the likelihood that $\vf$ is null outside this range. Concretely, on the first round of optimization you would have the approximation:

\begin{align*}
    \text{E}\left(\vy^T(\vx * \vf + \mathbbm{1}d) - \Delta\mathbbm{1}^T\exp(d+\vy * \vf)\right) \approx \vy^T \Bigg(\vX \mu & + \mathbbm{1}d \Bigg)
    - \Delta \mathbbm{1}^T \exp\left(d +  \vX \mu
    + \frac{1}{2}\mathrm{diag}(\vX \Sigma \vX^T) \right)
\end{align*}

\begin{figure}[t!b]
    \centering
    \includegraphics[width=1.0\textwidth]{neurips_2020_inducing_plot.pdf}
    \caption{
        Example showing the effects of having to choose a finite window over the filter to perform convolutions.  The finite window impedes inference due to a misinformed likelihood, however, using the span dictated by the inducing points we can realize a mechanism that adaptively windows more of the relevant space.  In this way computation time can be saved by choosing a smaller initial window and allowing it to adapt over time.}
    \label{fig:inducing_intuition}
\end{figure}

where $X_i \in \mathbb{R}^{K \times 50}$ since we chose 50 $ms$ as the span of the inducing points and elected to make the cutoff based on their span.  Thus, the approximation to the likelihood will be off the mark.  What can save us is the fact that the prior will be in disagreement and so some inducing points will be 'pulled' beyond the imposed cutoff.  If this is the case the window can be increased and consequently $X$ will be expanded. Note how in Fig~\ref{fig:inducing_intuition} the range of the inducing points dictate the effective filter seen by the likelihood --- as they extend further the truncated convolution more closely matches that of one performed using the full GP posterior.

\begin{figure}[t!b]
    \centering
    \includegraphics[width=0.9\textwidth]{figure1_mp}
    \caption{
        \textbf{NPGLM validation with known ground truth.}
        (\textbf{A})
        NPGLM recovers the span of the filter even with a bad initialization.
        NPGLM with a stimulus filter compromising 20 inducing points is initialized to have a temporal span of 50~ms.
        We plot the posterior at different optimization iterations.
        Note how not the entire span of the filter need be used when evaluating the convolutions in order to achieve a satisfactory approximation.
        Shading indicates the 95\% credible interval of the filter.
        (\textbf{B})
        NPGLM is robust to small sample size.
        Two Bayesian methods, NPGLM and GLM-ARD are trained on spike trains of different duration.
        Sample size from left to right: 500 (7 spikes), 5,000 (72 spikes), 15,000 (215 spikes) time bins.
        In left we see that posterior fit offered by GLM-ARD has little resemblance to the ground truth, whereas NPGLM appears more in line with the ground truth and achieves higher normalized log-likelihood.
        In middle and right the fit models both more closely resemble the ground truth with NPGLM achieving higher log-likelihood on the hold out set.
    }
    \label{fig:toy_evolution}
\end{figure}

\section{Experiments}
To examine the methodology described in the paper we examine its performance on synthetic data as well as real neural recordings. In all experiments we measure performance relative to standard GLM approaches where the filters are parameterized by a set of fixed basis functions.  One GLM we compare to places an automatic relevance determination (ARD)~\cite{qi2004predictive} prior over the weights of the basis functions -- we refer to this as GLM-ARD. For the other GLM no prior is set over the weights and all inference is performed using the MLE -- we refer to this one as GLM-MLE.  The methodology described in the paper is referred to as Non-Parametric GLM (NPGLM).

\subsection{Synthetic experiment 1}
We create a stimulus filter that is a sum of two Gaussians, specifically
\begin{equation}
    s(t)=3\,\mathcal{N}(t \mid 75, 12.5^2) - 1.8\,\mathcal{N}(t \mid 62.5, 12.5^2)
\end{equation}
and a history filter that is supposed to mimic a short refractory period followed by a transient self-excitation.

\begin{figure}[t!b]
    \centering
    \includegraphics[width=0.9\textwidth]{rebuttal_fig1_arxiv.pdf}
    \caption{
        Inference over a synthetic neuron with several tuned stimuli responses, a history filter, and positive coupling from two other neurons, A and B.  NPGLM is able to accurately capture the structure of the ground truth filters.  In many cases we see the basis function approach has significant non-zero values after the ground truth has subsided to zero.  Additionally, GLM-ARD is more confident even when inference does not capture the ground truth sufficiently.  Shading represents the 95\% credible region for both approaches.  NPGLM outperforms GLM-ARD  when examing log-likelihood on a held out test set (History filter panel top right).  }
    \label{fig:synthetic_expt_2}
\end{figure}
We use this toy example to quantify A) how well the adaptive truncation window mechanism functions and B) how inference scales with the amount of relevant data. Figure~\ref{fig:toy_evolution}B shows the resulting stimulus filters after fitting NPGLM and GLM-ARD for various length spike trains.  When there are only 7 spikes in the entire sequence the posterior under GLM-ARD does not bear much resemblance to the ground truth used for generating the data.  NPGLM also has trouble capturing the ground truth but qualitatively speaking the fit appears much closer. For the two other examples where 72 and 215 spikes were used when fitting the models both capture the ground truth reasonably well. However the wider tails present in the GLM-ARD fit do end up adversely affecting its evaluation on the test set.

To illustrate the ability of NPGLM to adaptively learn the temporal span of the filter in Figure~\ref{fig:toy_evolution} we initialize the stimulus filter so that the most extreme inducing point is right on the cusp of an influential part of the filter.  Specifically 10 inducing points are initialized uniformly from 0~ms to 50~ms. As described in Algorithm~\ref{algo} migration of the inducing points towards extreme locations in space leads to adaption of the design matrix in turn leading to a larger effective window to evaluate the convolutions.  Following this strategy, it is not long before the inducing points have migrated far enough that the posterior bears a strong resemblance to the ground truth.

\subsection{Toy Experiment 2}
Now, we examine a more complicated example.  This time three stimuli and two coupling filters are used to generate the synthetic data.  Each of the stimuli are associated with a distinct neural response dictated by the appropriate filter.  Further, two neurons 'A' and 'B' are coupled positively at different latencies to our 'reference' neuron.  NPGLM is initialized with inducing points spaced 15 $ms$ apart for the stimuli filters, 3 $ms$ apart for the history filter, and as before the DSE kernel is chosen for all GPs involved.  GLM-ARD is fit using the raised cosine basis and the temporal span is chosen to adequately cover each of the filters.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.8\textwidth]{rebuttal_fig2_arxiv.pdf}
    \caption{Run times as a function of the number of inducing points.  NPGLM is robust against the number of inducing points chosen. We can see using additional inducing points does not increase run times drastically due to the round-robin style of optimization. In contrast quality of inference when using a set of bases functions will depend significantly on their parameters.}
    \label{fig:runtimes}
\end{figure}

In Fig.~\ref{fig:synthetic_expt_2}  we can see for each of the synthetic filters that NPGLM is able to adequately capture their identifying features.  The alternative, GLM-ARD is able to as well although it appears to struggle maintaining smoothness in the tails of the coupling filters.  Additionally, we can note some spurious sharp features inferred in some of the stimuli filters where they were predominantly smooth.

We then used this same example to quantify the run time of NPGLM as a function of the number of inducing points.  For each of the stimuli filters inducing points varied from spacings of 5 $ms$ apart (more inducing points) to spacings of 60 $ms$ apart (less inducing points),  for example, in the extreme case of 60 $ms$ spacing the filter for Stim3 uses only 7 inducing points. Examining Fig.~\ref{fig:runtimes} run times are longer than GLM-ARD, however, this is to be expected and it is clear that full optimization is not required for a satisfactory level of inference.  Furthermore, even though NPGLM is not convex in it's parameters it achieves consistent, better inference, and is robust to the number of inducing points chosen.

\begin{figure}[tb]
    \includegraphics[width=\textwidth]{summary_data_expt3_fig2_ink}
    \caption{
        \textbf{NPGLM inference on six individual ON/OFF retinal ganglion cells is consistently superior to baseline models.}
        (\textbf{A}) Filters inferred by NPGLM using data from~\cite{Pillow2005}.
        For the most part the filters fitted to the NPGLM model strongly resemble the GLM-MLE and GLM-ARD fits.
        For certain neurons we see that different features may make themselves more well known i.e. neuron 4 and neuron 2.
        See Fig.~\ref{fig:RGC05:summary} for further details.
        (\textbf{B}) GLM-ARD fits.
        (\textbf{C}) GLM-MLE fits.
        (right) normalized test log-likelihood comparison shows that NPGLM consistently outperforms.
    }
    \label{fig:RGC05:all}
\end{figure}

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{summary_data_expt3_fig4_mp}
    \caption{
        \textbf{Spike trains generated from NPGLM outperforms baseline.}
        Two example retinal ganglion cells (\textbf{A} and \textbf{B}) fit with NPGLM and GLM with basis functions (ARD w/b).
        Spike trains are simulated using ancestral sampling.
        \textbf{(1)} Raster plot of data and model simulation (top: NPGLM; middle: data; bottom: ARD). The identical random stimulus was repeated for multiple trials.
        \textbf{(2)} Average firing rates estimated from data and simulations.
        \textbf{(3)} Empirical cumulative distribution of inter-spike intervals.
        \textbf{(4)} Goodness-of-fit measures. From left to right: normalized log-likelihood (nats/bin) for training, test and repeat datasets, and $R^2$ of firing rate for the repeat dataset.
    }
    \label{fig:RGC05:summary}
\end{figure}
\subsection{Retinal ganglion data}
We examine \textit{in vitro} multielectrode extracellular recordings from the retina of macaque monkeys~\cite{Pillow2005}. The monkeys were presented a spatially uniform cathode ray tube display refreshing at 120~Hz. Each frame took on one of two intensity values.  Previous analyses served to quantify the statistical differences in behavior between ``ON'' and ``OFF'' -- so the data set offers an established testbed for model comparison.  For this example NPGLM, GLM-ARD, and GLM-MLE are all fit to the recorded data and the bases used for both filters in the latter approaches is a series of raised cosine bumps scaled logarithmically in time.  For the stimuli filter 15 such bases spanning 300~ms are used and for the history filter 15 such basis spanning 100~ms are used.  NPGLM uses 15 inducing points spaced uniformly in time for the stimulus filter and 25 inducing points spaced logarithmically in time for the history filter.

For all six neurons the fit under each model is shown in Figure~\ref{fig:RGC05:all}.  The normalized log-likelihood for a novel repeat stimulus was also calculated to quantiatively compare how well each model performed.  For all of these cases except one we saw that the NPGLM approach performed better.  In the case of Neuron 2 it appears that NPGLM was unable to capture a sharper rise in the history filter occurring at 10~ms.  In practice the history filter is the more difficult of the two to capture because of its non-smooth behavior.

In Figure~\ref{fig:RGC05:summary} we examine one ``ON'' cell and one ``OFF'' cell in more depth.  For Neuron 2 we see that the NPGLM approach performs slightly better on the set of data used to fit the model but noticeably better when it comes to generalizing to the repeat stimulus.  This highlights the ability of such a nonparametric approach to generalize well -- even if the data is sparse in the number of spikes as is the case with Neuron 2.  Contrast this to Neuron 4 where while still outperforming the standard approach on the test the difference is not as pronounced.

\section{Discussion}
In this paper we proposed NPGLM, a nonparametric scheme for making automated inferences for temporal filters that characterizes the influence of external event and stimuli as well as self-generated cognitive decisions and behavior to neural spikes.
NPGLM infers time convolutional filters in autoregressive point process models by utilizing sparse Gaussian process prior and variational inference in a computationally efficient manner.
We verified that NPGLM inference is robust against variations in the sparsity of the GP approximation, and it produces reasonable posterior inference even in the small data regime.
Contrast to previous approaches that model the filters with basis functions, we sidestep the burdensome task of basis function selection and rigorously defining their fixed temporal properties.
We verified the competency of our method on retinal ganglion cell recordings, where it was able to infer filters consistent with previous approaches using basis functions.

While we used binned data in this study, we posit that using spike times at a higher time resolution may reveal structures or features not readily seen using other methods thanks to the infinite time resolution of GPs.
Further extensions can include tractable methods of inference for determining a more optimal convolution window, pruning inducing points that are not necessary, and adding more inducing points when appropriate.

\section*{Acknowledgement}
This work is supported by NSF CAREER IIS-1845836, IIS-1734910, and NIH UF1-NS115779.

\bibliographystyle{hunsrt}


\section{Introduction}

Symbolic dynamics, number theory and theoretical computer science share many connections, but arguably none so pervasive as their focus on understanding patterns, and specifically their desire to understand how relatively simple procedures---e.g., doubling, recurrence, iteration---can take the most mundane objects and produce rich structures with long range order. An interval becomes the Cantor set, two tiles arrange to give Penrose's tiling and a pair of simple functions lead to a strange attractor. While several notions and concepts have been introduced to describe the structure and symmetry of these limiting objects, the most successful of these has been the harmonic analysis of measures and the related theory of (fractal) dimension---for instance, via thermodynamic formalism; see Baake et al.~\cite{BGKS2019} and Kesseb\"ohmer and Strattmann \cite{KS2007}.

A central object in this context, ubiquitous to each of the above-mentioned areas, are substitutions systems. The long range order of substitution systems has been well studied; there is an abundance of literature on this still very active area of research going back to the seminal works of Wiener \cite{W1927} and Mahler \cite{M1927}. The monographs by Queff\'elec \cite{Qbook} and Baake and Grimm \cite{BGbook} contain details about both the tiling and symbolic pictures of these systems as well as the associated diffraction and spectral measures---the modern-classical means of examining the long range order of these systems. A canonical example is the Thue--Morse substitution on two letters given by
\begin{equation*}
   \varrho_{_{\rm TM}}: \, \begin{cases} a \mapsto ab \\
     b\mapsto ba\, . \end{cases}
\end{equation*} Applying $\varrho_{_{\rm TM}}$ iteratively to $a$, one obtains the fixed point $$\varrho_{_{\rm TM}}^{\infty}(a)=abbabaabbaababbabaababbaabbabaab\cdots.$$
Researchers in different areas tend to associate different objects to such substitutions. In symbolic dynamics, one is interested in the orbit closure of $\varrho_{_{\rm TM}}^{\infty}(a)$ under the shift map and in the induced dynamical system, the Thue--Morse dynamical system. In theoretical computer science, combinatorics on words and number theory, one is interested in---among other things---subword complexity and related questions of transcendence and Diophantine approximation. In our context, we choose values for $a$ and $b$, say $a=-b=1$, and define the Thue--Morse sequence ${t}$ as the $n$th term of the fixed point $\varrho_{_{\rm TM}}^\infty(1)$. The substitution then implies that both $t(2n)=t(n)$ and $t(2n+1)=-t(n).$ Analogous to the Thue--Morse word, here we consider the $2$-kernel of ${t}$, that is, the set of subsequences ${\rm ker}_2 ({t}):=\{(t(2^\ell n+r))_{n\geqslant 0}: \ell\geqslant 0, 0\leqslant r<2^\ell\}.$ In this particular case, the set ${\rm ker}_2 ({t})$ is finite; it contains only the sequences $(t(n))_{n\geqslant 0}$ and $(-t(n))_{n\geqslant 0}$.

In general, for an integer sequence ${f}$ and an integer $k\geqslant 2$, we consider the $k$-kernel of ${f}$, that is, the set of arithmetic subsequences $${\rm ker}_k ({f}):=\left\{(f(k^\ell n+r))_{n\geqslant 0}: \ell\geqslant 0, 0\leqslant r<k^\ell\right\}.$$ Like the example of the Thue--Morse sequence, if ${\rm ker}_k ({f})$ is a finite set, there is a related substitution system, and the sequence is called {\em $k$-automatic}---these sequences are codings of sequences that are the output of a deterministic finite automaton, which takes in the base-$k$ expansion of $n$ and outputs $f(n)$. We are interested in the generalisation wherein ${\rm ker}_k ({f})$ generates a finitely generated $\mathbb{Q}$-vector space $$\mathcal{V}_k(f):=\langle{\rm ker}_k ({f})\rangle_\mathbb{Q}.$$ In this case, the sequence $f$ is called  {\em $k$-regular} (or just {\em regular} when $k$ is understood). These sequences were introduced by Allouche and Shallit \cite{AS1992} in the early nineties as a natural generalisation of automatic sequences---a regular sequence that takes only finitely many values is automatic. The sequence $f$ and the space $\mathcal{V}_k(f)$ are the central foci of this paper. In particular, {\em we study the long range order of both $f$ and $\mathcal{V}_k(f)$ by treating them as dynamical objects and constructing a systematic measure-theoretic framework surrounding them, as regular sequences.}

In comparison to automatic sequences, the study of the long range order of unbounded regular sequences $f$, and so also the related spaces $\mathcal{V}_k(f)$, is not so straight-forward; neither diffraction nor spectral measures can be associated to them in a natural way.

As a first step of addressing the long range order of such objects, Baake and Coons~\cite{BC2018} introduced a natural probability measure associated to Stern's diatomic sequence, the canonical example of an unbounded regular sequence. Stern's diatomic sequence ${s}$ is given by $s(0)=0$, $s(1)=1$, and for $n\geqslant 1$ by the recurrences $s(2n)=s(n)$ and $s(2n+1)=s(n)+s(n+1)$. Its $2$-kernel ${\rm ker}_2 ({s})$ is infinite, but induces a finitely generated vector space; $$V_2(s)=\langle{\rm ker}_2 ({s})\rangle_{\mathbb{Q}}=\langle\{(s(n))_{n\geqslant 0},(s(n+1))_{n\geqslant 0}\}\rangle_{\mathbb{Q}}.$$ Stern's diatomic sequence has many interesting properties, e.g., the sequence of ratios $(s(n)/s(n+1))_{n\geqslant 0}$ enumerates the non-negative rationals in reduced form and without repeats, and the associated measure is singular continuous. Baake's and Coons's result rests on the fact that the sequence ${s}$ satisfies certain self-similar type properties: it has a fundamental region of recursion---between consecutive powers of two---and the sum of ${s}$ over this fundamental region is linearly recurrent, enabling the use of a volume-averaging process.

In this paper, we provide the complete generalisation of the result of Baake and Coons~\cite{BC2018}. This entails not only taking a $k$-regular sequence $f$ and associating a natural probability measure $\mu_f$ to $\mathcal{V}_k(f)$, but also determining the class of regular sequences $f$ for which this is possible. Indeed, there are regular sequences for which such measures do not exist, and for which the sequence of approximants is not eventually periodic.

The measure $\mu_f$ we construct is analogous to mass distributions supported on attractors of certain iterated function systems---those constructed by repeated subdivision; see Falconer~\cite{Fbook} for related definitions and Coons and Evans~\cite{CEpre} for a family of extended examples related to generalised Cantor sets. Here, rather than iterating intervals under a finite number of maps, the finite approximants are pure point measures on the unit interval, whose corresponding weights possess the recursive structure.
In this way, the existence of a natural measure associated to $f$ and $\mathcal{V}_k(f)$ provides a path to viewing the space $\mathcal{V}_k(f)$ of regular sequences associated to $f$ as a dynamical system and opens the possibility of considering these sequences and spaces as (fractal) geometric structures.

This article is organised as follows. In Section~\ref{sec:notation} we introduce notation and assumptions, and state our main two theorems. Section~\ref{sec:V} contains the proof of our first theorem, the existence of a measure associated to $\mathcal{V}_k(f)$, and Section~\ref{sec:f} contains a technical general result on the existence of a measure associated to $f$ that implies our second theorem. In Section~\ref{nonprim}, we discuss the necessity of our assumptions and provide witnessing examples. Finally, we offer some concluding remarks and open questions in Section~\ref{sec:cr}.

\section{Notation, assumptions and statement of results}\label{sec:notation}

The contents of the following section are technical, but they formalise a natural idea. We start with an integer-valued $k$-regular sequence $f$ and obtain a basis for $\mathcal{V}_k(f)$, the $\mathbb{Q}$-vector space generated by the $k$-kernel of $f$, which is maximal with respect to asymptotic growth over certain fundamental regions. In particular, let $k\geqslant 2$ be an integer, $f$ be a $k$-regular sequence and let the set of integer sequences \begin{equation}\label{basis}\{f=f_1,f_2,\ldots,f_d\}\subseteq {\rm ker}_k ({f})
\end{equation} be a basis for $\mathcal{V}_k(f):=\langle{\rm ker}_k ({f})\rangle_\mathbb{Q}$. Set ${\bf f}(m)=(f_1(m),f_2(m),\ldots,f_d(m))^T$. For each $a\in\{0,\ldots,k-1\}$ let ${\bf B}_a$ be the $d\times d$ integer matrix such that, for all $m\geqslant 0$,
\begin{equation}\label{eq:k kernel}
{\bf f}(km+a)={\bf B}_a\, {\bf f}(m).
\end{equation}
We refer the reader to the seminal paper of Allouche and Shallit \cite{AS1992} and Nishioka's monograph \cite[Ch.~15]{N1996} for details on existence and the finer definitions. Note that there is ${\bf w}\in\mathbb{Z}^{d\times 1}$ such that for each $i\in\{1,\ldots,d\}$ and $n>0$, we have $$f_i(m)=e_i^T {\bf B}_{(m)_k} {\bf w}=e_i^T{\bf B}_{i_0}{\bf B}_{i_1}\cdots {\bf B}_{i_s}{\bf w},$$ where $e_i$ is the $i$th elementary column vector, $(m)_k=i_s\cdots i_1i_0$ is the base-$k$ expansion of $m$ and ${\bf B}_{(m)_k}:={\bf B}_{i_0}{\bf B}_{i_1}\cdots {\bf B}_{i_s}.$ Set $${\bf B}:=\sum_{a=0}^{k-1}{\bf B}_a.$$

As we start to view the vector space $\mathcal{V}_k(f)$ dynamically, we note that the matrix ${\bf B}$ is analogous to the substitution matrix ${\bf M}_\varrho$ of a substitution $\varrho$ on a finite alphabet. In fact, if $\varrho$ is a constant-length substitution of length $k$ on the $d$ letters $0,1,\ldots,d-1$, then the matrices ${\bf B}_a$ are the so-called digit (instruction) matrices and ${\bf M}_\varrho={\bf B}$; compare \cite{BGM2019,Qbook}. When considering the dynamical properties of substitution systems it is common to assume that the substitution is primitive; that is, the non-negative substitution matrix ${\bf M}_\varrho$ is primitive. To continue our analogy with substitutions, we make similar assumptions. At first glance, it is reasonable to restrict ourselves to the assumption that ${\bf B}$ is primitive. But if this is the case for the $k$-regular sequence $f$, then we can consider $f$ as a $k^j$-regular sequence with $j$ being the smallest positive integer for which ${\bf B}^j$ is positive. Thus, in the context of $k$-regular sequences the distinction between primitivity and positivity is somewhat blurred. Note that this is also tacitly done for substitutions, where one normally chooses an appropriate power $j$ such that ${\bf M}_{\varrho}^{j}>0$, or equivalently, $\varrho^{j}(a)$ contains all the letters of the alphabet $\mathcal{A}$, for all $a\in \mathcal{A}$. Hence we make the following definition.

\begin{definition}\label{def:prim} We call a $k$-regular sequence $f$ {\em primitive} provided $f$ takes non-negative integer values, is not eventually zero, each of the $k$ digit matrices ${\bf B}_a$ are non-negative and the matrix ${\bf B}$ is positive.
\end{definition}

\noindent In Section \ref{nonprim}, we discuss the necessity of these assumptions. In particular, we present a sequence whose associated matrix is not primitive and show that our construction does not converge to a measure for this sequence.

We are close to being able to state our result. The final ingredients are the introduction of fundamental regions and the definition of the related pure point measures. To this end, we note that, in analogy to the lengths of the iterates of a substitution applied to an initial seed, a $k$-regular sequence $f$ has a fundamental region, the interval $[k^n,k^{n+1})$. Over these regions, the sum of a regular sequence is linearly recurrent, a property which can be thought of as a sort of self-similarity property for regular sequences. This property provides a structured volume to average by, and using it, we can construct a probability measure. Formally, for each $i\in\{1,\ldots,d\}$, set \begin{equation}\label{eq:Sigmai}\varSigma_i(n):=\sum_{m=k^n}^{k^{n+1}-1}f_i(m)\end{equation} and define
\begin{equation}\label{eq:def-meas}
   \mu^{}_{n,i} \, := \, \frac{1}{\varSigma_i(n)} \sum_{m=0}^{k^{n+1}-k^n - 1}
   f_i(k^n + m) \, \delta^{}_{m / k^n(k-1)},
\end{equation}
where $\delta_x$ denotes the unit Dirac measure at $x$. We can view
$(\mu^{}_{n,i})^{}_{n\in\mathbb{N}_0}$ as a sequence of probability measures on
the $1$-torus, the latter written as $\mathbb{T}=[0,1)$ with addition modulo
$1$. Here, we have simply re-interpreted the (normalised) values of
the sequence $(f_i(m))_{m\geqslant 0}$ between $k^n$ and $k^{n+1}-1$ as the weights of a pure point
probability measure on $\mathbb{T}$ supported on the set $\big\{ \frac{m}{k^n(k-1)} : 0 \leqslant m <
k^n(k-1) \big\}$. Set \begin{equation}\label{eq:mun}\boldsymbol{\mu}_{n}=\boldsymbol{\mu}_{f,n}:=(\mu_{n,1},\ldots,\mu_{n,d})^T.\end{equation}

Our first result, Theorem \ref{main1} below, provides a unique probability measure on $\mathbb{T}$ associated to the space $\mathcal{V}_k(f)$.

\begin{theorem}\label{main1} Let $f$ be a primitive $k$-regular sequence. If $f_1,\ldots,f_d$ form the basis of $\mathcal{V}_k(f)$ associated to ${\bf B}$, then the vectors $\boldsymbol{\mu}_{f,n}$ of pure point measures converges weakly to a vector of probability measures $\boldsymbol{\mu}_f=(\mu_f,\ldots,\mu_f)^T$ on $\mathbb{T}$.
\end{theorem}

Of course, one may also wish to consider the measure associated to $f$ apart from the full considerations surrounding that of the space $\mathcal{V}_k(f)$. As it turns out, with a little more specificity one can prove a stronger result.

\begin{theorem}\label{main2} Let $f$ be a positive integer-valued $k$-regular sequence, such that the spectral radius of ${\bf B}$ is the unique dominant eigenvalue of ${\bf B}$. Then the measure $\mu_f$ exists and is continuous.
\end{theorem}

While, in Theorems \ref{main1} and \ref{main2} we have assumed that $f$ is non-negative, this assumption can be traded for a more general condition in Theorem \ref{main2}. In Section \ref{sec:f}, we prove a slightly stronger result, where we replace this non-negativity assumption with two conditions both having to do with a related dilation equation. The first condition is  a natural non-vanishing condition and the second is that the solution of the associated dilation equation is of bounded variation. This more general situation may not result in a probability measure, but possibly, a signed measure.

\section{A natural probability measure associated to $\mathcal{V}_k(f)$}\label{sec:V}

The aim of this section is to prove Theorem \ref{main1}. To do this, we attempt to mimic the ideas behind the establishment of spectral measures associated to substitution dynamical systems. Fortunately, the generating power series of $k$-regular sequences satisfy functional equations that can be thought of as taking the place of a substitution---the generating power series of the sequences $f_i$ in \eqref{basis} are Mahler functions. That is, for each $i\in\{1,\ldots,d\}$ setting $F_i(z):=\sum_{{m}\geqslant 0}f_i({m}) z^{{m}}$ and ${\bf F}(z)=(F_1(z),\ldots,$ $F_d(z))^T$, we have (e.g., see Nishioka \cite[p.~153]{N1996}) that ${\bf F}(z)$ satisfies the Mahler functional equation \begin{equation}\label{MFE}{\bf F}(z)={\bf B}(z)\, {\bf F}(z^k),\end{equation} where ${\bf B}(z):=\sum_{a=0}^{k-1}{\bf B}_a\, z^a.$ The matrix-valued function ${\bf B}(z)$ is analogous to the Fourier cocycle in the renormalisation theory of substitution and inflation systems, which carries a lot of information about features of the underlying diffraction and the spectral measures; see Bufetov and Solomyak \cite{BS2014} and Baake, G\"ahler and Ma\~nibo~\cite{BGM2019}. Note that ${\bf B}(1)={\bf B}$; this specialisation will be discussed more below. Equation \eqref{MFE} shows that the functions $F_i(z)$ behave well under the Frobenius map $z\mapsto z^k$. The functional equation \eqref{MFE} is analogous to a substitution with repeated application mirroring the iterated composition of a substitution. This property is essentially what allows us to form certain cocycles (e.g., see \eqref{AmSm}) that, under the primitivity assumption above, have convergence properties which provide for the existence of the limit measures we desire. To achieve our goal, we require a few preliminary results.

For $i\in\{0,\ldots, k-1\}$ let $\varSigma_i(n)$ be as defined in \eqref{eq:Sigmai} and set $${\bf\Sigma}(n):=(\varSigma_1(n),\varSigma_2(n),\ldots,\varSigma_d(n))^T.$$ In Section \ref{sec:notation}, we stated that the sequences $\varSigma_i(n)$ are linearly recurrent; we prove this result as the following lemma.

\begin{lemma}\label{lem:Sigma n}
If $n\geqslant 1$, then ${\bf\Sigma}(n)={\bf B}\cdot {\bf\Sigma}(n-1).$
\end{lemma}

\begin{proof} Here we are considering the sums of each $f_i$ over all integers in the interval $[k^n,k^{n+1}-1]$. These are precisely all of the integers that have $n+1$ digits in their $k$-ary expansions. Noting that the $k$-ary expansion of a non-zero integer cannot begin with a zero, we thus have that $\varSigma_i(n)$ satisfies  $$\sum_{m=k^n}^{k^{n+1}-1}f_i(m) =\sum_{m=k^n}^{k^{n+1}-1}e_i^T {\bf B}_{(m)_k}{\bf w}=e_i^T \left({\bf B}_0+{\bf B}_1+\cdots+{\bf B}_{k-1}\right)^n\sum_{a=1}^{k-1}{\bf B}_{a} {\bf w}.$$
Thus \begin{align*}{\bf\Sigma}(n)&=\left({\bf B}_0+{\bf B}_1+\cdots+{\bf B}_{k-1}\right)^n\sum_{a=1}^{k-1}{\bf B}_{a} {\bf w}\\ &=\left({\bf B}_0+{\bf B}_1+\cdots+{\bf B}_{k-1}\right){\bf\Sigma}(n-1),\end{align*} where the second equality follows from the first by writing $n=1+(n-1)$.
\end{proof}

Now, define the polynomials $b_{ij}(z)$ by ${\bf B}(z)=\left(b_{ij}(z)\right)_{1\leqslant i,j\leqslant d}$ and for each $n\geqslant 1$ define the matrix
\begin{equation}\label{eq:A matrix}
{\bf A}_n(z):=\left(\frac{\varSigma_j(n-1)\, b_{ij}(z)}{\varSigma_i(n)}\right)_{1\leqslant i,j\leqslant d}.
\end{equation} Note that since $f$ is primitive, the $\varSigma_i(n)$ cannot vanish. The matrices ${\bf A}_n(z)$ are normalised versions of the matrix ${\bf B}(z)$. In particular, they allow us to lift the result of Lemma \ref{lem:Sigma n} to the level of measures; see Proposition \ref{prop:recursion mu} below. Before proving that result, we note the following corollary of Lemma \ref{lem:Sigma n}.

\begin{corollary} The matrix ${\bf A}_n(1)$ is a Markov matrix for all $n\geqslant 1$.
\end{corollary}

For a Dirac measure $\delta_{x}$, we let $(\delta_x)^r=\delta_{rx}$ be its $r$-fold convolution product with itself. This allows one to define the value of $p(\delta_x)$, for any polynomial $p\in\mathbb{C}[z]$.

\begin{proposition}\label{prop:recursion mu}
For $n\geqslant 1$, one has
\begin{equation}\label{eq:recursion mu}
\boldsymbol{\mu}_n={\bf A}_n\big(\delta_{1/k^{n}(k-1)}\big)\, \boldsymbol{\mu}_{n-1}
\end{equation}
where ${\bf A}_n(z)$ is the matrix-valued function defined in Eq.~\eqref{eq:A matrix}.
\end{proposition}

\begin{proof}
Fix an $i\in\{1,\ldots,d\}$ and consider the $i$th entry of the vector on the right-hand side of \eqref{eq:recursion mu}. For this entry, we have
\begin{align*}
e_{i}^T{\bf A}_n\big(&\delta_{1/k^{n}(k-1)}\big)\, \boldsymbol{\mu}_{n-1}=\sum_{j=1}^d \frac{\varSigma_j(n-1)\, b_{ij}(\delta_{1/k^{n}(k-1)})}{\varSigma_i(n)} \mu_{n-1,j}\\
&=\frac{1}{\varSigma_i(n)}\sum_{j=1}^d  \sum_{\ell=0}^{k^{n}-k^{n-1}-1} b_{ij}(\delta_{1/k^{n}(k-1)}) f_{j}(k^{n-1}+\ell)\,\delta_{\ell/k^{n-1}(k-1)}\\
&=\frac{1}{\varSigma_i(n)}  \sum_{\ell=0}^{k^{n}-k^{n-1}-1}\sum_{a=0}^{k-1} \sum_{j=1}^d ({\bf B}_a)_{ij}\, (\delta_{a/k^{n}(k-1)}) f_{j}(k^{n-1}+\ell)\,\delta_{\ell/k^{n-1}(k-1)}\\
&=\frac{1}{\varSigma_i(n)}  \sum_{\ell=0}^{k^{n}-k^{n-1}-1}\sum_{a=0}^{k-1} \delta_{(k\ell+a)/k^{n}(k-1)} \sum_{j=1}^d ({\bf B}_a)_{ij} \, f_{j}(k^{n-1}+\ell)\\
&=\frac{1}{\varSigma_i(n)}  \sum_{\ell=0}^{k^{n}-k^{n-1}-1}\sum_{a=0}^{k-1} \delta_{(k\ell+a)/k^{n}(k-1)}\, f_{i}(k^{n}+k\ell+a)\\
&=\frac{1}{\varSigma_i(n)}  \sum_{m=0}^{k^{n+1}-k^{n}-1} f_{i}(k^{n}+m)\,\delta_{m/k^{n}(k-1)}\\
&=\mu_{n,i}=e_{i}^T \boldsymbol{\mu}_{n}.
\end{align*}
Here, the third equality follows using the definition of $b_{ij}(z)$, the fifth equality follows by invoking \eqref{eq:k kernel} and the sixth step is just a change of index.
\end{proof}

The following is an immediate consequence of Proposition~\ref{prop:recursion mu}.

\begin{corollary}\label{coro:recursion mu hat}
For $n\geqslant 1$, the Fourier coefficients $\widehat{\mu}_n(t)$ satisfy \begin{equation}\label{eq: recursion mu hat} \widehat{\boldsymbol{\mu}_n}(t)= {\bf A}_n\big(\mathrm{e}^{-\frac{2\pi \mathrm{i} t}{k^n(k-1)}}\big)\, \widehat{\boldsymbol{\mu}_{n-1}}(t).
\end{equation}
for all $t\in\mathbb{Z}$.
\end{corollary}

\begin{remark}
While the convergents $\mu_n$ are pure point probability measures on $\mathbb{T}$, one can also consider the measures $\nu_n=\delta_\mathbb{Z}\ast \mu_n$, which are $\mathbb{Z}$-periodic measures in $\mathbb{R}$. Here, one has $\widehat{\nu}_n\colon \mathbb{R}\to \mathbb{C}$, with $\widehat{\mathbb{R}}=\mathbb{R}$. In the discussion below, we only carry out the analysis for $\mu_n$, but the (vague) convergence of the relevant matrix products also hold for $\nu_n$. \hfill $\Diamond$
\end{remark}

Using Equation~\eqref{eq:recursion mu}, one can construct the infinite matrix convolution and hope for the existence of the limit vector
\begin{equation}\label{eq: limit measure}
\boldsymbol{\mu}:=\bigg(\DOTSB\bigop{\ast}_{n=1}^{\infty} {\bf A}_{n} \big(\delta_{1/k^n(k-1)}\big)\bigg)\,\boldsymbol{\mu}_0.
\end{equation}
The existence of the limit vector $\boldsymbol{\mu}$ depends on the compact convergence of the Fourier coefficients, which boils down to the compact convergence of the analytic matrix product
\begin{equation*}
\prod_{n=1}^{\infty}  {\bf A}_n\big(\mathrm{e}^{-\frac{2\pi \mathrm{i} t}{k^n(k-1)}}\big),
\end{equation*} a property we demonstrate below.

We require the following lemma that the positivity of ${\bf B}$ implies the convergence of the quotients $\varSigma_i(n-1)/\varSigma_j(n)$ to a positive value.

\begin{lemma}\label{lem: primitive quotients} Let $f$ be a primitive $k$-regular sequence.  Then for each $i,j\in\{1,\ldots,d\}$, the limit of the quotient $\varSigma_j(n-1)/\varSigma_i(n)$ exists and is non-zero.
\end{lemma}

\begin{proof} This follows from the fact that ${\bf B}$ is positive and that $f$ is non-negative (and not trivial), so that given an $i\in\{1,\ldots,d\}$ there is a constant $c_i>0$ such that \begin{equation}\label{cisigma}\varSigma_i(n)\sim c_i\cdot\rho_{_{\rm PF}}({\bf B})^n,\end{equation} as $n\to \infty$, where the positive real number $\rho_{_{\rm PF}}({\bf B})$ is the Perron--Frobenius eigenvalue of ${\bf B}$. To see this, let $\rho_{_{\rm PF}}({\bf B})$ be the Perron--Frobenius eigenvalue of ${\bf B}$ and let ${\bf v}$ be the corresponding positive Perron--Frobenius eigenvector. Let $j$ be a positive integer so that ${\bf B}^j {\bf w}$ is entry-wise greater than ${\bf v}$. Then for $n>j$, \begin{multline*}\varSigma_i(n)=e_i^T {\bf B}^n\, {\bf w}=e_i^T {\bf B}^{n-j}\, {\bf B}^{j}\, {\bf w}=e_i^T {\bf B}^{n-j}\, {\bf v}+e_i^T {\bf B}^{n-j}\,  ({\bf B}^{j}\,{\bf w}-{\bf v})\\ \geqslant e_i^T {\bf B}^{n-j}\, {\bf v}=e_i^T \rho_{_{\rm PF}}({\bf B})^{n-j}\, {\bf v}=(e_i^T{\bf v}/\rho_{_{\rm PF}}({\bf B})^j)\cdot \rho_{_{\rm PF}}({\bf B})^n.\end{multline*} Since $\rho_{_{\rm PF}}({\bf B})$ is a simple eigenvalue of ${\bf B}$ with maximal modulus, comparing the eigenvalue expansion of the linear recurrent sequence $\varSigma_i(n)$ with this inequality proves \eqref{cisigma}. The lemma follows immediately.
\end{proof}

For the proof of the following proposition, we adapt a technique used by Baake and Grimm in determining the intensities of Bragg peaks of Pisot substitutions via the internal Fourier cocycle \cite[Thm.~4.6]{BG2019}.

\begin{proposition}\label{prop:conv general}
Let $f$ be a primitive $k$-regular sequence. Then the matrix product $$\prod_{n=1}^{\infty} {\bf A}_n\big(\mathrm{e}^{-\frac{2\pi \mathrm{i} t}{k^n(k-1)}}\big),$$ as a function of $t$, converges uniformly on $\mathbb{T}$ and compactly on $\mathbb{R}$.
\end{proposition}

\begin{proof}
It follows from the convergence of the quotients that ${\bf A}_n(z)\to {\bf A}(z)$ as $n\to\infty$, with $({\bf A}(z))_{ij}=c_{ij}b_{ij}(z)$. In particular, one has ${\bf A}_n(1)\to {\bf A}(1)$, where every ${\bf A}_n(1)$ is a Markov matrix.
Setting ${\bf S}_{n}(t):={\bf A}_n\big(\mathrm{e}^{\frac{-2\pi it}{(k-1)}}\big)$, one gets the convergence
${\bf S}_n(t)\to {\bf S}(t):={\bf A}(\mathrm{e}^{\frac{-2\pi it}{(k-1)}})$.
One has
\begin{equation}\label{AmSm}{\bf A}^{(m)}\big(\mathrm{e}^{\frac{-2\pi it}{(k-1)}}\big)={\bf S}^{(m)}(t):={\bf S}_m\bigg(\frac{t}{k^m}\bigg)\, {\bf S}_{m-1}\bigg(\frac{t}{k^{m-1}}\bigg)\cdots\, {\bf S}_1\bigg(\frac{t}{k}\bigg).\end{equation}
Note that ${\bf S}^{(m)}(0)={\bf S}_m(0)\, {\bf S}_{m-1}(0)\cdots{\bf S}_1(0)$ is also a Markov matrix.

For $0\leqslant r\leqslant \ell$ define the product
\[
{\bf S}^{(\ell,r)}(t):={\bf S}_{\ell}\bigg(
\frac{t}{k^\ell}\bigg)\cdots\, {\bf S}_{r}\bigg(
\frac{t}{k^r}\bigg).
\]
These products satisfy the identity
${\bf S}^{(\ell,r)}(t)={\bf S}^{(\ell,c)}(t)\, {\bf S}^{(c,r)}(t)$,
for $0\leqslant r\leqslant c\leqslant \ell$.
Note that ${\bf S}^{(\ell,0)}(t)={\bf S}^{(\ell)}(t)$. Moreover,
${\bf S}^{(\ell,r)}(0)$ is a Markov matrix for any such $\ell,r$, being a finite product of Markov matrices.
The first goal is to show that ${\bf S}^{(m,\ell)}(t)$ is equicontinuous at $t=0$ for all $0\leqslant\ell\leqslant m$.

Employing a telescoping argument, one obtains the equality
\[
{\bf S}^{(m,\ell)}(t)-{\bf S}^{(m,\ell)}(0)=
\sum_{j=\ell-1}^{m-1} \left({\bf S}^{(m,j+2)}(0)\, {\bf S}^{(j+1,\ell)}(t)-
{\bf S}^{(m,j+1)}(0)\, {\bf S}^{(j,\ell)}(t)\right)
\]
This implies the following
\begin{align*}
\|{\bf S}^{(m,\ell)}(t)&-{\bf S}^{(m,\ell)}(0)\|_{\infty}\\ &\leqslant \sum_{j=0}^{m-1}
\|{\bf S}^{(m,j+2)}(0)\,{\bf S}^{(j+1,\ell)}(t)-
{\bf S}^{(m,j+1)}(0)\, {\bf S}^{(j,\ell)}(t)\|_{\infty}\\
&\leqslant \sum_{j=0}^{m-1} \|{\bf S}^{(m,j+2)}(0)\|_{\infty}\cdot
\|{\bf S}^{(j+1,\ell)}(t)-{\bf S}_{j+1}(0)\, {\bf S}^{(j,\ell)}(t) \|_{\infty}\\
&\leqslant \sum_{j=0}^{m-1} \|{\bf S}^{(m,j+2)}(0)\|_{\infty}\cdot
\|{\bf S}_{j+1}(\tfrac{t}{k^{j+1}})-{\bf S}_{j+1}(0)\|_{\infty}\cdot \|{\bf S}^{(j,\ell)}(t)\|_{\infty}\\
&\leqslant \sum_{j=0}^{m-1}
\|{\bf S}_{j+1}(\tfrac{t}{k^{j+1}})-{\bf S}_{j+1}(0)\|_{\infty},
\end{align*}
where the last step follows from the properties $\|{\bf S}^{(m,j+2)}(0)\|_{\infty}=1$ and $$\|{\bf S}^{(j,\ell)}(t)\|_{\infty}\leqslant \|{\bf S}^{(j,\ell)}(0)\|_{\infty}=1,$$ since the matrices ${\bf B}_a$ are non-negative and since both ${\bf S}^{(j,\ell)}(t)$ and ${\bf S}^{(j,\ell)}(0)$ are Markov.

Now, let $\varepsilon>0$ be given and choose $\delta>0$ such that
\[\|{\bf S}_{j+1}(\tfrac{t}{k^{j+1}})-{\bf S}_{j+1}(0)\|_{\infty}<\frac{\varepsilon}{k^{j+1}}\]
holds whenever $|\tfrac{t}{k^{j+1}}|<\delta$. This yields
\begin{equation}\label{eq:equicont general}
\|{\bf S}^{(m,\ell)}(t)-{\bf S}^{(m,\ell)}(0)\|_{\infty}
 \leqslant \sum_{j=0}^{m-1}
\|{\bf S}_{j+1}(\tfrac{t}{k^{j+1}})-{\bf S}_{j+1}(0)\|_{\infty}
< \sum_{j=0}^{\infty} \frac{\varepsilon}{k^{j}}
=\varepsilon\big(\frac{k}{k-1}\big)
\end{equation}
which proves that ${\bf S}^{(m,\ell)}(t)$ is equicontinuous at $t=0$.

To show compact convergence, we prove that ${\bf S}^{(m)}(t)={\bf S}^{(m,0)}(t)$ is uniformly Cauchy for any given compact set $K\subseteq \mathbb{Z}_{>0}$. To this end, choose $p$ such that $|\tfrac{t}{k^p}|<\delta$ for all $t\in K$, which implies $|\tfrac{t}{k^j}|<\delta$ for all $j\geqslant p$ and $t\in K$.  One then has
\begin{align*}
\|{\bf S}^{(p+q+r)}(t)-{\bf S}^{(p+q)}(t)\|_{\infty} & \leqslant
\|{\bf S}^{(p+q+r,p+1)}(t)-{\bf S}^{(p+q,p+1)}(t)\|_{\infty}\cdot\|{\bf S}^{(p)}(t)\|_{\infty}\\
& \leqslant \|{\bf S}^{(p+q+r,p+1)}(t)-{\bf S}^{(p+q,p+1)}(t)\|_{\infty},
\end{align*}
since $\|{\bf S}^{(p)}(t)\|\leqslant 1$. Thus, using the triangle inequality, we obtain
\begin{multline*}
\|{\bf S}^{(p+q+r,p+1)}(t)-{\bf S}^{(p+q,p+1)}(t)\|_{\infty}
\leqslant\|{\bf S}^{(p+q+r,p+1)}(t)-{\bf S}^{(p+q+r,p+1)}(0)\|_{\infty}\\
+
\|{\bf S}^{(p+q+r,p+1)}(0)-{\bf S}^{(p+q,p+1)}(0)\|_{\infty}+
\|{\bf S}^{(p+q,p+1)}(t)-{\bf S}^{(p+q,p+1)}(0)\|_{\infty}
\end{multline*}
where the first and the third summands are strictly less than $\varepsilon\cdot k/(k-1)$
by invoking Eq.~\eqref{eq:equicont general} for large enough $p$.
The second summand splits further, into
\begin{multline}
\|{\bf S}^{(p+q+r,p+1)}(0)-{\bf S}^{(p+q,p+1)}(0)\|_{\infty}\leqslant
\|{\bf S}^{(p+q+r,p+1)}(0)-{\bf A}(1)^{q+r}\|_{\infty}\\+\|{\bf A}(1)^{q+r}-{\bf A}(1)^{q}\|_{\infty}+\|{\bf S}^{(p+q,p+1)}(0)-{\bf A}(1)^q\|_{\infty}, \label{eq: sum S(0)}
\end{multline}
Since ${\bf S}_{\ell}(0)$ converges to ${\bf A}(1)$, one can choose $p^{\prime}$ large enough such that
\[\|{\bf S}_{\ell}(0)-{\bf A}(1)\|_{\infty}<\frac{\varepsilon}{k^\ell}\] for all $\ell\geqslant p^{\prime}$. One can then invoke a similar argument to show that the first and the third summands in Eq.~\eqref{eq: sum S(0)} are strictly less than $\varepsilon\cdot k/(k-1)$ whenever $p\geqslant p^{\prime}$. Finally, since ${\bf A}(1)$ is Markov, ${\bf A}(1)$ converges to a steady state matrix whence one can choose $q$ large enough such that
the second summand is less than $\varepsilon$.

Since $r$ is arbitrary, this means for all $m,n\geqslant \max\left\{p,p^{\prime}\right\}+q$ one has that
\[
\|{\bf S}^{(m)}(t)-{\bf S}^{(n)}(t)\|_{\infty}<\varepsilon\left(\frac{5k-1}{k-1}\right),
\]
and thus, ${\bf S}^{(n)}(t)$ converges uniformly on $\mathbb{T}$ and compactly on $\mathbb{R}$.
\end{proof}

The following result is now an immediate corollary  of Proposition \ref{prop:conv general} by L\'evy's continuity theorem \cite[Thm.~3.14]{BFbook}.

\begin{theorem}\label{thm:muexists}
Let $f$ be a primitive $k$-regular sequence. Then the weak limit measure vector $\boldsymbol \mu$ in Eq.~\eqref{eq: limit measure} exists.
\end{theorem}

\begin{remark}
As stated previously, the non-negativity assumptions on the $k$ matrices ${\bf B}_a$ and the positivity (primitivity, initially) assumption on ${\bf B}$ are natural, especially if one views it as the corresponding analogue of the substitution matrix ${\bf M}_\varrho$ for shift spaces arising from a substitution $\varrho$. For the latter, the primitivity of ${\bf M}_\varrho$ implies that the hull $\mathbb{X}_{\varrho}$ defined by $\varrho$ is strictly ergodic, and hence the diffraction measure $\widehat{\gamma}_{w}$ is constant on $\mathbb{X}_{\varrho}$ given an arbitrary choice of weight function $w$; see Baake's and Grimm's monograph \cite{BGbook} for details and definitions.
Alternatively, one can consider matrix Riesz products of measures, where the approximants are absolutely continuous with respect to Lebesgue measure, and weak convergence to the matrix of correlation measures is likewise guaranteed by primitivity; see Bufetov and Solomyak~\cite[Lem.~2.2]{BS2014} and Queff\'elec \cite[Thm.~8.1]{Qbook}. \hfill $\Diamond$
\end{remark}

In what follows, we recover an analogous uniqueness result for the basis $\left\{f_1,\ldots,f_d\right\}$ of the $\mathbb{Q}$-vector space $\mathcal{V}_k(f)$.

\begin{proposition}\label{prop: uniqueness}
Let $f$ be a primitive $k$-regular sequence. Then $\boldsymbol{\mu}=(\mu_f,\ldots,\mu_f)^{T}$ where $\mu_f$ is a probability measure on $\mathbb{T}$. That is, for each $i\in\{1,\ldots,d\}$, the weak limit of $\mu_{n,i}$ is $\mu_f$.
\end{proposition}

\begin{proof} For a fixed $t\in\mathbb{R}$, one has $\widehat{\boldsymbol{\mu}}(t)={\bf L}(t)\, \widehat{\boldsymbol{\mu}}_0(t)$, where $${\bf L}(t):=\prod_{n=1}^{\infty}  {\bf A}_n\big(\mathrm{e}^{-\frac{2\pi \mathrm{i} t}{k^n(k-1)}}\big).$$ From the convergence in Lemma~\ref{lem: primitive quotients}, ${\bf A}_{n}(1)\to {\bf A}(1)$ as $n\to \infty$, where ${\bf A}(1)$ is a primitive Markov matrix, where the primitivity of ${\bf A}(1)$ follows from the positivity of ${\bf B}$ and Lemma~\ref{lem: primitive quotients}.
This means ${\bf A}(1)^n$ converges to the rank-$1$ projector $ {\bf P}_{\bf A}$ corresponding to the eigenvector ${\bf 1}=(1,\ldots,1)^{T}$. One then gets the equality
$${\bf A}(1)\, {\bf L}(t)={\bf L}(t)= {\bf P}_{\bf A}\, {\bf L}(t),$$ which implies $\widehat{\boldsymbol{\mu}}(t)={\bf P}_{\bf A}\, \widehat{\boldsymbol{\mu}}(t)=c(t)\, {\bf 1}$.
This means that the limit measures $\mu_i$ have the same Fourier coefficients for all $t$, and hence must correspond to the same measure $\mu=\mu_f$ for all $1\leqslant i\leqslant d$.
\end{proof}

Theorem~\ref{main1} follows by combining Theorem~\ref{thm:muexists} and Proposition \ref{prop: uniqueness}.

\section{A natural probability measure associated to $f$}\label{sec:f}

In this section we will prove Theorem \ref{main2}, via the general situation described at the end of Section \ref{sec:notation}. To set up, we recall the following notation, now for only a single $k$-regular sequence $f$, $$\varSigma_f(n):=\sum_{m=k^n}^{k^{n+1}-1}f(m)$$ and $$
   \mu^{}_{f,n} \, := \, \frac{1}{\varSigma_f(n)} \sum_{m=0}^{k^{n+1}-k^n - 1}
   f(k^n + m) \, \delta^{}_{m / k^n(k-1)},
$$
where $\delta_x$ denotes the unit Dirac measure at $x$. Also as previously, let $f$ be defined by the matrices ${\bf B}_0,\ldots,{\bf B}_{k-1}$ and the vector ${\bf w}\in\mathbb{Z}^{d\times 1}$ such that $f(m)=e_1^T {\bf B}_{(m)_k} {\bf w}$, where $(m)_k=i_s\cdots i_1i_0$ is the base-$k$ expansion of $m$ and ${\bf B}_{(m)_k}:={\bf B}_{i_0}{\bf B}_{i_1}\cdots {\bf B}_{i_s}.$ As before, set ${\bf B}:=\sum_{a=0}^{k-1}{\bf B}_a.$ Let $\rho({\bf M})$ denote the spectral radius of the matrix ${\bf M}$ and denote the {\em joint spectral radius} of a finite set of matrices $\{{\bf M}_1,{\bf M}_2,\ldots, {\bf M}_{\ell}\}$, by the real number $$\rho^*(\{{\bf M}_1,{\bf M}_2,\ldots, {\bf M}_{\ell}\})=\limsup_{n\to\infty}\max_{1\leqslant i_0,i_1,\ldots,i_{n-1}\leqslant \ell}\left\| {\bf M}_{i_1}{\bf M}_{i_2}\cdots{\bf M}_{i_{n}}\right\|^{1/n},$$ where $\|\cdot\|$ is any (submultiplicative) matrix norm. This quantity was introduced by Rota and Strang \cite{RS1960} and has a wide range of applications. For an extensive treatment, see Jungers's monograph \cite{J2009}.

Unlike Theorem \ref{main1}, so also unlike the previous section, we do not yet assume that $f$ is non-negative, nor that ${\bf B}$ is positive. However, to avoid degeneracies (discussed in the next section), we assume that the spectral radius $\rho({\bf B})$ is the unique dominant eigenvalue of ${\bf B}$, that $$\rho:=\rho({\bf B})>\rho^*(\{{\bf B}_0,\ldots,{\bf B}_{k-1}\})=:\rho^*,$$ that for $n$ large enough $\varSigma_f(n)\neq 0$ and that the asymptotical behaviour of $\varSigma_f(n)$ is determined by $\rho({\bf B})$. This last assumption will be highlighted and made explicit below.

To exploit the asymptotical behaviour of $\varSigma_f(n)$, we use a result of Dumas \cite[Thm.~3]{D2013} on the asymptotic nature of the partial sums $\sum_{m\leqslant x}f(m)$. To set up Dumas' result, we require some further notation. Let ${\bf J}_\rho=\rho{\bf I}_v+{\bf Z}$ be the Jordan block associated to $\rho$ from the Jordan normal form of ${\bf B}$, where ${\bf I}_v$ is the $v\times v$ identity matrix and ${\bf Z}$ is the nilpotent matrix of index $v$ with ones on the superdiagonal and zeros elsewhere.  Let ${\bf V}_\rho$ be the $d\times v$ matrix whose columns are the elements associated to $\rho$ from some Jordan basis of ${\bf B}$. Let the vector ${\bf w}_\rho$ be such that ${\bf V}_\rho {\bf w}_\rho$ equals the component of ${\bf w}$ in the invariant subspace of ${\bf B}$ associated to $\rho$. Finally, we define the matrix-valued function ${\bf F}_\rho:\mathbb{R}\to\mathbb{C}^{d\times v}$ by \begin{equation}\label{dilation}{\bf F}_{\rho}(x) \cdot {\bf J}_{\rho} = \sum_{a=0}^{k-1}{\bf B}_{a} \cdot {\bf F}_{\rho}(k x-a),\end{equation} with the boundary conditions $${\bf F}_{\rho}(x)=\begin{cases}0 & \text{ for } x \leqslant 0\\ {\bf V}_{\rho} &\text{ for } x \geqslant 1.\end{cases}$$ The function ${\bf F}_\rho$ exists and is unique since $\rho>\rho^*$. Moreover, the function ${\bf F}_\rho$ is H\"older continuous with exponent $\alpha$ for any $\alpha<\log_k(\rho/\rho^*).$ Functional equations such as \eqref{dilation} are known as \emph{dilation equations} or \emph{two-scale difference equations} in the literature; seminal work on these was done by Daubechies and Lagarias \cite{DL1991,DL1992}. See also Micchelli and Prautzsch \cite{MP1989}. A point of interest for our context is that both the above-mentioned papers of Daubechies and Lagarias as well as the seminal paper of Allouche and Shallit \cite{AS1992} introducing $k$-regular sequences were published within the span of one year. Twenty years later, Dumas \cite{D2013} connected these ideas by showing explicitly how one can use a dilation equation to determine the  asymptotic growth of the partial sums of a regular sequence. We record his result here in the special case fit for our purpose.

\begin{theorem}[Dumas]\label{thm:Dumas} Suppose that the spectral radius $\rho=\rho({\bf B})$ is the unique dominant eigenvalue of ${\bf B}$ and that $\rho({\bf B})>\rho^*(\{{\bf B}_0,\ldots,{\bf B}_{k-1}\})$. Then $$\sum_{m\leqslant x}f(m) = e_1^T\, {\bf E}_\rho(\log_k(x))\, {\bf M}_\rho + O(x^{\log_k (r)}),$$ where $r$ is any positive number strictly between $\rho$ and the modulus of the next-largest eigenvalue of ${\bf B}$ and \begin{multline*}{\bf E}_\rho(x):=({\bf I}_d-{\bf B}_0)\, {\bf V}_\rho\, ({\bf I}_v-{\bf J}_\rho)^{-1}\\+\big(-({\bf I}_d-{\bf B}_0)\, {\bf V}_\rho\, ({\bf I}_v-{\bf J}_\rho)^{-1}+{\bf F}_\rho(k^{\{ x\}-1})\big)\cdot\rho^{\lfloor x\rfloor+1}\cdot {\bf P}_\rho(\lfloor x\rfloor),\end{multline*} where ${\bf P}_\rho(\lfloor x\rfloor)=({\bf I}_v+(1/\rho){\bf Z})^{\lfloor x\rfloor}$. \end{theorem}

A point to be made here is that, while Theorem \ref{thm:Dumas} is certainly technical, in the cases we apply it, for integers $x\in[k^n, k^{n+1})$, the integer part $\lfloor x\rfloor=n$ is constant, so only the dependence on the fractional parts $\{x\}$ will need to be dealt with. For a detailed example of how this theorem  can be applied to give a distribution function, see Baake and Coons~\cite[Sec.~3]{BC2018}, where they give an account concerning the Stern sequence.

Now, applying Theorem \ref{thm:Dumas} to the complete sums $\varSigma_f(n)$ gives $$\varSigma_f(n)=\sum_{m=k^n}^{k^{n+1}-1}f(m)
=\rho^{n+1} e_1^T \left({\bf F}_\rho\bigg(\frac{k-1/k^n}{k}\bigg)-{\bf F}_\rho\bigg(\frac{1}{k}\bigg)\right){\bf P}_\rho(n)\, {\bf M}_\rho+O(r^n),$$ for some $r<\rho$. As $n\to\infty$, since ${\bf F}_\rho$ is H\"older continuous, $${\bf F}_\rho\bigg(\frac{k-1/k^n}{k}\bigg)={\bf V}_\rho+o(1),$$ so that \begin{equation}\label{eq:denom}\varSigma_f(n)
=\rho^{n+1} e_1^T \left({\bf V}_\rho-{\bf F}_\rho\bigg(\frac{1}{k}\bigg)\right){\bf P}_\rho(n)\, {\bf M}_\rho+o(\rho^n).\end{equation}
For $o(\rho^n)$ to be the true error term in Eq.~\eqref{eq:denom}, for sufficiently large $n$ we must have that \begin{equation}\label{condition} e_1^T \left({\bf V}_\rho-{\bf F}_\rho\bigg(\frac{1}{k}\bigg)\right){\bf P}_\rho(n)\, {\bf M}_\rho \neq 0.\end{equation} Since $${\bf P}_\rho(n)=({\bf I}_v+(1/\rho){\bf Z})^{n}=\sum_{j=0}^{v-1}{n\choose j}\rho^{j}{\bf Z}^j,$$ the non-zero condition \eqref{condition} holds precisely when there is a $j\in\{0,\ldots,v-1\}$ such that \begin{equation}\label{zjcond} e_1^T \left({\bf V}_\rho-{\bf F}_\rho\bigg(\frac{1}{k}\bigg)\right){\bf Z}^j\, {\bf M}_\rho \neq 0.\end{equation} This being the case, if $\ell$ is the largest $j\in\{0,\ldots,v\}$ for which \eqref{zjcond} holds, then \begin{equation}\label{sigfnl} \varSigma_f(n)
=\rho^{n-\ell+1}{n\choose \ell} e_1^T \left({\bf V}_\rho-{\bf F}_\rho\bigg(\frac{1}{k}\bigg)\right){\bf Z}^\ell\, {\bf M}_\rho+o(\rho^nn^\ell).\end{equation}

Arguing as in the previous paragraph allows us to prove the following result.

\begin{theorem}\label{main2explicit} Let $f$ be an integer-valued $k$-regular sequence. Suppose that $\rho({\bf B})$ is the unique dominant eigenvalue of ${\bf B}$, $\rho({\bf B})>\rho^*(\{{\bf B}_0,\ldots,{\bf B}_{k-1}\})$ and that Eq.~\eqref{zjcond} holds. Then, the limit $\mu_{f}([0,x])$ of the sequence $\mu_{f,n}([0,x])$ exists. Moreover, the function $\mu_{f}([0,x])$ is H\"older continuous with exponent $\alpha$ for any $\alpha<\log_k(\rho/\rho^*)$.
\end{theorem}

\begin{proof} Let $x\in\mathbb{T}$ and consider the sequence of functions $\mu_{f,n}([0,x])$. Then, applying the argument of the above paragraph with $(k-1/k^n)/k$ replaced by $1+(k-(1+1/k^n))x$, we have \begin{multline*}\varSigma_f(n)\cdot \mu_{f,n}([0,x])=\sum_{m=k^n}^{k^n(1+(k-(1+1/k^n))x)}f(m)\\
=\rho^{n+1} e_1^T \left({\bf F}_\rho\bigg(\frac{1+(k-(1+1/k^n))x}{k}\bigg)\right.\\
\left.-{\bf F}_\rho\bigg(\frac{1}{k}\bigg)\right){\bf P}_\rho(n)\, {\bf M}_\rho+o(\rho^n),
\end{multline*} for some $r<\rho$. As before, using the H\"older continuity of ${\bf F}_\rho$ and letting $\ell$ be the largest value for which \eqref{zjcond} holds, we obtain that $\varSigma_f(n)\cdot \mu_{f,n}([0,x])$ equals $$
\rho^{n-\ell+1} {n\choose \ell}e_1^T \left({\bf F}_\rho\bigg(\frac{1+(k-1)x}{k}\bigg)-{\bf F}_\rho\bigg(\frac{1}{k}\bigg)+o(1)\right){\bf Z}^\ell\, {\bf M}_\rho+o(\rho^n n^\ell).
$$ Using the asymptotic \eqref{sigfnl} yields \begin{equation*} \mu_{f,n}([0,x])=\frac{e_1^T \left({\bf F}_\rho\big(\frac{1+(k-1)x}{k}\big)-{\bf F}_\rho\big(\frac{1}{k}\big)+o(1)\right){\bf Z}^\ell\, {\bf M}_\rho+o(1)}{e_1^T \left({\bf V}_\rho-{\bf F}_\rho\big(\frac{1}{k}\big)\right){\bf Z}^\ell\, {\bf M}_\rho+o(1)}.\end{equation*} By our assumption that \eqref{zjcond} holds, the denominator limits to a nonzero value, and so the point-wise limit of $\mu_{f,n}([0,x])$ exists for all $x$; explicitly, \begin{equation*} \mu_{f}([0,x])=\lim_{n\to\infty}\mu_{f,n}([0,x])=\frac{e_1^T \left({\bf F}_\rho\big(\frac{1+(k-1)x}{k}\big)-{\bf F}_\rho\big(\frac{1}{k}\big)\right){\bf Z}^\ell\, {\bf M}_\rho}{e_1^T \left({\bf V}_\rho-{\bf F}_\rho\big(\frac{1}{k}\big)\right){\bf Z}^\ell\, {\bf M}_\rho}.\end{equation*} Finally, we note that the function $\mu_{f}([0,x])$ is H\"older continuous with exponent $\alpha$ for any $\alpha<\log_k(\rho/\rho^*)$, a property it inherits directly from ${\bf F}_\rho$.
\end{proof}

While we have used the notation $\mu_{f}([0,x])$, the assumptions of Theorem \ref{main2explicit} are not strong enough to guarantee the existence of a measure $\mu_f$ for which $\mu_{f}([0,x])$ is a distribution function. The following result contains a sufficient condition and is a generalisation of the result stated in Section~\ref{sec:notation} as Theorem \ref{main2}.

\begin{theorem}\label{BVexplicit} Let $f$ be an integer-valued $k$-regular sequence. Suppose that $\rho({\bf B})$ is the unique dominant eigenvalue of ${\bf B}$, $\rho({\bf B})>\rho^*(\{{\bf B}_0,\ldots,{\bf B}_{k-1}\})$ and that Eq.~\eqref{zjcond} holds. Suppose that $\mu_{f}([0,x])$, as provided by Theorem \ref{main2explicit}, is a function of bounded variation. Then $\mu_{f}([0,x])$ is the distribution function of a measure $\mu_f$, which is continuous with respect to Lebesgue measure.
\end{theorem}

\begin{proof} We start with the function $\mu_{f}([0,x])$ provided by Theorem \ref{main2explicit}. Assuming that $\mu_{f}([0,x])$ is of bounded variation, we form the (possibly) signed measure $\mu_f$ via a Riemann--Stieltjes integral, assigning the value $$\int_a^b {\rm d} \mu_f([0,x]):=\mu_{f}([0,b])-\mu_{f}([0,a])$$ to any interval $(a,b]\subseteq\mathbb{T}$; see Lang \cite[Chp.~X]{Lang} for details on Riemann--Stieltjes integration and measure. This is a Borel measure and has distribution function $\mu_{f}([0,x])$. Since the distribution functions $\mu_{f,n}([0,x])$ of the measures $\mu_{f,n}$ are converging point-wise to the continuous distribution function $\mu_{f}([0,x])$ of the measure $\mu_{f}$, the measures $\mu_{f,n}$ are converging weakly to $\mu_{f}$. Moreover, the measure $\mu_{f}$ has no pure points since $\mu_{f}([0,x])$ is continuous, thus $\mu_{f}$ is continuous with respect to Lebesgue measure.
\end{proof}

This result gives Theorem \ref{main2} as a corollary.

\begin{proof}[Proof of Theorem \ref{main2}] This follows immediately from Theorem \ref{main2explicit} using the fact that the partial sums $\sum_{m\leqslant y}f(m)$ are  increasing with $y$ to show that $\mu_f([0,x])$ is increasing, so of bounded variation.
\end{proof}

\section{Some comments on assumptions}\label{nonprim}

Theorem \ref{main1} relies on one main assumption, the primitivity of $f$, and Theorem~\ref{BVexplicit} relies on a two key assumptions. In this section, we consider these assumptions, highlighting their necessity with examples.

The assumption of primitivity of the $k$-regular sequence $f$ in Theorem \ref{main1} implies that the matrix ${\bf B}$ is positive; recall, this positivity was gained by starting with a primitive matrix and considering $f$ as a $k^j$-regular sequence, where $j$ was the minimal positive integer such that ${\bf B}^j$ is positive. Sometimes such a choice is not immediately possible. For example, consider the Josephus sequence $J$, which is $2$-regular and determined by $J(0)=0$ and the recursions $J(2n)=2J(n)-1$ and $J(2n+1)=2J(n)+1$. These recursions imply that the sequences $J$ and ${\bf 1}$ (the constant sequence) form a basis for $\mathcal{V}_2(J)$. With this basis, we arrive at the linear representation $J(n)=e_1^T\, {\bf B}_{(n)_2} (0,1)^T,$ where $${\bf B}_0 =\begin{pmatrix}
2 &-1 \\
0 & 1
\end{pmatrix}\quad\mbox{and}\quad {\bf B}_1 =\begin{pmatrix}
2 &1 \\
0 & 1
\end{pmatrix},\quad\mbox{so that}\quad{\bf B}=\begin{pmatrix}
4 &0 \\
0 & 2
\end{pmatrix}.$$ In this case, ${\bf B}$ is neither positive nor primitive, so Theorem \ref{main1} cannot be applied. However, according to Theorem \ref{main2}, one can apply the same limiting measures $\mu_{n,i}$ via the results in Section \ref{sec:f}. In the limit one gets $\mu_{J}=h(x)\cdot{\lambda}$ whereas $\mu_{\bf 1}={\lambda}$, where $\lambda$ is normalised Haar measure on $\mathbb{T}$ and $h(x)=2x$ is the Radon--Nikodym density for the limit measure associated to
$J$.

But, there is hope here---a change of basis allows the use of Theorem \ref{main1}. We need only note that if we conjugate ${\bf B}_0$ and ${\bf B}_1$ by the matrix $$\begin{pmatrix}
1 &-1 \\
1 & 1
\end{pmatrix},$$ we arrive at a new linear representation of the sequence $J$, and a new basis for $\mathcal{V}_2(J)$, where the matrix ${\bf B}$ is replaced by $$\begin{pmatrix}
1 &-1 \\
1 & 1
\end{pmatrix} {\bf B}\begin{pmatrix}
1 &-1 \\
1 & 1
\end{pmatrix}^{-1}=\begin{pmatrix}
1 &-1 \\
1 & 1
\end{pmatrix} \begin{pmatrix}
4 &0 \\
0 & 2
\end{pmatrix}\begin{pmatrix}
1 &-1 \\
1 & 1
\end{pmatrix}^{-1}=\begin{pmatrix}
3 &1 \\
1 & 3
\end{pmatrix},$$ which is positive. Thus, we can apply Theorem \ref{main1} using this new basis resulting in the vector of measures $\boldsymbol{\mu}_J=\mu_J(1,1)^T.$ This method has forgotten the measure $\mu_{\bf 1}$. This process, in general (and if possible), picks out the maximal measure, in the sense of growth of elements of $\mathcal{V}_k(f)$. This maximality becomes evident noticing that one can choose a basis for $\mathcal{V}_k(f)$ from ${\rm ker}_k(f)$; see Allouche and Shallit \cite[Thm.~2.2(b)]{AS1992}. For more information on the Josephus sequence and for a study on measures associated to affine $2$-regular sequences, see Evans \cite{Epre}.

Shifting now to Theorem \ref{BVexplicit}, consider the assumption in that $\rho({\bf B})$ is the unique dominant eigenvalue of ${\bf B}$ (which holds if ${\bf B}$ were say, primitive). There are plenty of examples of regular sequences where this is not the case. Dumas \cite[Ex.~5]{D2013} gave an interesting example of a $2$-regular sequence, which we denote by $D$, whose associated matrix ${\bf B}$ has a negative integer entry, is not primitive and does not have a maximal eigenvalue. Dumas' sequence $D$ is defined  by setting $${\bf B}_0=\left(\begin{matrix} 1 & 0\\ 0& 1\end{matrix}\right)\quad\mbox{and}\quad {\bf B}_1=\left(\begin{matrix} 3 & -3\\ 3 &3\end{matrix}\right),$$ and \begin{equation}\label{D}D(m)=e_1^T\, {\bf B}_{(m)_2} e_1=e_1^T\, {\bf B}_{1}^{s_2(m)} e_1,\end{equation} where $s_2(m)$ is the sum of the bits of $m$. As in the previous section, we consider the sequence of pure point measures  \begin{equation}\label{muDn}\mu_{D,n}:=\frac{1}{\varSigma_D(n)}\sum_{m=2^{n}}^{2^{n+1}-1}D(m)\, \delta_{m/2^{n+1}},\end{equation} where \begin{equation}\label{SigmaDn}\varSigma_D(n):=\sum_{m=2^{n}}^{2^{n+1}-1}D(m)=\frac{3}{2}\cdot\big((1+\mathrm{i})(4+3\mathrm{i})^n+(1-\mathrm{i})(4-3\mathrm{i})^n\big),\end{equation} since in this case, the matrix ${\bf B}={\bf B}_0+{\bf B}_1$ has eigenvalues $4+3\mathrm{i}$ and $4-3\mathrm{i}$, each of modulus $5$. To show that the `limit' $\mu_D$ does not exist, it is enough to prove that the sequence $\big(\mu_{D,n}([0,1/2))\big)_{n\geqslant 0}$ does not have a limit. To see why this is enough, note that if the measure $\mu_{D}$ did exist, then so would its distribution function $\mu_D([0,x)):[0,1]\to\mathbb{R}$, and necessarily we would have that $\lim_{n\to\infty}\mu_{D,n}([0,1/2))=\mu_D([0,1/2))$.

\begin{proposition}\label{diverge} The sequence $\big(\mu_{D,n}([0,1/2))\big)_{n\geqslant 0}$ does not converge. Moreover, it is not eventually periodic.
\end{proposition}

\begin{proof} To prove this result, we show that there is a subsequence of this sequence which is unbounded. To this end, note that \eqref{muDn} gives $$\mu_{D,n}([0,1/2))=\frac{1}{\varSigma_D(n)}\sum_{m=2^{n}}^{2^{n}+2^{n-1}-1}D(m).$$ Using \eqref{D} we have that \begin{multline*}\sum_{m=2^{n}}^{2^{n}+2^{n-1}-1}D(m)=e_1^T\Bigg(\sum_{\substack{w\in\{0,1\}^*\\ |w|=n-1}}{\bf B}_w\Bigg){\bf B}_0{\bf B}_1 e_1\\ =e_1^T\left({\bf B}_0+{\bf B}_1\right)^{n-1}{\bf B}_1 e_1=\sum_{m=2^{n-1}}^{2^{n}-1}D(m)=\varSigma_D(n-1),\end{multline*} where we have used that ${\bf B}_0$ is the identity matrix to obtain the middle equality. Thus $$\mu_{D,n}([0,1/2))=\frac{\varSigma_D(n-1)}{\varSigma_D(n)}.$$ By a simple calculation applying \eqref{SigmaDn} to both the numerator and denominator, with some rearrangement, we have $$\mu_{D,n}([0,1/2))=\frac{1}{4+3\mathrm{i}}\cdot\frac{1-\mathrm{e}^{-\mathrm{i}(2 \vartheta\cdot(n-1)-\pi/2)}}{1-\mathrm{e}^{-\mathrm{i}(2 \vartheta\cdot(n)-\pi/2)}},$$ where $\vartheta\approx 0.6435$ is the solution of $\cos\vartheta=4/5$. Now since $\vartheta$ is irrational and not a rational multiple of $\pi$, we have that the sequence of fractional parts $$\big(\{2 \vartheta\cdot(n)-\pi/2\}\big)_{n\geqslant 0}$$ is equidistributed in $[0,1)$. In particular, let $M>0$ be a positive integer. Then there is an $\varepsilon>0$ satisfying $$0<\varepsilon<\frac{\left|1-\mathrm{e}^{-2\mathrm{i} \vartheta}\right|}{5M+1},$$ and there are infinitely many $n$ such that $$\Big|1-\mathrm{e}^{-\mathrm{i}(2 \vartheta\cdot(n)-\pi/2)}\Big|<\varepsilon.$$ For these infinitely many $n$, $$\big|\mu_{D,n}([0,1/2))\big|>\frac{1}{5}\cdot\frac{\left|1-\mathrm{e}^{-2\mathrm{i} \vartheta}\right|-\varepsilon}{\varepsilon}=\frac{1}{5}\cdot\frac{\left|1-\mathrm{e}^{-2\mathrm{i} \vartheta}\right|}{\varepsilon}-\frac{1}{5}>M.$$ Since $M>0$ can be chosen arbitrarily large, $\big(\mu_{D,n}([0,1/2))\big)_{n\geqslant 0}$ is unbounded, which is the desired result.
\end{proof}

To observe large values of $\big(\mu_{D,n}([0,1/2))\big)_{n\geqslant 0}$ one must be patient. Figure \ref{fig:Dnoper} shows values of $\mu_{D,n}([0,1/2))$ for $n$ from $0$ to $100000$.
Modifying the proof of Proposition \ref{diverge}, \emph{mutatis mutandis}, one can show that the sequences $\big(\mu_{D,n}([0,1/2^k))\big)_{n\geqslant 0}$ are unbounded for every $k\geqslant 1$.

\begin{figure}[ht]
\begin{center}
  \includegraphics[width=.32\linewidth]{NoPeriod0to200.eps}
  \includegraphics[width=.32\linewidth]{NoPeriod201to10k.eps}
  \includegraphics[width=.32\linewidth]{NoPeriod10kplus1to100k.eps}\end{center}
\caption{The values of $\mu_{D,n}([0,1/2))$ for (left) $n=0,\ldots,200$, (middle) $n=201,\ldots,10000$ and (right) $n=10001,\ldots,100000$.}
\label{fig:Dnoper}
\end{figure}

\begin{corollary} The sequence of the pure point measures $\mu_{D,n}$ does not converge weakly to a finite Borel measure.
\end{corollary}

We now consider the dependence of Theorem \ref{BVexplicit} on the non-vanishing condition in Eq.~\eqref{zjcond}. While the dependence on this condition adds to the technicality of Theorem~\ref{BVexplicit} and to the difficulty in applying the result in full generality, seen from the view of some of common situations, the condition is quite natural. When ${\bf B}$'s dominant eigenvalue is simple, the non-vanishing condition in Eq.~\eqref{zjcond} is required to avoid certain degenerate situations. For example, consider the positive $k$-regular sequence $f$, where $f(m)=e_1^T\, {\bf B}_{(m)_k}\, {\bf w}$, with
$${\bf B}_0={\bf B}_1=\cdots={\bf B}_{k-1}=\left(\begin{matrix} * & 0\\ 0 & c \end{matrix}\right), \quad\mbox{and}\quad{\bf w}=\left(\begin{matrix} *\\ 0 \end{matrix}\right),$$ where the `$*$' indicates any positive values. For large integers $c$ the dominant eigenvalue of ${\bf B}$ will be $k\cdot c$, but none of the matrices ${\bf B}_a$ nor ${\bf w}$ have any component that interacts with this part of the matrix ${\bf B}$, and so this growth is not reflected in the behaviour of the sequence. Of course, the linear representation of the sequence is not unique, so this problem is not intrinsic to the sequence. More probable is that one has made a sub-optimal choice of the matrices ${\bf B}_a$---as is readily apparent, one can fix this example by deleting the un`$*$' parts.

With these examples in mind, assuming that degenerate situations can be avoided, the sufficient assumptions of Theorem \ref{BVexplicit} seem close to necessary for the existence of such measures. Of course, on a case-by-case basis, the measure could exist without satisfying these conditions, but it is certainly not guaranteed in general.

\section{Concluding remarks}\label{sec:cr}

In this paper, we have constructed measures associated to $k$-regular sequences and the vector spaces generated by their $k$-kernels in an effort to view these structures as a new source of dynamical objects, whose recursive properties are reminiscent of fractals.

In working with examples, curiosities  having to do with spectral type are evident, and natural examples with certain properties are elusive.

A case in point is the seeming scarcity of regular sequences yielding pure point measures. Of course, Theorem~\ref{main2} concludes that if $f$ is primitive (and other conditions), then $\mu_f$ is continuous. So to find examples, it seems wise to study sequences $f$ which have many zeros. For a trivial example, if one takes $f=\chi_2$ to be the characteristic sequence of the powers of two, then we easily compute that $\mu_{\chi_2}=\delta_0$. Staying with binary sequences for the moment, as soon as the density of ones in the fundamental regions $[2^n,2^{n+1})$ is bounded below, we are spreading the mass equally to these points, so the mass allocated to any one point in $\mathbb{T}$ must go to zero, and so assuming some regularity on the distribution of the non-zero values of the sequence the measure will be continuous; see Coons and Evans \cite{CEpre} for a family of singular continuous examples related to iterated function systems. Some non-trivial examples of sequences with pure-point measures exist, see Evans \cite{Epre}. However, these still occur in a relatively trivial way, where certain terms of the sequence contain a non-zero fraction of its entire mass. It seems likely that one could produce pure point measures differently, by instead having many small point masses concentrate in a single location, but it is unclear if this can be done in a non-trivial way within the confines of regular sequences. For a non-regular example, consider the values of the binomial coefficients, say $2^n\choose m$ between powers of two and follow our process. The measure obtained will be $\delta_{1/2}$. This happens because, even though the ratio of every value to the sum is going to zero, the sequence is dominated by the central binomial coefficients and the limiting Gaussian curve, when scaled to the interval $[0,1)$, is collapsing to its mean.

Our paradigm contrasts starkly with the situation for diffraction measures, where there is a plethora of examples with pure point spectral type, e.g. those arising from model sets (both in the regular and weak sense) and Toeplitz systems ; see Keller \cite{Keller} and Keller and Richard \cite{Keller3}.
For example, the diffraction measure of the paperfolding sequence is pure point \cite[p.~380]{BGbook}, but in our construction the produced measure is Lebesgue measure.
This is similar to the contrast between the spectral theory of Schr\"odinger operators and diffraction measures where absence of eigenvalues is also prevalent in the substitution setting; see Damanik, Embree and Gorodetski \cite{DEG2015}.

The question of spectral purity also arises. In most of the examples we have produced, the measure $\mu_f$ is of pure type. As suggested by Theorem~\ref{main2}, focussing solely on continuous measures, is there a natural example of a regular sequence with continuous measure of mixed spectral type, that is, having both absolutely continuous and singular continuous components in its Lebesgue decomposition? While examples can be constructed---using that the set of $k$-regular sequences is a group under point-wise addition \cite{AS1992}---we have yet to find a natural non-trivial example or a ``good'' sufficient criterion for purity.

As a final comment, we note that the H\"older exponent  of the distribution function $\mu_f([0,x])$ is related to (possibly fractal) geometric properties of the measure $\mu_f$. To formalise this, let $\mu$ be a probability measure on $\mathbb{T}$. The \emph{lower local dimension} of $\mu$ at $x\in\mathbb{T}$ is given by
\[
\underline{{\rm dim}}_{\rm loc} \, \mu(x)\colon=\liminf_{r\to 0} \frac{\log\mu(B_r(x))}{\log r}.
\]
One can also view $\underline{{\rm dim}}_{\rm loc} \, \mu(x)$ as the ``best possible local H\"older exponent,'' that is,
\[
\underline{{\rm dim}}_{\rm loc} \, \mu(x)=\sup \left\{\alpha\geqslant 0\colon \mu\big(B_r(x)\big)=O(r^{\alpha})\ \mbox{as}\ r\to 0\right\},
\]
see Bufetov and Solomyak \cite[Sec.~3.2]{BS2018}, Baake et al.~\cite{BGKS2019} and Mattila, Mor\'an and Rey~\cite{MMR2000}. If $\mu$ is absolutely continuous, one has $\underline{{\rm dim}}_{\rm loc} \, \mu(x)=1$ for Lebesgue almost all $x$, whereas $\underline{{\rm dim}}_{\rm loc} \, \mu(x)=0$ for all $x$ when $\mu$ is pure point.

\begin{proposition}
Let $\mu$ be a probability measure on $\mathbb{T}$, with distribution function $F$. Suppose $F$ is H\"older continuous with exponent $\alpha$. Then $\underline{{\rm dim}}_{\rm loc} \, \mu(x)\geqslant \alpha$ for all $x\in\mathbb{T}$.
\end{proposition}

\begin{proof}
One has, assuming $r<1$,
$$
\underline{{\rm dim}}_{\rm loc} \, \mu(x) =\liminf_{r\to 0}  \frac{\log\mu(B_r(x))}{\log r}= \liminf_{r\to 0}  \frac{\log|F(x+\frac{r}{2})-F(x-\frac{r}{2})|}{\log r}\geqslant \alpha,
$$
where the inequality follows from the H\"older continuity of $F$.
\end{proof}

We have the following immediate consequence of this for positive $k$-regular sequences.

\begin{corollary}
Let $f$ be a $k$-regular sequence which satisfies the conditions of Theorem~\ref{main2} and let $\mu_f$ be the corresponding continuous measure. Then, one has
\[
\underline{{\rm dim}}_{\rm loc} (\mu_f,x)\geqslant \log_{k}\big(\rho/\rho^{\ast}\big)
\]
uniformly on $\mathbb{T}$.
\end{corollary}

For example, considering the Josephus sequence, the distribution function is Lipschitz, so that $\underline{{\rm dim}}_{\rm loc}(\mu_J,x)\geqslant 1$ for all $x\in\mathbb{T}$. In particular, the local dimension exists and is one. Possibly a more interesting example is the Stern sequence \cite{BC2018}. Here one has a H\"older exponent strictly less than one. In this case    $\underline{{\rm dim}}_{\rm loc}(\mu_s,x)\geqslant \log_2({3}/{\tau})\approx 0.890721$ for all $x\in\mathbb{T}$ where $\tau$ is the golden mean.

\section*{Acknowledgements}

It is our pleasure to thank Michael Baake for helpful conversations and comments, and for suggesting important references. M.~Coons and J.~Evans acknowledge the support of the Commonwealth of Australia and N.~Ma\~nibo acknowledges the support of the German Research Foundation (DFG) via the CRC 1283.

\bibliographystyle{amsplain}
\providecommand{\bysame}{\leavevmode\hbox to3em{\hrulefill}\thinspace}
\providecommand{\MR}{\relax\ifhmode\unskip\space\fi MR }
\providecommand{\MRhref}[2]{%
  \href{http://www.ams.org/mathscinet-getitem?mr=#1}{#2}
}
\providecommand{\href}[2]{#2}