A Compact Orbital Angular Momentum Spectrometer Using Quantum Zeno
Interrogation

Paul Bierdz, Hui Deng*

Department of Physics, University of Michigan, Ann Arbor, USA

<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a scheme to measure the orbital angular momentum spectrum of light
using a precisely timed optical loop and quantum non-demolition measurements.
We also discuss the influence of imperfect optical components.

OCIS codes: 050.4865 Optical vortices, 120.4290 Nondestructive testing.

## References and links

  * [1] L. Allen, M. W. Beijersbergen, R. J. C. Spreeuw, and J. P. Woerdman, “Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes,” Physical Review A 45, 8185 (1992).
  * [2] G. Molina-Terriza, J. P. Torres, and L. Torner, “Twisted photons,” Nat Phys 3, 305–310 (2007).
  * [3] S. Franke-Arnold, L. Allen, and M. J. Padgett, “Advances in optical angular momentum,” Laser & Photonics Review 2, 299–313 (2008).
  * [4] G. Molina-Terriza, J. P. Torres, and L. Torner, “Management of the angular momentum of light: Preparation of photons in multidimensional vector states of angular momentum,” Physical Review Letters 88, 013,601 (2001).
  * [5] A. Mair, A. Vaziri, G. Weihs, and A. Zeilinger, “Entanglement of the orbital angular momentum states of photons,” Nature 412, 313–316 (2001).
  * [6] G. Molina-Terriza, A. Vaziri, R. Ursin, and A. Zeilinger, “Experimental quantum coin tossing,” Physical Review Letters 94, 040,501 (2005).
  * [7] D. Kaszlikowski, P. Gnacinacuteski, M. Zdotukowski, W. Miklaszewski, and A. Zeilinger, “Violations of local realism by two entangled N-Dimensional systems are stronger than for two qubits,” Physical Review Letters 85, 4418 (2000).
  * [8] N. J. Cerf, M. Bourennane, A. Karlsson, and N. Gisin, “Security of quantum key distribution using d-Level systems,” Physical Review Letters 88, 127,902 (2002).
  * [9] B. P. Lanyon, M. Barbieri, M. P. Almeida, T. Jennewein, T. C. Ralph, K. J. Resch, G. J. Pryde, J. L. O’Brien, A. Gilchrist, and A. G. White, “Simplifying quantum logic using higher-dimensional hilbert spaces,” Nat Phys 5, 134–140 (2009).
  * [10] J. T. Barreiro, T. Wei, and P. G. Kwiat, “Beating the channel capacity limit for linear photonic superdense coding,” Nat Phys 4, 282–286 (2008).
  * [11] J. M. Hickmann, E. J. S. Fonseca, W. C. Soares, and S. Chávez-Cerda, “Unveiling a truncated optical lattice associated with a triangular aperture using light’s orbital angular momentum,” Physical Review Letters 105, 053,904 (2010).
  * [12] Z. Wang, Z. Zhang, and Q. Lin, “A novel method to determine the helical phase structure of Laguerre–Gaussian beams,” Journal of Optics A: Pure and Applied Optics 11, 085,702 (2009).
  * [13] V. Bazhekov, M. Vasnetsov, and M. S. Soskin, “Laser-beams with screw dislocation in their wave-fronts,” JETP Letters 52, 429–431 (1990).
  * [14] N. R. Heckenberg, R. McDuff, C. P. Smith, and A. G. White, “Generation of optical phase singularities by computer-generated holograms,” Optics Letters 17, 221–223 (1992).
  * [15] J. Arlt, K. Dholakia, L. Allen, and M. J. Padgett, “The production of multiringed Laguerre–Gaussian modes by computer-generated holograms,” Journal of Modern Optics 45, 1231 (1998).
  * [16] J. Leach, J. Courtial, K. Skeldon, S. M. Barnett, S. Franke-Arnold, and M. J. Padgett, “Interferometric methods to measure orbital and spin, or the total angular momentum of a single photon,” Physical Review Letters 92, 013,601 (2004).
  * [17] G. C. G. Berkhout, M. P. J. Lavery, J. Courtial, M. W. Beijersbergen, and M. J. Padgett, “Efficient sorting of orbital angular momentum states of light,” Physical Review Letters 105, 153,601 (2010).
  * [18] A. Peres, “Zeno paradox in quantum theory,” American Journal of Physics 48, 931 (1980).
  * [19] A. C. Elitzur and L. Vaidman, “Quantum mechanical interaction-free measurements,” Foundations of Physics 7, 987–997 (1993).
  * [20] P. G. Kwiat, A. G. White, J. R. Mitchell, O. Nairz, G. Weihs, H. Weinfurter, and A. Zeilinger, “High-Efficiency quantum interrogation measurements via the quantum Zeno effect,” Physical Review Letters 83, 4725 (1999).
  * [21] S0 and S1 are optical switches that can be switched from been transmittive to reflective. S0 needs to transmit the initial incident light, and be switched to be reflective by the end of the first outer-loop cycle. A high repetition rate is not required and it can be implemeted either mechanically or opto-electrically. Alternatively it could also be a static high-reflectance mirror, if the one-time transmission loss at the very beginning can be tolerated. S1 needs to be switched every QZI loop cycle ($\Delta T$) and need to be polarization insensitive. The fastest switches are Pockels cells, which can operate at 10-100 GHz with $99\%$ transmission. One scheme, similar to the one implemented in Ref. [20], is to an interferometer with Pockels cells in one of the arms. The Pockels cells introduce $\pi$ phase shift in the arm when activated, and thus switch the beam between two output ports of the interferometer. Using two Pockels cells rotated relative to each other will cancel the birefringent effect.
  * [22] M. W. Beijersbergen, R. P. C. Coerwinkel, M. Kristensen, and J. P. Woerdman, “Helical-wavefront laser beams produced with a spiral phaseplate,” Optics Communications 112, 321–327 (1994).
  * [23] B. Misra and E. C. G. Sudarshan, “The Zeno’s paradox in quantum theory,” Journal of Mathematical Physics 18, 756 (1977).
  * [24] Technically, each $|\alpha|^{2}$ is slightly different, but the difference is well within 1%. The $|\alpha|^{2}$ that matters most is the one corresponding to the QZI loop (including S1), which, without any optimization, consists of 4 beam splitters, 5 mirrors, 1 waveplate and up to three Pockels cells. Assuming all optics are anti-reflection coated so that loss is 1% at each Pockels cell and 0.1% at each other component, we have $|\alpha|^{2}\geq 0.96$.
  * [25] J. Jang, “Optical interaction-free measurement of semitransparent objects,” Physical Review A 59, 2322 (1999).

Light carries quantized orbital angular momentum (OAM) of $l\hbar$ per photon
when the electric field has an overall azimuthal dependence on phase,
$e^{-il\phi}$, where $\phi$ is the azimuthal angle about the beam propagation
axis [1, 2, 3]. The OAM quantum number $l$ takes integer values from $-\infty$
to $+\infty$. With an infinite number of states available, OAM can be utilized
for qudit ($d$-dimensional, $d>2$) systems [4] that allow, for example, higher
dimensional entanglement [5], quantum coin-tossing [6], increased violations
of local realism [7], improved security for quantum key distribution [8],
simplified quantum gates [9] and superdense coding for quantum communications
with increased channel capacity [10].

Determining OAM of light in an unknown state, however, is more challenging
than measuring different polarizations or frequencies. Eigenstates of OAM can
be deduced from the diffraction interference pattern with judiciously chosen
apertures [11, 12]. But the method requires a large collection of photons to
develop the pattern. Moreover it may become very complicated for
superpositions of OAM states. An $l$-fold fork diffraction grating [13, 14,
15] can separate a pre-determined OAM component from others at the single
photon level, but cannot be readily applied to determine arbitrary OAM state
of light. A cascade Mach-Zehnder interferometer setup, using Dove prisms to
introduce an $l$-dependent phase shift, can separate different OAM components
of light into different output ports of the interferometers, even at the
single-photon level [16]. But the setup requires $N-1$ mutually stabilized
interferometers to detect $N$ OAM-modes. A recently proposed scheme [17]
requires only a spatial light modulator (or a specially designed hologram),
which converts the twisting phase structure of OAM states into a linear phase
gradient, and a lens, which focuses different OAM components to different
spatial locations on the focal plane. However, the method relies on intricate
spatial modulation of the phase, and thus may have limited applicability to
broadband ultrafast pulses. Moreover, the extinction ratio between different
OAM states is limited to about $10$. To increase the extinction ratio or to
detect higher order OAM states will require larger and increasingly complex
holograms.

In this paper, we present a compact OAM-spectrometer comprising of only one
interferometer nested within an optical loop (Fig. 1). It uses a Quantum Zeno
Interrogator (QZI) [18, 19, 20] (shaded region in Fig. 1) to perform
counterfactual measurements on the OAM state, and thus maps different OAM
components of an arbitrary input light pulse into different time bins at the
output. It can achieve very high extinction ratios between different OAM
states and can work for arbitrarily high OAM orders limited mainly by optical
losses.

Fig. 1: A schematic of the compact OAM spectrometer. The Quantum Zeno
Interrogator (shaded region) distinguishes between zero and nonzero OAM
states. The outer loop decreases the OAM value of light by one per round trip.
All the beam splitters are polarizing beam splitters (PBSs) that transmits
horizontally polarized light and reflects vertically polarized light. The OAM
filter transmits states with zero OAM, but blocks states with non-zero OAM. S0
and S1 are switching mirrors that either transmits or reflects incident light
[21]. R1 and R2 are fixed polarization rotators, which can be half wave
plates. P1 and P2 are fast polarization switches, such as Pockels cells. When
activated, P1 and P2 switches horizontal polarization to vertical and vice
versa. When de-activated, they are transparent to light. The shaded region is
a Quantum Zeno Interrogator [20] which separates OAM components with $l=0$ and
$l\neq 0$ into different polarizations. Hence at PBS3, zero OAM component is
sent to the detector while the none-zero OAM component is sent back into the
outer-loop. The outer loop decreased OAM by one per round trip via, for
example, a vortex phase plate (VPP) [22].

We illustrate now how the spectrometer works by tracing, as an example, a
horizontally polarized input pulse with an OAM value $l=l_{0}\geq 0$, noted as
$|\psi^{(0)}\rangle=|H,l_{0}\rangle$. The input pulse first transmits through
optical switches S0 and S1 [21], and enters the QZI. The polarization rotator
R1 rotates its polarization by $\Delta\theta=\pi/(2N)$, and the state becomes
$|\psi^{(0)}_{1}\rangle=\cos\left(\frac{\pi}{2N}\right)|H,l_{0}\rangle+\sin\left(\frac{\pi}{2N}\right)|V,l_{0}\rangle$.
If $l=0$, the horizontal and vertical components of $|\psi^{(0)}_{1}\rangle$
passes through the lower and upper arms of the interferometer, respectively.
They recombine into the same state $|\psi^{(0)}_{1}\rangle$ at the polarizing
beam splitter PBS2 (neglecting an overall phase factor). S1 is switched to be
reflective at the end of the first QZI loop, and the combined beam continues
to loop in the QZI. The polarization is rotated by $\Delta\theta=\pi/(2N)$
each loop. After $N$ loops, the light becomes vertically polarized and enters
only the upper path of the interferometer. At this point, the polarization
switch P1 is activated and switches the polarization into horizontal. Hence
the light transmits through both PBS2 and PBS3, and arrives at the detector at
time $T_{0}$.

If $l_{0}\neq 0$, however, the vertical component is sent to the upper path at
PBS1, and is then blocked by the OAM filter. Only the horizontal component
emerges after PBS2, the state collapses into $|H,l_{0}\rangle$ with a
probability $\cos^{2}\left(\frac{\pi}{2N}\right)$. After $N$ loops, a fraction
$p=\cos(\pi/(2N))^{2N}$ of the light remains in the horizontal polarization in
the lower arm of the interferometer, while a fraction $1-p$ of the light is
lost (blocked by OAM filter). At this point, the polarization switch P2 is
activated and switches the polarization to vertical, and the light reflects
off both PBS2 and PBS3, and enters the outer loop. By this time, S0 is
switched to be reflective. As the light cycles in the outer loop, the
polarization in rotated back to horizontal by R2, and the OAM value is
decreased by $\Delta l=1$ per cycle by a vortex phase plate (VPP) [22]. After
$l_{0}$ cycles, $l=0$. When the light enters the QZI again, it will exit the
spectrometer to the detector, at a time
$T(l_{0})=T_{0}+l_{0}(NL_{QZI}+L_{out})/c$. Here $L_{QZI}$ and $L_{out}$ are
the optical path lengths of the QZI loop (from S1 to PBS2 back to S1) and the
outer-loop (from S1 to PBS2, to PBS3, to S0, back to S1). The detected
fraction of the light intensity is
$P(l_{0})=p^{l_{0}}=\cos\left(\frac{\pi}{2N}\right)^{2Nl_{0}}$.

In short, the OAM spectrometer sorts different OAM components into different
time intervals separated by $\Delta T=(NL_{QZI}+L_{out})/c$ with a _perfect_
extinction ratio. The total transmission efficiency of the spectrometer is
$P(l_{0})$ for the component with OAM of $l_{0}\hbar$. $P(l_{0})\rightarrow 1$
for all $l_{0}$ as $N\rightarrow\infty$ due to the quantum Zeno effect [23],
as shown in the Fig. 2(a).

Fig. 2: The probability of detecting the correct OAM value as a function of
the number of loops ($N$) in the QZI using a perfect OAM filter. (a) Neglect
optical loss. (b) Assume $|\alpha|^{2}=0.96$ based on commercially available
optics. When optical loss is included, there exists an optimal $N$ for higher
order OAM states, due to the compromise between the quantum Zeno enhancement
and optical loss.

In practice, optical components introduce loss. Assuming high quality, but
commercially available optical components, we estimate a round trip
transmission of $|\alpha|^{2}\sim 0.96$ [24] per cycle for both the outer loop
($\alpha_{out}$), the QZI loop ($\alpha_{QZI}$) and initial and final optics
($\alpha_{init,final}$). Hence the total transmission efficiency of the OAM
spectrometer becomes
$P(l_{0})\approx\alpha^{(2N+2)(l_{0}+1)}\cos\left(\frac{\pi}{2N}\right)^{2Nl_{0}}$
for the $l_{0}$-th order OAM component. We plot in Fig. 2(b) the $P(l_{0})$
vs. $N$ for OAM components $l_{0}=0-10$. With increasing $N$, the quantum Zeno
effect leads to an increase in $P(l_{0})$, while loss leads to a decrease in
$P(l_{0})$. As a result, an optimal $N$ is found at about $7-8$ for high order
OAM components. Note that the extinction ratio between different OAM states
remains infinite even in the presence of loss. Crosstalk would only take place
when the OAM filter is not completely opaque to nonzero OAM states.

To take into account imperfect OAM filters, we derive below the general
expression for the transmission efficiency and extinction ratio, with finite
$N$ and optical loss. We consider the OAM filter having a complex transmission
coefficient $\sqrt{T(l)}e^{i\phi(l)}$ for the $l$th OAM component. If the
state $|\psi\rangle=|H,l_{0}\rangle$ enters the QZI, after $N$ cycles, it
exits the QZI loop in a polarization superposition state
$p_{H}|H\rangle+p_{V}|V\rangle$ [25], where $p_{H}$ and $p_{V}$ are given by:

$\left(\begin{array}[]{c}p_{H}\\\ p_{V}\\\
\end{array}\right)=\alpha_{QZI}^{N}\left[\left(\begin{array}[]{cc}1&0\\\
0&\sqrt{T(l)}e^{i\phi(l)}\\\
\end{array}\right)\left(\begin{array}[]{cc}\cos\left(\frac{\pi}{2N}\right)&\sin\left(\frac{\pi}{2N}\right)\\\
\sin\left(\frac{\pi}{2N}\right)&-\cos\left(\frac{\pi}{2N}\right)\\\
\end{array}\right)\right]^{N}\left(\begin{array}[]{c}1\\\ 0\\\
\end{array}\right).$ (1)

The pulse re-enters the outer loop at PBS3 with probability $|p_{H}|^{2}$,
corresponding to a successful interrogation by the QZI (if $l_{0}\neq 0$).
With probably $|p_{V}|^{2}$, the pulse exits toward the detector,
corresponding to an error (if $l_{0}\neq 0$). The total loss of this QZI
interrogation is $|loss|^{2}=1-|p_{H}|^{2}-|p_{V}|^{2}$. In the outer loop,
the OAM value of the pulse is lowered by $1$ via the VPP, and the intensity of
the pulse is reduced by a factor $|\alpha_{out}|^{2}$ per loop. Therefore, the
probability of detecting the OAM eigenstate $l_{0}$ in the $l$th time interval
(or, measured as with OAM $l\hbar$) is given by:

$P(l;l_{0})=|\alpha_{init,final}|^{2}|p_{V}(l_{0}-l)|^{2}\prod_{m=l_{0}-l+1}^{l_{0}}\left(|\alpha_{out}|^{2}|p_{H}(m)|^{2}\right).$
(2)

And we define the extinction ratio $\eta$ as:

$\eta(l_{0})=P(l_{0};l_{0})/\sum_{l\neq l_{0}}P(l;l_{0}).$ (3)

With an imperfect OAM filter, light with nonzero OAM has a finite probability
of transmitting through the filter in vertical polarization after the $N$th
QZI-loop. It will then be switched to horizontal polarization by P1 and exit
at a time interval corresponding to components with a lower OAM. Consequently,
the extinction ratio is reduced. If the light is transmitted through the
filter before the $N$th loop, it will results in a larger loss. An imperfect
OAM filter may also partially block light with zero OAM, which which also
results in loss.

Figure 3(a) shows, per quantum Zeno interrogation of light with OAM of
$l\hbar$, the probabilities of the light exiting toward the detector
($|p_{V}|^{2}$), re-entering the outer-loop ($|p_{H}|^{2}$) and being lost
($1-|p_{V}|^{2}-|p_{H}|^{2}$). These probabilities are plotted as a function
of transmission $T(l,a_{0})$ and $N$. The crossing of $|p_{H}|^{2}$ and
$|p_{V}|^{2}$ separates the regimes when the interrogation result is more
likely (to the right side) or less likely (to the left side) to be correct
than incorrect.

Fig. 3: The probabilities of different outcomes of a QZI interrogation as a
function of the transmission of the OAM filter, neglecting optical loss. The
blue solid line represents detecting OAM=0, the red dashed line is detecting
OAM$\neq 0$, and the orange dotted line, loss. (a) $N=8$. (b) $N=2-10$. Fig.
4: Transmission of the pinhole spatial filter (a) as a function of the
normalized aperture size $a_{0}$, for OAM components with $l_{0}=0-3$ and (b)
as a function of $l_{0}$ with $a_{0}=0.8$. Fig. 5: (a) Extinction ratio $\eta$
as a function of the number of loops $N$ for various losses $|\alpha|^{2}$.
Solid symbols are for $l_{0}=1$ and open symbols are for $l_{0}=3$. $l_{0}>3$
are essentially indistinguishable from $l_{0}=3$. For the $l_{0}=0$ case, the
extinction ratio is over a 1000 for all $|\alpha|^{2}$ values because no
premature measurements are possible. The additional green crosses labeled as
$|\alpha|^{2}=0.95^{*}$ represents $|\alpha|^{2}=0.96$ but including
misalignment of the OAM filter and VPP as discussed in the text. (b)
Extinction ratio $\eta$ as a function of the normalized aperture size $a_{0}$
for $l_{0}=6$, $\Delta l=1-3$, $N=8$, and $|\alpha|^{2}=0.96$. Skipping OAM
states increases the extinction ratio by orders of magnitude.

As a practical example of an imperfect OAM filter, we consider a pinhole
spatial filter. Light with OAM of $l\hbar\neq 0$ has zero intensity at the
center of the beam, while light without OAM has maximum intensity at the
center. Hence a very simple pinhole efficiently distinguishes light with and
without OAM. The intensity distribution of a Laguerre-Gaussian beam, a
paraxial beam possessing OAM $l\hbar$, is given by [1]:

$I_{LG}(l;\rho)=\frac{I_{0}}{\int_{0}^{\infty}duu^{|l|}e^{-u}L_{|l|}(u)}\left(\frac{\sqrt{2}\rho}{w_{0}}\right)^{|l|}L_{|l|}\left(\frac{2\rho}{w_{0}^{2}}\right)e^{-\frac{\rho^{2}}{w_{0}^{2}}}$
(4)

Where $L_{l}(x)$ is the $l$th order Laguerre Polynomial. Thus the transmission
$T(l)$ through a pinhole with a radius $a_{0}$ (normalized by the waist of the
$l_{0}=0$ Gaussian beam) is:

$T(l,a_{0})=\left.{\int_{0}^{a_{0}}\int_{0}^{2\pi}\rho d\rho d\phi
I_{LG}(l;\rho)}\right/{\int_{0}^{\infty}\int_{0}^{2\pi}\rho d\rho d\phi
I_{LG}(l;\rho)}$ (5)

Figure 4(a) shows $T(l,a_{0})$ vs. $a_{0}$ for $l=0-3$. The transmission
decreases sharply with increasing $l$ when $a_{0}$ is smaller than $\sim 0.8$.
Choosing $a_{0}=0.8$, we show in Fig. 4(b) the nearly exponential decrease of
$T(l,a_{0})$ with $l$. These values are also marked by the red vertical lines
in Fig. 3(a). Due to the fast decrease of $T(l,a_{0})$ from $l=0$ to $l\geq
1$, a large extinction ratio is readily achieved, which is very well
approximated by:

$\eta(l_{0})=P(l_{0};l_{0})/\sum P(l\neq
l_{0};l_{0})\approx\frac{\alpha_{out}|p_{V}(0)|^{2}|p_{H}(1)|^{2}}{|p_{V}(1)|^{2}}.$
(6)

$\eta$ is essentially the same for all OAM components, and it is mainly
determined by how well the QZI can distinguish between states with OAM values
$l_{0}=0$ and $l_{0}=1$. We plot in Fig. 5(a) $\eta$ vs. $N$ for
$|\alpha|^{2}=0.9-1$. In general, $\eta$ increases with $N$ but decreases with
$|\alpha|^{2}$, resulting in an optimal $N$ for each $|\alpha|^{2}<1$. Even
for $|\alpha|^{2}=0.9$, $\eta>70$ can be reached with $N=7$. For
$|\alpha|^{2}=0.96$, $\eta$ peaks at $\sim 180$.

Fig. 6: (a)The probability of measuring an OAM value $l$ for a given input
state $l_{0}$ (Equation 2), using pinhole as the OAM filter, $N=8$,
$|\alpha|^{2}=0.96$, and misalignment of 10% and 1%, respectively, at the
pinhole filter and VPP. Despite the decrease in probability for the diagonal
elements at large $l_{0}$, the off diagonal elements decrease much faster, as
implied by the large extinction ratios. (b) The diagonal elements of (a) as a
function of $N$ for $l_{0}=0-10$.

An additional source of error is due to the misalignment of the beam through
two OAM-sensitive components: the OAM filter (e.g. a pinhole) and the VPP.
Misalignment at the pinhole filter leads to reduced coupling efficiency of the
zero OAM state, increased transmission of non-zero OAM states, and thus
reduced extinction ratio. Misalignment on the VPP changes the desired OAM
state into a superposition with neighboring OAM orders. However, these
neighboring orders have very small amplitudes (e.g. $<1\%$ with $1\%$
misalignment) [4], and they are further filtered out through the QZI loop,
resulting in negligible reduction in the extinction ratio. The main effect of
misalignment at VPP is the slightly reduced transmission of the correct OAM
state, hence reduced overall detection probability. We illustrate the effects
of misalignment on the extinction ratio in Fig. 5(a) (the green crosses),
assuming conservatively $10\%$ misalignment of the focused beam waist at the
pinhole and $1\%$ misalignment of the collimated beam waist at the VPP.
Extinction ratios over 100 are still readily achieved.

The extinction ratio can be increased by many orders of magnitude if we only
need to measure every other order, or every third order of OAM (Figure 5(b)).
Correspondingly, we can choose smaller aperture sizes and a VPP that reduces
the $l$ by $\Delta l=2$ or $3$ per passing. A smaller aperture size also
introduces extra loss, but only in the final QZI on the zero OAM state, and
thus only decreases the detection probability by about a factor of two.

To evaluate the overall performance of the OAM spectrometer, we plot in Fig.
6(a) $P(l;l_{0})$ vs. $l$ and $l_{0}$ on the log scale, including loss and
misalignment. The diagonal elements $P(l_{0};l_{0})$ correspond to correctly
detecting an OAM component. They are two orders of magnitude higher than
neighboring off diagonal elements, consistent with the high extinction ratios
calculated before. In Fig. 6(b), we show $P(l_{0};l_{0})$ as a function of $N$
for different $l_{0}$. $N\sim 8$ gives the highest probability for detecting
high order OAM components, while still maintaining an extinction ratio of
above 100.

In summary, we present a compact OAM spectrometer that disperses light of
different OAM values in time. Loss is significant for high order OAM
components with commercially available optical components. However, the high
loss doesn’t have an appreciable effect on the signal to noise ratio;
extinction ratios of $>\\!100$ are readily achieved even after taking into
account optical loss and misalignments. The extinction ratio can be further
improved by many orders of magnitude by skipping OAM orders, or by using a
better OAM filter than a simple pinhole.


11institutetext: Naver Labs Europe
6 chemin Maupertuis, 38240 Meylan, France
<EMAIL_ADDRESS>
# Multi-task learning of time series and its application to the travel demand

Boris Chidlovskii

# Multi-task learning of time series and its application to the travel demand

Boris Chidlovskii

###### Abstract

We address the problem of modeling and prediction of a set of temporal events
in the context of intelligent transportation systems. To leverage the
information shared by different events, we propose a multi-task learning
framework. We develop a support vector regression model for joint learning of
mutually dependent time series. It is the regularization-based multi-task
learning previously developed for the classification case and extended to time
series. We discuss the relatedness of observed time series and first deploy
the dynamic time warping distance measure to identify groups of similar
series. Then we take into account both time and scale warping and propose to
align multiple time series by inferring their common latent representation. We
test the proposed models on the problem of travel demand prediction in Nancy
(France) public transport system and analyze the benefits of multi-task
learning.

##  1 Introduction

Time-dependent entities and their analysis represent a great challenge in
multiple domains, including finance, science, medicine, entertainment, etc. In
the domain of transport, intelligent transportation systems have been designed
to track and analyze a large number of time-dependent events, such as vehicle
position, road load, travel time and demand, service reliability, traffic
density, etc. Efficient planning and management of transportation network, in
particular, an adequate response to changing traffic conditions, require an
accurate modeling and real-time prediction of all these time-dependent
entities.

It is common to represent time-dependent entities as time series. In this
paper we address the problem of simultaneous modeling and prediction for a set
of time series. Our work is motivated by the problem of travel demand
prediction in a public transport system. Figure 1 shows the public transport
network of Nancy city (France), Figure 2 reports three time series
representing the daily passenger travel demand in different network locations.
In these series, one item refers to the number of passengers boarding public
transport vehicles (buses and trams) during a certain period of time. Similar
and dissimilar patterns are easily recognizable in the series, they are often
relevant to peak and off-peak times.

Figure 1: Nancy public transport network.

In statistics and signal processing, autoregressive models have been
intensively used to study time-varying processes. An autoregressive (AR) model
specifies that the output variable depends linearly on its own previous
values. Combined with the integrated moving average (IMA) it forms a more
general ARIMA model of time series; the later can be further generalized to a
SARIMA model that takes into account seasonal effects [11]. Multiple
techniques developed by the machine learning community, such as neural
networks, perceptrons, support-vector machines and regressions, have been also
adapted to model time series [17].

In the transportation domain, the time series has been a subject of multiple
studies in the context of traffic flow forecast [20]. The most recent study
presents an extensive experimental protocol allowing to compare both
autoregressive models and machine learning techniques, and to infer a general
picture of the advantages and limitations of each approach [9].

In the most conventional setting, all time series are modeled individually.
However, if they are related and share important information, a multi-task
learning approach can be proposed to model all of them in a joint way.

Complex dependencies between multiple time series is of high interest for
anomaly detection [14, 4]. Causality analysis discovers the temporal
dependencies between a large number of time series data coming from different
source, like sensors, logs or physical measurements. Inferring the Granger
causality between series became particularly popular, for its scalability and
interpretable insights on the possible reasons of anomalies [14].

The multi-task learning is aimed at leveraging the information of multiple,
mutually related learning tasks to make more accurate predictions for the
individual tasks [3]. Related information contained in tasks can be exploited
to mutually increase the quality of predictions. Examples of successful
application domains for multi-task learning include computational biology,
natural language processing, computer vision [8], where multiple biological,
textual and visual object classes may share some of the relevant features.

In the multi-task learning, the prediction accuracy in each task is leveraged
by making use of data from the other tasks. Previous work has focused on
relating the tasks via regularization [5], mutualization [6], covariance
sharing [2] and mainly addresses the classification case. In this paper we
extend the regularization-based technique [5] to the time series.

Figure 2: Three daily travel demand time series.

Our motivation example for developing a multi-task learning approach is the
travel demand prediction in a public transport network. We integrate a number
of relevant features in prediction models.

We first collect historical data of the travel demand in a public transport
system. It is extracted from a massive collection of individual passenger
boarding events in the network observed over some period of time. Second, we
provide the calendar data, including weekends, school vacation, public
holidays; these data are often critical to estimate the level of travel
demand. Third, we integrate weather information, such as temperature,
precipitation and wind speed, as factors that may also impact the demand.
Lastly, we can benefit from the traffic information, describing the traffic
status in a network segment relevant to a given task.

Unlike multi-task learning for the classification, the mutual relatedness is
harder to assess in time series, as they are often a subject of time shifting,
distortion and re-scaling effects. To address this problem, we first use the
dynamic time warping distance measure, it helps identify similar series and
group them together. Then we develop an alternative approach which attempts to
infer a latent representation common to all observed time series. Such a
common profile is used to align all series in Viterbi-like manner used in
Markovian models.

The contribution of this work can be summarized as follows:

  1. 1.

We develop a multi-task learning approach to the joint learning for a set of
time series.

  2. 2.

For a set or related time series, we design a support vector regression model
derived from regularization-based techniques developed for the classification
case.

  3. 3.

For a given set of time series, we first use the dynamic time warping to
measure their similarity; this helps identify clusters of related time series.

  4. 4.

We then propose another approach to address both time and scale warping of a
set of time series, it infers their common latent profile and use this profile
to align them among themselves.

  5. 5.

We test both methods on the travel demand prediction in Nancy public transport
system and analyze advantages of each method.

The remainder of the paper is organized as follows. The problem of time series
prediction in single and multi-task setting is introduced in Section II. In
Section III, we describe the regularization-based support vector regression
for multi-task learning. Section III.A addresses the series relatedness and
describes how the dynamic time warping distance can be used to identify
similar time series. Section IV proposes an alternative approach to the series
relatedness problem; it introduces the generative model for inferring a common
latent representation of a set of time series and aligns them by both time and
scale warping. Section V reports evaluation results of experiments with multi-
task learning on a travel demand data collected in real public transport
network and analyses the results. Finally Section VI concludes the paper.

##  2 Time Series Prediction

Time series analysis and prediction theory have been intensively studied over
last 40 years [17, 11]. Well known applications are financial market
prediction, environment and weather forecasting, biology, control systems, and
many other applications involving non-linear processes [13].

The goal of time series prediction is to estimate value $y$ at time $i$ based
on past values $y_{i-1}$, $y_{i-2},\ldots$. It is common to limit the depth to
the last $k$ values: $y_{i}=f(y_{i-k},y_{i-k+1},...,y_{i-1}),i=k,\ldots,n,$
where $k\geq 1$ is the modeling depth (lag). Moreover, it can take into
account a number of external characteristics ${\bf v}_{i}$ at time $i$ which
may have impact on the series values:

$y_{i}=f({\bf v}_{i},y_{i-k},y_{i-k+1},...,y_{i-1}),i=k,\ldots,n.$ (1)

The objective of time series prediction is to find a function $f({\bf
x}_{i})$, where ${\bf x}_{i}=({\bf v}_{i},y_{i-k},y_{i-k+1},...,y_{i-1})$ such
that ${\hat{y}}_{i}$, the predicted value of the time series at a future point
in time, is consistent and minimizes a regularized fit function.

The underlying system models and time series data generating processes are
generally complex and the models for these systems are usually not known in
advance. Accurate and unbiased estimation of the time series data produced by
these systems cannot always be achieved using linear techniques, and thus the
estimation process requires more advanced time series prediction algorithms.

###  2.1 SVM for Time Series Prediction

The Support Vector Machine (SVM) is machine learning technique which has been
successfully used for different tasks such as pattern recognition, object
classification, and in the case of time series prediction and regression
analysis. In Support Vector Regression (SVR), a function for a time series is
estimated using observed data. Being a data-driven methodology, SVMs differ
from more traditional time series prediction methodologies, like ARIMA or
Kalman Filter [17].

The underlying motivation for using SVMs for time series is the ability to
make accurate prediction when the underlying system processes are typically
nonlinear, non-stationary and not defined a-priori. In a number of
applications, SVMs have been proven to outperform other non-linear techniques
including neural-network based non-linear prediction techniques such as multi-
layer perceptrons [9, 20, 12].

Given a set of time series data $({\bf x}_{i},y_{i})$, $i=1,\ldots,n$, a time
series prediction method defines a function $f({\bf x})$ that will have an
output equal to the predicted value for some prediction horizon. By using
regression analysis, the prediction function for the linear regression is
defined as $f({\bf x})={\bf w}^{T}{\bf x}+b$. If the data is not linear in its
input space, the data ${\bf x}$ are mapped into a higher dimension space, via
a kernel function $\phi({\bf x})$, then perform a linear regression in the
higher dimensional feature space [12, 16]:

$f({\bf x})={\bf w}^{T}\phi({\bf x})+b.$

The goal is therefore to find optimal weights ${\bf w}$ and threshold $b$, as
well as to define the criteria for finding an optimal set of weights.

Given training data $({\bf x}_{i},y_{i}),i=1,\ldots,n$, ${\bf x}_{i}\in{\cal
R}^{d}$, $y_{i}\in{\cal R}$, SVM first maps input vectors ${\bf x}$ onto the
feature space $\Phi,\phi({\bf x})\in\Phi$ and then approximates the regression
by a linear function $f({\bf x})={\bf w}^{T}\phi(x)+b$. This is obtained by
solving the following optimization problem in the $\varepsilon$-insensitive
tube [19]:

$\begin{array}[]{ll}{\rm min}&\frac{1}{2}||{\bf
w}||^{2}+C\sum_{i=1}^{n}(\xi_{i}+\xi^{*}_{i})\\\ {\rm
s.t.}&\left\\{\begin{array}[]{l}y_{i}-({\bf w}^{T}\phi({\bf
x}_{i})-b)\leq\varepsilon+\xi_{i}\\\ ({\bf w}^{T}\phi({\bf
x}_{i})+b)-y_{i}\leq\varepsilon+\xi_{i}^{*},i=1,\ldots,n\\\ \xi_{i}^{*}\geq
0,\xi_{i}\geq 0,i=1,\ldots,n\\\ \end{array}\right.\end{array}$ (2)

where $\xi_{i},\xi_{i}^{*},i=1,\ldots,n$ are slack variables, measuring the
deviation from $\varepsilon$-insensitive tube, $C$ is the regularization
parameter.

##  3 Multi-task SVM Regression

Suppose we dispose a set of $m$ related time series, $m>1$. Training samples
are represented as a set $\\{({\bf
x}_{ri},y_{ri}),r=1,\ldots,m,i=1,\ldots,n\\}$. We denote the indices from
series $r$ by $T_{r}=\\{i_{r1},\ldots,i_{rn}\\},r=1,\ldots,m.$ 111Without loss
of generality, we make assumption of the same length $n$ for all series.

Making an assumption of relatedness between the tasks serves us to design the
multi–task learning method. We assume that the tasks are related in a way that
the true models are all close to some common model ${\bf w}$ [5] and every
task model ${\bf w}_{r}$ can be written as ${\bf w}+{\bf w}_{r}$,
$r=1,\ldots,m$, where the individual vectors ${\bf w}_{r}$ are small when the
tasks are similar to each other.

We adapt the regularization approach [5], according to which vectors of each
task ${\bf x}_{i},i\in T_{r}$ are mapped into two different spaces. One is a
space $\Phi$ common to all tasks, $\phi({\bf x}_{i})\in\Phi$; another is a
correction space $\Phi_{r}$, specific to every task $r$, $\phi_{r}({\bf
x}_{i})\in\Phi_{r}$.

Individual task information is accounted in the slack variables, defined as
follows:

$\xi_{i}$ | = | ${\bf w}_{r}^{T}\phi({\bf x}_{i})+b_{r},i\in T_{r},r=1,\ldots,m$
---|---|---
$\xi_{i}^{*}$ | = | ${{\bf w}_{r}}^{T}\phi({\bf x}_{i})+b^{*}_{r},i\in T_{r},r=1,\ldots,m.$
(3)

All slack variables are non-negative, $\xi_{r}({\bf x}_{i}),\xi^{*}_{r}({\bf
x}_{i})\geq 0,r=1,\ldots,m$. Thus samples maped in the correcting space have
to lie on one side of the corresponding correcting function. The correcting
function also has to pass through some points with slack variables being zero.

The goal of the multi-task version of SVR is to estimate $m$ regression
models, one model per task. Multi-task learning SVR incorporates individual
task information into estimated regression functions through the slack
variables. We specify the following parameterized functions for $m$ regression
models:

$f({\bf x})+f_{r}({\bf x})={\bf w}^{T}\phi({\bf x})+b+{\bf
w}_{r}^{T}\phi_{r}({\bf x})+b_{r},r=1,\ldots,m,$ (4)

where $f({\bf x})={\bf w}^{T}\phi({\bf x})+b$ is the common estimation
function and $f_{r}({\bf x})={\bf w}_{r}^{T}\phi_{r}({\bf x})+b_{r}$ is a
correction function for series $r$, $r=1,\ldots,m$. Therefore the MTL SVR
formulation solves the following optimization problem:

$\begin{array}[]{ll}{\rm min}&\frac{1}{2}||{\bf
w}||^{2}+\mu\sum_{i=1}^{t}||{\bf
w}_{r}||^{2}+C\sum_{i=1}^{n}(\xi^{r}_{i}+\xi^{r*}_{i})\\\ {\rm
s.t.}&\left\\{\begin{array}[]{l}y^{r}_{i}-({\bf w}^{T}\phi({\bf
x}_{i})+b+{{\bf w}_{r}}^{T}\phi_{r}({\bf
x}_{i})+b_{r})\leq\varepsilon+\xi^{r}_{i}\\\ ({\bf w}^{T}\phi({\bf
x}_{i})+b+{{\bf w}_{r}}^{T}\phi_{r}({\bf
x}_{i})+b_{r})-y^{r}_{i}\leq\varepsilon+\xi^{r*}_{i}\\\ \xi^{r*}_{i}\geq
0,\xi^{r}_{i}\geq 0,i=1,\ldots,n,r=1,\ldots,m.\\\
\end{array}\right.\end{array}$ (5)

where $\mu$ is a regularization parameter for the individual correction
functions. Using the dual optimization technique (similar to standard SVM),
the dual form of the above optimization problem is as follows:

$\begin{array}[]{ll}{\rm
max}&-\varepsilon\sum_{i=1}^{n}({\alpha}_{i}^{*}+{\alpha}_{i})+\sum_{i=1}^{n}({\alpha}_{i}^{*}-{\alpha}_{i})y_{i}\\\
&-\frac{1}{2}\sum_{i,j=1}^{n}({\alpha}_{i}^{*}-{\alpha}_{i})({\alpha}_{j}^{*}-{\alpha}_{j})\phi({\bf
x}_{i})^{T}\phi({\bf x}_{j})\\\ &-\frac{1}{2\mu}\sum_{r=1}^{t}\sum_{i,j\in
T_{r}}({\alpha}_{i}^{*}-{\alpha}_{i})({\alpha}_{j}^{*}-{\alpha}_{j}){\phi_{r}({\bf
x}_{i})}^{T}\phi_{r}({\bf x}_{j})\\\ {\rm
s.t.}&\left\\{\begin{array}[]{ll}\sum_{i\in
T_{r}}({\alpha}_{i}^{*}-{\alpha}_{i})=0,r=1,\ldots,m,\\\
0\leq{\alpha}_{i},{\alpha}_{i}^{*}\leq C,i=1,\ldots,n.\\\
\end{array}\right.\end{array}$ (6)

Then vectors ${\bf w}$, ${\bf w}_{r}$ can be expressed in terms of training
samples:

$\begin{array}[]{rll}{\bf
w}&=&\sum_{i=1}^{n}({\alpha}_{i}^{*}-{\alpha}_{i})\phi({\bf x}_{i}),\\\ {\bf
w}_{r}&=&\frac{1}{\mu}\sum_{i\in
T_{r}}({\alpha}_{i}^{*}-{\alpha}_{i})\phi_{r}({\bf x}_{i}),r=1,\ldots,m.\\\
\end{array}$ (7)

Note that in addition to usual parameters of individual SVR, $C$ and
$\varepsilon$, the MTL extension requires also to tune parameter $\mu$.

###  3.1 Multiple Series Relatedness

In the previous section, multi-task learning was assumed to be applied to a
set of mutually related time series. Ensuring the task relatedness is crucial;
applying the multi-task learning to unrelated series may have a negative
impact on the prediction quality.

A pair of time series is said to be mutually related if they behave similarly
and show up the same tendencies. To address the relatedness between time
series, we first consider conventional measures of similarity which is of
fundamental importance for time series analysis [7]. The straightforward
approach to measure the similarity between two time series is to evaluate the
Euclidean distance on their normalized representation. However, Euclidean
distance ignores any shifting and distortion in time series and has been found
unsuitable in different domains [18].

Besides the Euclidean-based distance measure, the time-warping distance has
become one of the most popular and field-tested similarity measures. The
dynamic time warping (DTW) distance is efficient as the time-series similarity
measure for minimizing the effects of shifting and distortion in time; it
allows elastic transformation of time series in order to detect similar shapes
with different phases. Given two time series, DTW algorithm yields the optimal
solution in the $O(n^{2})$ time.

Since many time series mining algorithms use similarity search as a core
subroutine, To avoid that the time taken for DTW similarity becames a
bottleneck, multiple efforts have been undertaken to improve the complexity of
the DTW distance algorithms. In recent work [15], a novel DTW similarity
algorithm for time series is proposed in order to search and mine massive time
series. It permits to search with DTW similarity faster than with the
Euclidean distance algorithm. For such problems as motif discovery and time
series clustering, it scales up to very large datasets accounting for billions
of time series.

Below we use the DTW similarity to identify groups of similar time series and
then apply the multi-task learning to each of the identified clusters.
Algorithm 1 below summarizes the multi-task learning on a given set of time
series.

0: Set of time series ${\bf y}_{1},\ldots,{\bf y}_{m}$, ${\bf
y}_{r}=(y_{r1},\ldots,y_{rn})$, $r=1,\ldots,m$

0: Number of clusters $L$

0: Observations $({\bf v}_{r},\ldots,{\bf v}_{rn})$, including weather,
calendar and traffic data, for every series ${\bf y}_{r}$.

1: for every pair $({\bf y}_{r},{\bf y}_{t}),1\leq r<t\leq m$ do

2: Measure the similarity between series ${\bf y}_{r}$ and ${\bf y}_{t}$ by
the dynamic time warping (DTW) distance.

3: If series are too long (more 10,000 elements), uniformly sample series
before applying the DTW algorithm.

4: end for

5: Apply agglomerative clustering algorithm to form $L$ clusters
$G_{1},\ldots,G_{L}$ of series using DTW distance as metric. Any series ${\bf
y}_{r}$ belongs to exactly one group.

6: for each group $G$ of relates tasks do

7: Apply the multi-task learning described in Section 3

8: end for

9: return Multi-task models for all groups and tasks

Algorithm 1 Multi-task learning for related and unrelated time series.

##  4 Time series common profile

In the previous section, we grouped time series by using their DTW distance as
similarity metric. This helps identify similar tasks but does not guarantee
their relatedness.

In this section, we make a step further in the analysis of time series
relatedness, by taking into account both time warping and scale warping when
aligning two or more time series. Our approach is based on an assumption that
all series in the set are noisy and rescaled replicates of one common and
latent representation [10]. To leverage the information contained in such
noisy replicates, such a common representation can be inferred and used to
align the series among themselves in an appropriate way.

The latent trace is modeled as an underlying, noiseless representation of the
set of observed noisy time series. Output time series are generated from this
model by moving through a sequence of hidden states and emitting an observable
value at each step, as in Hidden Markov Machines (HMMs). To account for
changes in the value amplitude within and across replicates, the latent time
states are augmented by a set of scale states, which control how the emission
signal will be scaled relative to the value of the latent trace. The latent
trace is learned in unsupervised manner, as well as the transition
probabilities controlling the Markovian evolution of the scale and time states
and the overall noise level observed data.

The Continuous Profile Model (CPM) [10] is a generative model for a set of $m$
time series, ${\bf y}_{1},{\bf y}_{2},...,{\bf y}_{m}$. Time series in the set
are noisy replicas of a common and latent representation, defined as a trace
${\bf z}=(z_{1},\ldots,z_{M})$. Any given time series ${\bf y}_{r}$ is a sub-
sampled version of ${\bf z}$ after applying local scale transformations. Size
$M$ of ${\bf z}$ is larger than time series size $n$ so that any experimental
data could be mapped precisely to the correct underlying trace point. In
practice, $M=\lambda n$ with $1<\lambda<3$.

The resolution of the latent trace ${\bf z}$ is higher with respect to the
observed time series, so the experimental time can be made effectively to
speed up or slow down by advancing along the latent trace in larger or smaller
jumps. The sub-sampling and local scaling used during the generation of each
observed time series are determined by a sequence of hidden state variables.
The state sequence for series ${\bf y}_{r}$ is denoted ${\pi^{r}}$. Each state
in the sequence maps to a pair (time state, scale state),
$\pi^{r}_{i}\rightarrow\\{\tau_{i}^{r},\psi_{i}^{r}\\}$, where time state
$\tau_{i}^{r}\in(\tau_{1},\ldots,\tau_{M}$), and scale states
$\psi_{i}^{r}\in(\psi_{1},\ldots,\psi_{Q}$), $Q$ is the number of scale
states.

The probability of emitting observed value $y_{i}^{r}$ in state $\pi_{i}^{r}$
is defined by the emission probability distribution
$A_{\pi_{i}^{r}}(y^{r}_{i}|{\bf z})\equiv p(y_{i}^{r}|\pi_{i}^{r},{\bf
z},\sigma,u^{r})$ which is defined by probability distribution ${\cal
N}(y_{i}^{r};z_{\tau_{i}^{r}}\psi_{i}^{r}u^{r},\sigma)$, where $\sigma$ is the
noise level of the observed data, ${\cal N}(a;b,c)$ denotes a Gaussian
probability density for $a$ with mean $b$ and standard deviation $c$. The
$u_{r}$ are real-valued scale parameters, one per observed time series, that
correct for any overall scale difference between time series $r,r=1,\ldots,m$,
and the latent trace ${\bf z}$.

The transitions between time states are independent from transitions between
scale states, so
$p(\pi_{i}^{r}|\pi^{r}_{i-1})=p(\psi_{i}|\psi_{i-1})p^{r}(\tau_{i}|\tau_{i-1})$.
To enforce that time must move forward, any time state can jump ahead no more
than $J_{\tau}$ states. Similarly, scale state transitions are allowed between
neighboring scale states only, this disallows arbitrary jumps on the local
scale. These constraints keep the number of legal transitions to a tractable
computational size. Each observed time series ${\bf y}^{r}$ has its own time
transition probability distribution to account for experiment-specific
patterns. The time transition probability distribution for ${\bf y}^{r}$ is
given by multinomial distribution as follows

$p^{r}(\tau_{i}=a|\tau_{i-1}=b)=\left\\{\begin{array}[]{ll}d_{l}^{r},&{\rm if\
}a-b=l,l=1,\ldots,J_{\tau},\\\ 0,&{\rm otherwise.}\\\ \end{array}\right.$ (8)

where $\sum_{i=1}^{J_{\tau}}d^{r}_{i}=1$. The time scale probability
transition for all series is given by

$p(\psi_{i}=a|\psi_{i-1}=b)=\left\\{\begin{array}[]{ll}s_{0},&{\rm if\
}D(a,b)=0,\\\ s_{1},&{\rm if\ }D(a,b)=\pm 1,\\\ 0,&{\rm otherwise,}\\\
\end{array}\right.$ (9)

where $D(a,b)=\pm 1$ means that $a$ is one scale state larger (+1) or smaller
(-1) than $b$, and $D(a,b)=0$ means that $a=b$ and the scale state does not
change. The distribution $J_{\tau}$ are constrained by $2s_{1}+s_{0}=1$.

$J_{\tau}$ determines the maximum allowable speedup of one portion of a time
series relative to another portion, within the same series or across different
series. However, the length of time for which any series can move so rapidly
is constrained by the length of the latent trace; thus the maximum overall
ratio in speeds achievable by the model between any two entire time series is
given by $min(J_{\tau},\lambda)$.

After training, one may examine either the latent trace or the alignment of
each observed time series to the latent trace. Such alignments can be achieved
by several methods, including use of the Viterbi algorithm to find the highest
likelihood path through the hidden states [1], or sampling from the posterior
over hidden state sequences. Viterbi alignment is found to work well in our
experiments.

Training the CPM with the EM algorithm is similar to HHMs with Baum-Welch
algorithm. The $E$-Step is computed exactly using the Forward-Backward
algorithm, which provides the posterior probability over states for each time
point of every observed time series and the pairwise state posteriors. The
$M$-Step consists of a series of analytical updates to the various parameters,
see [10] for more detail.

Figure 3 shows an example of alignment for five travel demand times series.

Figure 3: Time and scale alignment of 5 travel demand time series.

##  5 Evaluation

We run experiments to test learning methods described in previous sections on
the task of travel demand prediction in Nancy (France) public transport
network. Two alternative approaches to model a set time series are tested in
different settings and scenarios. One is a single-task learning approach where
a separate support vector regression model is estimated for each task
independently. Another is the multi-task learning approach where one common
and a number of individual correction models are jointly learned to estimate
related regression models. A set of time series representing the travel demand
are analyzed and pre-processed. In one instantiation, they are clustered by
using the DTW distance as described in Algorithm 1. In another, a latent
profile common to all the series is inferred and then all series are aligned
to this common profile, as described in Section IV.

Public is GDP time series include annual per capita GDP time series for
several countries. One interest lies in studying the ”periodic” behaviour of
such series in connection with understanding business cycles. Another lies in
forecasting turning points. Another lies in comparisons of characteristics of
the series across national economies.

We dispose the following dataset to evaluate different prediction models:

  * •

Historical data of the travel demand, which accounts for the total of 5,29 M
passenger boarding events registered at 1,067 stops and 28 tram and bus lines
in the city network; the dataset covers the period of 90 days from March to
May 2012, including 13 weekends, 2 weeks of school vacations and three public
holidays. Each boarding transaction includes the route, stop, direction and
the timestamp of the event. Time series describing the travel demand are then
generated by the aggregation of boarding timestamps at a regular time period.
In the tests, two different periods are used to produce the time series, of 10
and 60 minutes.

  * •

Weather data include the air temperature (minimum and maximum), air humidity,
precipitation, wind speed, etc. All data are registered hourly 222Available
from weather station Nancy-Ochey at http://www.meteociel.fr/temps-reel/
obs_villes.php?code2=7181..

  * •

Calendar data, including the binary tags for weekends, school vacations and
public holidays.

  * •

Traffic information on the network segment relevant to a given task. The
average bus adherence delay accumulated between timestamps $i-1$ and $i$ for
each series is used.

####  5.0.1 Protocol and experimental setup

We compare the performance of different methods for time series by using the
standard measure of Mean Square Error (MSE). It is defined as follows

$E_{MSE}=\frac{1}{n}\sum\limits_{i=1}^{n}\left(y_{i}-\hat{y}_{i}\right)^{2},$
(10)

where $y_{i}$ is the true value at time point $i$ in the test set while
$\hat{y}_{i}$ is the prediction.

Experiments for the travel demand prediction are conducted on two levels of
the public transport network; one on city zones, another for individual stops.
For the zone demand experiment, a preliminary analysis of the travel dataset
has been done, with the goal to extract all travel patterns required for the
inference of origin-destination (OD) matrices. The city network is then
segmented into zones, where one zone is defined as a compact city region with
stops showing up a high similarity of their travel patterns. Figure 5
demonstrates the segmentation of Nancy city network into 15 zones, where blue
points refer to the bus and tram stops and each zone is shown as a convex hull
built around geographical locations of all stops included in the zone. For
this case, we generate 15 times series describing the travel demand in all
zones.

For the stop demand experiment, one tram and four bus lines have been
selected. For each line, a representation set of 10 stops has been selected to
produce the set of travel demand time series.

The evaluation prototype has been implemented in MATLAB as extension of MALSAR
toolkit for multi-task learning via structured regularization [21].

Figure 4: GDP per capita: a) 9 series, b) 1T and MT learing. Figure 5: Nancy
city zoning by travel demand.

In all experiments, all time series are split into training and testing sets
in the same way. The training set includes starting $P$% items (days); the
testing set includes remaining $100-P$%, where $P$ varies from 40 to 90.

For each of two experiments, we evaluate three different learning approaches
in different scenarios:

  1. 1.

We compare the prediction performance of the three following methods:

     * •

Single-task (ST) learning, where individual SVR models are learned and tested
independently; the average of their MSE values is reported;

     * •

Multi-task (MT) learning, where one common SVR model and one correction model
per time series are learned; the series are clustered by their DTW distance,
as described in Section 3.

     * •

MT learning with CPM alignment, where the time series get aligned to their
common profile, as described in Section 4.

  2. 2.

We run experiments with two different demand aggregation periods, of 10
minutes and 1 hour.

All time series are normalized and scaled to the logarithmic space. In SVR,
the linear kernel function $\phi({\bf x})={\bf x}$ is preferred to the
Gaussian one. The depth $k=2$ for the past values (lag) was determined by
cross-validation and remains the same in all experiments. Parameters $C$ and
$\varepsilon$ are automatically tuned as in the standard SVM regression [19],
the optimal value of parameter $\mu$ is determined by re-sampling.

Algorithm 1 for clustering time series by their DTW distance takes as input
the number of clusters $L$. In the zone experiment, where we produce 15 time
series, the number of clusters to test is between 2 and 15. In the stop
experiment, 10 time series represent any selected line, and the number of
clusters is between 2 and 10. The optimal number of clusters is detected by
the cross-validation, the best results are reported for $L$ between 3 and 6 in
the stop experiment and 4 or 5 in zone experiment.

For the CPM inference and alignment, seven latent traces ${\bf
z}_{1},\ldots,{\bf z}_{7}$ are inferred from the daily times series in the
training set, one for each day of week. At the test time, a daily time series
${\bf y}_{i}$ as well as relevant observations ${\bf v}_{i}$ get aligned to
the same day trace ${\bf z}$ using the Viterbi algorithm. The size $M$ of the
latent for series of size $n$ is $\lambda n$ where $\lambda=2.5$. The scaling
space size $Q$ is limited to 5 states, they are evenly spaced on logarithmic
scale. The MSE values are then averaged over 7 day-specific experiments.
Finally, the number of iterations for EM algorithm was limited to 20.

In all figures presented below, the whisker plots are used to show the mean,
variance and outliers of MSE values.

####  5.0.2 Stop experiments

Figure 6 reports evaluation results of single- and multi-task learning for the
stop experiment with 1 hour period. In addition to the three referred methods,
it also includes MSE values for the single-task learning without external
observations ${\bf v}_{i}$ (weather, calendar and traffic features). As one
can see, taking external observations into account reduces the error values by
40%-60%; this measures the impact of external observations on the predictive
behavior of support vector regression models 333Analysis of individual
contribution of external observations is beyond the scope of this study.. The
multi-task learning with DTW outperforms the single-task learning by 2%-7%,
multi-task learning with CPM further reduces the error by 3%-12%.

Figure 7 reports the error values for the stop experiment, this time with 10
minute period. Shorter period of data aggregation leads to a higher variance
in time series; this results in higher error values. Both multi-task learning
methods improve the predictive performance with respect to single-task
learning, with multi-task learning with CPM being the best.

Figure 6: Stop demand prediction for 1 hour period.  Figure 7: Stop demand
prediction for 10 minute period.

####  5.0.3 Zone experiments

Figure 6 reports the evaluation results of the referred methods for the zone
experiment in 0-delay scenario, with 1 hour period. The multi-task learning
with DTW outperforms the single-task learning by 0%-9%, multi-task learning
with CPM reduces the error by 3%-12%.

Figure 6 reports results for the same experiment in 1-delay scenario.
Unavailability of the last value in series makes the prediction harder for all
methods, both in terms of the mean and the variance. The multi-task with DTW
plays the equal game with the single-task learning or outperforms it in some
cases.

In this experiment, multi-task learning with DTW performs better than the
single-task learning for few values of $L$; in other cases, it often shows a
worse performance (see Figure 9.a for cases $P$=40% and 90%). The multi-task
learning with CPM resists better to the data delay, it outperforms the single-
task learning by 3%-16%.

Figure 8: Zone demand prediction.

Figure 9: a) MT with DTW for different number of clusters. b) MT with CPM for
different numbers of stops.

####  5.0.4 Discussion

Analysis of the evaluation results in the previous subsections suggests that
multi-task learning can offer a comfortable performance improvement over
single-task setting. However, unlike the multi-task learning for
classification, the task relatedness raises a particular concern. A simple
normalization is insufficient and advanced techniques are required to leverage
the information shared by different time series.

Between two methods presented in Sections III and IV, the continuous profile
model (CPM) gives an advantage over the DTW distance based clustering.
Considered as a straightforward choice to identify related time series,
Algorithm 1 is sensible to the way the times series get clustered. Moreover,
DTW-based alignment does not guarantee the relatedness of selected tasks.

Time series alignment with CPM seems to solve the problem, the method appears
more stable, it yields a performance improvement in a majority of experiments
and scenarios. The number of time series we tested in both zone and stop
experiments is modest and do not raise particular issues when inferring daily
latent traces, even for 10 minute period. In Figure 9.b we relax the number of
stops in the stop experiment. Instead of 10 representative stops, we pick up
randomly 2 to 20 stops, and compare multi-task CPM to single-task learning,
for $P$=40% and 90%.

Computationally, multi-task DTW and CPM methods behave similarly, due to the
straightforward dynamic programming implementation of the DTW algorithm and a
reasonable convergence of the EM algorithm when inferring 7 daily latent
traces and Viterbi-based alignment of the observed time series.

##  6 Conclusion

In a public transport network, an accurate estimation of the number of
travelers entering the network at any period of time is crucial for efficient
planning of transport services. The estimation model is learned from
historical data and relevant observations. For a set of time series describing
the demand, the single-task learning can be further turned into by multi-task
learning in order to leverage the information shared by different time series.
In this paper, we describe the regularization-based support vector regression
for multi-task learning and propose two ways of addressing the series
relatedness. One way is to cluster times series by their dynamic time warping
distance. Another way is to infer a latent representation common to all time
series and then align the series by this common representation. We apply the
proposed methods on the travel demand prediction in Nancy public transport
network and report the evaluation results for single- and multi-task learning
methods.

## References

  * [1] L. E. Baum, T. Petrie, G. Soules, and Weiss N.  A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains.  Ann. Math. Statist., 41(1):164–171, 1970.
  * [2] Edwin Bonilla, Kian Ming Chai, and Christopher Williams.  Multi-task gaussian process prediction.  In Proc. Neural Information Processing Systems (NIPS) 2007, Vancouver, BC, Canada, 2007.
  * [3] R. Caruana.  Multitask Learning.  Machine learning, 28(1):41–75, 1997.
  * [4] S. Chawla, Yu Zheng, and Jiafeng Hu.  Inferring the root cause in road traffic anomalies.  In , 2012 IEEE 12th International Conference on Data Mining (ICDM), pages 141–150, 2012.
  * [5] Theodoros Evgeniou and Massimiliano Pontil.  Regularized Multi-Task Learning.  In Proc. 10th ACM SIGKDD Intern. Conference on Knowledge Discovery and Data Mining, KDD ’04, pages 109–117, New York, NY, USA, 2004. ACM.
  * [6] Jean Baptiste Faddoul, Boris Chidlovskii, Rémi Gilleron, and Fabien Torre.  Learning Multiple Tasks with Boosted Decision Trees.  In Proc. ECML/PKDD, pages 681–696, 2012.
  * [7] Tak-Chung Fu.  A review on time series data mining.  Engineering Applications of Artificial Intelligence, 24(1):164 – 181, 2011.
  * [8] J.M. Leiva-Murillo, L. Gomez-Chova, and G. Camps-Valls.  Multi-task Remote Sensing Data Classification.  IEEE Transactions on Geoscience and Remote Sensing, 51(1):151 –161, January 2013.
  * [9] Marco Lippi, Matteo Bertini, and Paolo Frasconi.  Short-term traffic flow forecasting: An experimental comparison of time-series analysis and supervised learning.  IEEE Transactions on Intelligent Transportation Systems, 14(2):871–882, 2013.
  * [10] Jennifer Listgarten, Radford M. Neal, Sam T. Roweis, and Andrew Emili.  Multiple alignment of continuous time series.  In Proc. NIPS, 2004.
  * [11] Terence C. Mills.  Time Series Techniques for Economists.  Cambridge University Press, 1990.
  * [12] Klaus-Robert Müller, Alex J. Smola, Gunnar Rätsch, Bernhard Schölkopf, Jens Kohlmorgen, and Vladimir Vapnik.  Predicting Time Series with Support Vector Machines.  In Proc. 7th Intern. Conference on Artificial Neural Networks, ICANN ’97, pages 999–1004, London, UK, UK, 1997. Springer-Verlag.
  * [13] S. J. Orfanidis.  Optimum Signal Processing: An Introduction, 2nd ed. New York: Mac-Millan, New York, 1988.
  * [14] Huida Qiu, Yan Liu, N.A. Subrahmanya, and Weichang Li.  Granger causality for time-series anomaly detection.  In Data Mining (ICDM), 2012 IEEE 12th International Conference on, pages 1074–1079, 2012.
  * [15] Thanawin Rakthanmanon, Bilson Campana, Abdullah Mueen, Gustavo Batista, Brandon Westover, Qiang Zhu, Jesin Zakaria, and Eamonn Keogh.  Searching and mining trillions of time series subsequences under dynamic time warping.  In Proc. 18th ACM SIGKDD Intern. Conference on Knowledge Discovery and Data Mining, KDD ’12, pages 262–270, New York, NY, USA, 2012.
  * [16] Stefan Rüping.  SVM Kernels for Time Series Analysis.  In Proc. of Tagungsband der GI Workshop-Woche, 2001.
  * [17] N. Sapankevych and R. Sankar.  Time Series Prediction Using Support Vector Machines: A Survey.  Computational Intelligence Magazine, IEEE, 4(2):24 –38, may 2009\.
  * [18] Pavel Senin.  Dynamic time warping algorithm review.  In Technical Report CSDL-08-04. Information and Computer Science Department. University of Hawaii at Manoa. Honolulu, USA, 2008.
  * [19] Alex J. Smola and Bernhard Schölkopf.  A Tutorial on Support Vector Regression.  Statistics and Computing, 14(3):199–222, 2004.
  * [20] Man-Chun Tan, S. C. Wong, Jian-Min Xu, Zhan-Rong Guan, and Peng Zhang.  An aggregation approach to short-term traffic flow prediction.  IEEE Transactions on Intelligent Transportation Systems, 10(1):60–69, 2009.
  * [21] J. Zhou, J. Chen, and J. Ye.  MALSAR: Multi-tAsk Learning via StructurAl Regularization.  Arizona State University, 2011.