We define a compact formalism for multiunit spike-train metrics by
opportunely characterising the space of multiunit feature vectors as a
tensor product. Previous results from Houghton and Kreuz (2012, On the
efficient calculation of van Rossum distances. Network: Computation in
Neural Systems, 2012, 23, 48-58) on a clever formula for multiunit Van
Rossum metrics are then re-derived within this framework, also fixing
some errors in the original calculations.
A compact formalism for kernel-based multiunit spike train metrics
Consider a network with \(C\) cells. Let
\[\mathcal{U}= \left\{ \boldsymbol{u}^1, \boldsymbol{u}^2, \ldots, \boldsymbol{u}^C \right \}\]
be an observation of network activity, where
\[\boldsymbol{u}^i = \left\{ u_1^i, u_2^i, \ldots, u_{N_{\boldsymbol{u}^i}}^i \right \}\]
is the (ordered) set of times of the spikes emitted by cell \(i\).
Let
\(\mathcal{V}= \left\{\boldsymbol{v}^1, \boldsymbol{v}^2, \ldots, \boldsymbol{v}^C\right\}\)
be another observation, different in general from \(\mathcal{U}\).
To compute a kernel based multiunit distance between \(\mathcal{U}\)
and \(\mathcal{V}\), we map them to the tensor product space
\(\mathcal{S} \doteq
\mathbb{R}^C\bigotimes L_2(\mathbb{R}\rightarrow\mathbb{R})\) by defining
\[|\mathcal{U}\rangle = \sum_{i=1}^C |i\rangle \otimes |\boldsymbol{u}^i\rangle\]
where we consider \(\mathbb{R}^C\) and
\(L_2(\mathbb{R}\rightarrow\mathbb{R})\) to be equipped with the
usual euclidean distances, consequently inducing an euclidean metric
structure on \(\mathcal{S}\) too.
Conceptually, the set of vectors
\(\left\{ |i\rangle \right\}_{i=1}^C
\subset \mathbb{R}^C\) represents the different cells, while each
\(|\boldsymbol{u}^i\rangle \in L_2(\mathbb{R}\rightarrow\mathbb{R})\)
represents the convolution of a spike train of cell \(i\) with a
real-valued feature function
\(\phi: \mathbb{R}\rightarrow\mathbb{R}\),
\[\langle t|\boldsymbol{u}\rangle = \sum_{n=1}^{N_{\boldsymbol{u}}}\phi(t-u_n)\]
In practice, we will never use the feature functions directly, but we
will be only interested in the inner products of the \(|i\rangle\)
and \(|\boldsymbol{u}\rangle\) vectors. We call
\(c_{ij}\doteq\langle i|j \rangle=\langle i|j \rangle_{\mathbb{R}^C}=c_{ji}\)
the multiunit mixing coefficient for cells \(i\) and \(j\),
and
\(\langle \boldsymbol{u}|\boldsymbol{v}\rangle=\langle \boldsymbol{u}|\boldsymbol{v}\rangle_{L_2}\)
the single-unit inner product,
\[\begin{split}\label{eq:singleunit_intprod}
\begin{split}
\langle \boldsymbol{u}|\boldsymbol{v}\rangle & = \langle \left\{ u_1, u_2, \ldots,
u_{N}\right\}|\left\{ v_1, v_2, \ldots, v_{M}\right\} \rangle = \\
&= \int\textrm{d }\!t \langle \boldsymbol{u}|t \rangle\langle t|\boldsymbol{v}\rangle = \int\textrm{d }\!t
\sum_{n=1}^N\sum_{m=1}^M\phi\left(t-u_n\right)\phi\left(t-v_m\right)\\
&\doteq \sum_{n=1}^N\sum_{m=1}^M \mathcal{K}(u_n,v_m)
\end{split}\end{split}\]
where \(\mathcal{K}(t_1,t_2)\doteq\int\textrm{d }\!t
\left[\phi\left(t-t_1\right)\phi\left(t-t_2\right)\right]\) is the
single-unit metric kernel, and where we have used the fact that the
feature function \(\phi\) is real-valued. It follows immediately
from the definition above that
\(\langle \boldsymbol{u}|\boldsymbol{v}\rangle=\langle \boldsymbol{v}|\boldsymbol{u}\rangle\).
Note that, given a cell pair \((i,j)\) or a spike train pair
\((\boldsymbol{u},\boldsymbol{v})\), \(c_{ij}\) does not depend
on spike times and \(\langle \boldsymbol{u}|\boldsymbol{v}\rangle\)
does not depend on cell labeling.
With this notation, we can define the multi-unit spike train distance
as
\[\label{eq:distance_as_intprod}
\left\Vert |\mathcal{U}\rangle - |\mathcal{V}\rangle \right\Vert^2 = \langle \mathcal{U}|\mathcal{U}\rangle + \langle \mathcal{V}|\mathcal{V}\rangle
- 2 \langle \mathcal{U}|\mathcal{V}\rangle\]
where the multi-unit spike train inner product
\(\langle \mathcal{V}|\mathcal{U}\rangle\) between
\(\mathcal{U}\) and \(\mathcal{V}\) is just the natural bilinear
operation induced on \(\mathcal{S}\) by the tensor product
structure:
\[\begin{split}\label{eq:multiunit_intprod}
\begin{split}
\langle \mathcal{V}|\mathcal{U}\rangle &= \sum_{i,j=1}^C
\langle i|j \rangle\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle
= \sum_{i,j=1}^C c_{ij}\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle\\
&= \sum_{i=1}^C\left[ c_{ii}
\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + c_{ij}
\left(\sum_{j<i}\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle +
\sum_{j>i}\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle \right) \right]
\end{split}\end{split}\]
But \(c_{ij}=c_{ji}\) and
\(\langle \boldsymbol{v}|\boldsymbol{u}\rangle=\langle \boldsymbol{u}|\boldsymbol{v}\rangle\),
so
\[\begin{split}\begin{split}
\sum_{i=1}^C\sum_{j<i}c_{ij}\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle &=
\sum_{j=1}^C\sum_{i<j}c_{ji}\langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle =
\sum_{i=1}^C\sum_{j>i}c_{ji}\langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle =
\sum_{i=1}^C\sum_{j>i}c_{ij}\langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle\\
&=\sum_{i=1}^C\sum_{j>i}c_{ij}\langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle
\end{split}\end{split}\]
and
\[\begin{split}\label{eq:multiprod_j_ge_i}
\langle \mathcal{V}|\mathcal{U}\rangle = \sum_{i=1}^C\left[ c_{ii}
\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + c_{ij}
\sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle +
\langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle \right) \right]\end{split}\]
Now, normally we are interested in the particular case where
\(c_{ij}\) is the same for all pair of distinct cells:
\[\begin{split}c_{ij} = \begin{cases}
1 & \textrm{if } i=j\\
c & \textrm{if } i\neq j
\end{cases}\end{split}\]
and under this assumption we can write
\[\begin{split}\label{eq:multiprod_constant_c}
\langle \mathcal{V}|\mathcal{U}\rangle = \sum_{i=1}^C\left[\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + c\sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle \right) \right]\end{split}\]
and
\[\begin{split}\begin{split}
\left\Vert |\mathcal{U}\rangle - |\mathcal{V}\rangle \right\Vert^2 =
\sum_{i=1}^C\Bigg\{&\langle \boldsymbol{u}^i|\boldsymbol{u}^i \rangle +
c\sum_{j>i}\left(\langle \boldsymbol{u}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{u}^i|\boldsymbol{u}^j \rangle
\right) +\\
+&\langle \boldsymbol{v}^i|\boldsymbol{v}^i \rangle + c\sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{v}^j \rangle +
\langle \boldsymbol{v}^i|\boldsymbol{v}^j \rangle \right) + \\
-2&\left[\langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle +
c\sum_{j>i}\left(\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{u}^i|\boldsymbol{v}^j \rangle
\right)\right] \Bigg\}
\end{split}\end{split}\]
Rearranging the terms
\[\begin{split}\begin{gathered}
\label{eq:multidist_constant_c}
\left\Vert |\mathcal{U}\rangle - |\mathcal{V}\rangle \right\Vert^2 =
\sum_{i=1}^C\Bigg[\langle \boldsymbol{u}^i|\boldsymbol{u}^i \rangle + \langle \boldsymbol{v}^i|\boldsymbol{v}^i \rangle
- 2 \langle \boldsymbol{v}^i|\boldsymbol{u}^i \rangle + \\
+ 2c\sum_{j>i}\left(\langle \boldsymbol{u}^i|\boldsymbol{u}^j \rangle + \langle \boldsymbol{v}^i|\boldsymbol{v}^j \rangle
-\langle \boldsymbol{v}^i|\boldsymbol{u}^j \rangle - \langle \boldsymbol{v}^j|\boldsymbol{u}^i \rangle\right)\Bigg]\end{gathered}\end{split}\]
Van Rossum-like metrics
In Van Rossum-like metrics, the feature function and the single-unit
kernel are, for \(\tau\neq 0\),
\[\begin{split}\begin{gathered}
\phi^{\textrm{VR}}_{\tau}(t) = \sqrt{\frac{2}{\tau}}\cdot e^{-t/\tau}\theta(t) \\
\mathcal{K}^{\textrm{VR}}_{\tau}(t_1,t_2) = \begin{cases}
1 & \textrm{if } t_1=t_2\\
e^{-\left\vert t_1-t_2 \right\vert/\tau} & \textrm{if } t_1\neq t_2
\end{cases}\end{gathered}\end{split}\]
where \(\theta\) is the Heaviside step function (with
\(\theta(0)=1\)), and we have chosen to normalise
\(\phi^{\textrm{VR}}_{\tau}\) so that
\[\left\Vert \phi^{\textrm{VR}}_{\tau} \right\Vert_2 = \sqrt{\int\textrm{d }\!t \left[\phi^{\textrm{VR}}_{\tau}(t)\right]^2} = 1 \quad.\]
In the \(\tau\rightarrow 0\) limit,
\[\begin{split}\begin{gathered}
\phi^{\textrm{VR}}_{0}(t) = \delta(t)\\
\mathcal{K}^{\textrm{VR}}_{0}(t_1,t_2) = \begin{cases}
1 & \textrm{if } t_1=t_2\\
0 & \textrm{if } t_1\neq t_2
\end{cases}\end{gathered}\end{split}\]
In particular, the single-unit inner product now becomes
\[\langle \boldsymbol{u}|\boldsymbol{v}\rangle = \sum_{n=1}^N\sum_{m=1}^M
\mathcal{K}^{\textrm{VR}}(u_n,v_m) = \sum_{n=1}^N\sum_{m=1}^M
e^{-\left\vert u_n-v_m \right\vert/\tau}\]
Markage formulas
For a spike train \(\boldsymbol{u}\) of length \(N\) and a time
\(t\) we define the index
\(\tilde{N}\left( \boldsymbol{u}, t\right)\)
\[\begin{split}\tilde{N}\left( \boldsymbol{u}, t\right) \doteq \max\{n | u_n < t\}\end{split}\]
which we can use to re-write
\(\langle \boldsymbol{u}|\boldsymbol{u}\rangle\) without the
absolute values:
\[\begin{split}\begin{split}
\langle \boldsymbol{u}|\boldsymbol{u}\rangle&= \sum_{n=1}^N \left(
\sum_{m|v_m<u_n}e^{-(u_n-v_m)/\tau} +
\sum_{m|v_m>u_n}e^{-(v_m-u_n)/\tau} +
\sum_{m=1}^M\delta\left(u_n,v_m\right)
\right)\\
&= \sum_{n=1}^N \left( \sum_{m|v_m<u_n}e^{-(u_n-v_m)/\tau} +
\sum_{m|u_m<v_n}e^{-(v_n-u_m)/\tau} +
\sum_{m=1}^M\delta\left(u_n,v_m\right) \right)\\
&= \sum_{n=1}^N \left[
\sum_{m=1}^{\tilde{N}\left( \boldsymbol{v}, u_n\right)}e^{-(u_n-v_m)/\tau} +
\sum_{m=1}^{\tilde{N}\left( \boldsymbol{u}, v_n\right)}e^{-(v_n-u_m)/\tau} +
\delta\left(u_n,v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)+1}\right)
\right]\\
&= \sum_{n=1}^N \Bigg[ e^{-(u_n-v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)})/\tau}
\sum_{m=1}^{\tilde{N}\left( \boldsymbol{v}, u_n\right)}e^{-(v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)}-v_m)/\tau} +
\\
&\phantom{ = \sum_{n=1}^N } + e^{-(v_n-u_{\tilde{N}\left( \boldsymbol{u}, v_n\right)})/\tau}
\sum_{m=1}^{\tilde{N}\left( \boldsymbol{u}, v_n\right)}e^{-(u_{\tilde{N}\left( \boldsymbol{u}, v_n\right)}-u_m)/\tau}
+ \\
&\phantom{ = \sum_{n=1}^N } + \delta\left(u_n,v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)+1}\right) \Bigg]\\
\end{split}\end{split}\]
For a spike train \(\boldsymbol{u}\) of length \(N\), we also
define the the markage vector \(\boldsymbol{m}\), with the same
length as \(\boldsymbol{u}\), through the following recursive
assignment:
\[\begin{split}\begin{aligned}
m_1(\boldsymbol{u}) &\doteq 0 \\
m_n(\boldsymbol{u}) &\doteq \left(m_{n-1} + 1\right) e^{-(u_n - u_{n-1})/\tau}
\quad \forall n \in \{2,\ldots,N\}\label{eq:markage_definition}\end{aligned}\end{split}\]
It is easy to see that
\[\begin{split}\begin{split}
m_n(\boldsymbol{u}) &= \sum_{k=1}^{n-1}e^{-(u_n - u_k)/\tau} =
\left(\sum_{k=1}^{n}e^{-(u_n - u_k)/\tau}\right) - e^{-(u_n - u_n)/\tau}\\
&= \sum_{k=1}^{n}e^{-(u_n - u_k)/\tau} - 1
\end{split}\end{split}\]
and in particular
\[\label{eq:markage_sum}
\sum_{n=1}^{\tilde{N}\left( \boldsymbol{u}, t\right)}e^{-(u_{\tilde{N}\left( \boldsymbol{u}, t\right)}-u_n)/\tau} = 1 + m_{\tilde{N}\left( \boldsymbol{u}, t\right)}(\boldsymbol{u})\]
With this definition, we get
\[\begin{split}\label{eq:singleunit_intprod_markage}
\begin{split}
\langle \boldsymbol{u}|\boldsymbol{v}\rangle &= \sum_{n=1}^N \Bigg[
e^{-(u_n-v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)})/\tau} \left(1 + m_{\tilde{N}\left( \boldsymbol{v}, u_n\right)}(\boldsymbol{v})\right)
+ \\
&\phantom{= \sum_{n=1}^N} + e^{-(v_n-u_{\tilde{N}\left( \boldsymbol{u}, v_n\right)})/\tau}
\left(1 + m_{\tilde{N}\left( \boldsymbol{u}, v_n\right)}(\boldsymbol{u})\right) + \\
&\phantom{= \sum_{n=1}^N} +
\delta\left(u_n,v_{\tilde{N}\left( \boldsymbol{v}, u_n\right)+1}\right) \Bigg]
\end{split}\end{split}\]
Finally, note that because of the definition of the markage vector
\[e^{-(u_n-u_{\tilde{N}\left( \boldsymbol{u}, u_n\right)})/\tau} \left(1 +
m_{\tilde{N}\left( \boldsymbol{u}, u_n\right)}(\boldsymbol{u})\right) = e^{-(u_n-u_{n-1})/\tau}\left(1+m(\boldsymbol{u})\right) = m_{n}(\boldsymbol{u})\]
so that in particular
\[\begin{split}\label{eq:singleunit_squarenorm_markage}
\begin{split}
\langle \boldsymbol{u}|\boldsymbol{u}\rangle &= \sum_{n=1}^N \left(1 + 2m_{n}(\boldsymbol{u})\right)
\end{split}\end{split}\]
A formula for the efficient computation of multiunit Van Rossum spike
train metrics, used by pymuvr
, can then be obtained by opportunely
substituting these expressions for the single-unit scalar products in
the definition of the multiunit distance.