# Field configurable skyrmion devices utilising fine planar motion control

Charles Kind  Department of Computer Science, University of Bristol, Bristol
BS8 1UB, UK

(September 3, 2024)

###### Abstract

Magnetic skyrmions are nano-scale spin textures whose stability, size and ease
of manipulation make them strong contenders for next generation data and logic
applications. Here we illustrate fine motion control of skyrmions, we show how
they can be moved in any direction and for any distance, in the plane. We
demonstrate, using simulations, two novel skyrmionic devices whose function is
predicated upon the fine planar motion control of skyrmions. The first device
is a field configurable logic gate that could be used as a component in a
larger skyrmionic processor. The second is a field configurable skyrmion
spiking neuron device that could be a component in larger skyrmionic neural
networks. Such devices are highly energy efficient and are potentially useful
in the field of neuromorphic computing, which seeks to build ’brainlike’
computers. Finally we explore the concept that this free planar motion, which
enables us to ignore the skyrmion Hall effect, offers an important new
paradigm in skyrmionic device theory moving away from the current applications
such as racetrack memory.

##  I Introduction

Magnetic skyrmions are topologically stable, particle-like, spin textures in
the magnetization of chiral magnets [1]. They have been proposed as a possible
route to a new generation of low power, high density memory [2, 3] and as
candidates for use in neuromorphic computing systems [4, 5, 6].

Many different types of skyrmionic device have been considered since the idea
of using them as media for data and logic applications was first introduced.
The predominant type of device considered is based upon the concept of
magnetic racetrack first introduced by Stuart Parkin et al. in 2008 [7] and
put a skyrmionic perspective by Albert Fert et al. in 2013 [2]. The basic idea
is that we push skyrmions down a thin wire exploiting the low electrical
currents required to drive them to maximise their velocities. Scientists have
considered ferrimagnetic [8] and antiferromagnetic [9] domains to overcome
potentially undesirable physical features of skyrmions such as the topological
Hall effect [10] which tends to push the skyrmion at an angle from the
direction of electron flow. More recently skyrmions have been considered in
the context of logic-in-memory architecture [11] and as artificial neurons
[4], though still using one dimensional (1D) current flow to drive motion. In
this work we consider full two dimensional (2D) skyrmion motion in the plane.
We demonstrate using analysis of the Thiele [12] equation and using
simulations that fine 2D motion control of skyrmions is possible and can
effectively ignore the Hall effect even in ferromagnetic systems. Furthermore,
we model two simple devices that exploit this mechanism demonstrating the
impact this shift in perspective will have on emerging spintronics
applications.

Magnetic skyrmions can be found in the magnetization of certain chiral magnets
[1]. Present in both bulk magnetic materials and interfacially in magnetic
multilayers, they are stabilised by the Dzyaloshinkii–Moriya interaction (DMI)
[13, 14]. The magnetisation vectors of a single skyrmion, in the two
dimensional (2D) continuous model, are a cover of the two-sphere and carry a
unitary degree, $|Q|=1$, defined as

$\displaystyle
Q=\frac{1}{4\pi}\int\mbox{\boldmath$n$}\cdot\left(\partial_{1}\mbox{\boldmath$n$}\times\partial_{2}\mbox{\boldmath$n$}\right)\>d^{2}x,$
(1)

where $\mbox{\boldmath$n$}\left(\mbox{\boldmath$x$}\right)$ is the unit vector
field of magnetization. Under this definition, a single stable skyrmion has
the degree $Q=-1$. This quantized degree is the reason for the topological
stability of the skyrmion. In the continuous mathematical model the energy
cost for moving away from this state is infinite but in real magnetic systems
the spins are defined on a discrete lattice at the atomic level and samples
are finite in extent meaning that the energy barrier is also finite.
Nonetheless the energy barrier is real and high enough that skyrmions in the
systems we consider can be thought of as stable with respect to perturbing
energies below that barrier.

Free planar motion is achieved through the pulsing of current at high
frequency in directions orthogonal to each other. This then forms a basis and
the almost linear response of skyrmion velocity to applied current is
exploited to derive the required vectors. The setup of the two devices is
shown in Fig. 1. The first device, Fig. 1 (a), is a skyrmion configurable
logic gate (SCLG) and the second, Fig. 1 (c), a skyrmion configurable neural
network (SCNN) node based on the work of Xing Chen et al. [15]. Detailed
explanations of their form and function will follow after we explore the
theory of 2D motion.

Figure 1: A skyrmion configurable logic gate and configurable neural network
node. (a) A sketch of a skyrmion configurable logic gate (SCLG). The dashed
yellow areas represent magnetic tunnel junctions where changes in topological
degree cause quantized changes in resistivity thus measuring skyrmion
location. Note the two pairs of electrical contacts, each can pair can swap
the source and sink. (b) A top view of the junction from (a) with planar
dimensions included. (c) A sketch of a skyrmion configurable neural network
(SCNN) node. The dashed yellow areas represent magnetic tunnel junctions.
There are two pairs of electrical contacts as in (a). Here we also see four
red arrows representing separate localised external magnetic fields
($B_{ext}$) applied to each arm of the device. (d) A top view of the junction
from (c) with planar dimensions included. (e) A schematic of a network of
SCNNs. The node Sigma $\left(\sum\right)$ accumulates voltage inputs from
presynaptic nodes and fires at a threshold. The external field and local
frequency of the node can be altered to change its firing characteristics. The
node ’fires’ when the skyrmion enters the end bulb.

##  II Skyrmion motion theory

We apply the Thiele [12] approach to the dynamics of our skyrmionic
magnetisation textures assuming the stationary limit, where the magnetisation
texture moves with constant velocity, and assuming that the texture does not
deform. We consider the current in plane (CIP) geometry using the Zhang-Li
model [16]. The travelling wave ansatz,
$\mbox{\boldmath$n$}\left(\mbox{\boldmath$r$},t\right)=\mbox{\boldmath$n$}\left(\mbox{\boldmath$r$}-\mbox{\boldmath$v$}_{d}t\right)$
, is then applied to the Landau-Lifshitz-Gilbert [17] dynamical equation with
torque to produce [18, 19, 20] Thiele’s equation [12].

$\displaystyle\mbox{\boldmath$G$}\times\left(\mbox{\boldmath$j$}-\mbox{\boldmath$v$}\right)+\mbox{\boldmath$D$}\left(\beta\mbox{\boldmath$j$}-\alpha\mbox{\boldmath$v$}\right)=0,$
(2)

here $\mbox{\boldmath$G$}=4\pi Q\mbox{\boldmath$e$}_{z}$ is the gyrocoupling
vector, $D$ is the dissipative tensor with components

$\displaystyle
D_{ij}=\int\partial_{i}\mbox{\boldmath$n$}\cdot\partial_{j}\mbox{\boldmath$n$}\>d^{2}x,$
(3)

$v$ is the velocity of the centre of mass of the skyrmion, $j$ is the current
vector, $\alpha$ is the Gilbert damping constant and $\beta$ is the
dissipative spin transfer torque parameter. For an axially symmetric texture
the dissipative tensor $D$ will have components $D_{1,1}=D_{2,2}$ and
$D_{i,j}=0$, for $i\neq j$ [21].

Solving Thieles equation, given
$\mbox{\boldmath$v$}=\left(\mbox{\boldmath$v$}_{x},\mbox{\boldmath$v$}_{y}\right)$,
we get

$\displaystyle\mbox{\boldmath$v$}_{x}=\frac{\left(\alpha\beta
D^{2}+G^{2}\right)\mbox{\boldmath$j$}_{x}-GD\left(\alpha-\beta\right)\mbox{\boldmath$j$}_{y}}{\left(\alpha
D\right)^{2}+G^{2}},$ (4)
$\displaystyle\mbox{\boldmath$v$}_{y}=\frac{GD\left(\alpha-\beta\right)\mbox{\boldmath$j$}_{x}+\left(\alpha\beta
D^{2}+G^{2}\right)\mbox{\boldmath$j$}_{y}}{\left(\alpha D\right)^{2}+G^{2}}.$
(5)

If we now solve for two orthogonal currents,
$\mbox{\boldmath$j$}_{1}=\left(\mbox{\boldmath$j$},0\right)$ and
$\mbox{\boldmath$j$}_{2}=\left(0,\mbox{\boldmath$j$}\right)$, it is clear that
the product of the resulting velocity vectors is zero.

We know that above the de-pinning threshold the current velocity relationship
of skyrmions is close to linear [22, 23] and that if we reverse the the
current flow we reverse the direction of motion [24]. We therefore measure the
skyrmion velocity for two different current magnitudes, here $3\times
10^{11}$A m-2 and $6\times 10^{11}$A m-2, in both the $e_{x}$ and $e_{y}$
directions and derive unit vectors and average velocity per $1\times 10^{10}$A
m-2. From these values, assuming linearity at the current levels applied, we
can calculate a current magnitude and time required in each direction of
current application, $e_{x}$ and $e_{y}$. The currents are then applied
divided by the time-steps taken, here we switch the current at 200 MHz [8], to
arrive at a prescription for moving a skyrmion too a specific point in the
plane. As the deflection caused by the skyrmion Hall effect is included in the
derived vectors it can effectively be ignored. This process will work for any
topological degree, Eq. (1), of skyrmion including skyrmion bags [25, 26] and
other higher degree skyrmions as the variation in Hall deflection makes no
difference for this closed system.

As the currents are pulsed in two separate directions the maximum attainable
velocity is approximately half of the velocity attainable with current pulsed
in a single direction. This can be reduced if steps are taken to optimise the
velocity by increasing pulse magnitude and decreasing pulse time for the
shorter vector direction however we do not explore this here.

We use the graphics processing unit (GPU) enabled micromagnetic solver Mumax3
[27] to test our model by applying current in the cardinal and intercardinal
directions, $\frac{2\pi}{n}|n\in\left\\{1...8\right\\}$, and for distances in
the range $\left[2,40\right]$ nm in steps of $2$ nm. Simulation parameters are
given in the appendix. We find that the predicted directions vary by less than
2% and the distances by less than 1%.

##  III Skyrmion configurable devices

To demonstrate the effectiveness of this paradigm we simulate two different
basic devices that use 2D skyrmion motion. These devices are designed such
that they could be built into digital or analogue architectures as the primary
required inputs are current magnitude and pulse duration and the outputs are
current from the activation of magnetic tunnel junction (MTJ) switches. MTJ’s
can detect quantised changes in resistivity due to the presence of skyrmions
[28] and are considered ideal for the electrical detection of skyrmions. The
ability to be deployed in analogue or digital architectures is important for
use in neuromorphic applications and in solely analogue for use in high
radiation environments such as fusion reactors and in space. We call them
configurable as both the SCLG and the SCNN could be configured in the field to
fulfil a number of roles.

Figure 2: A field configurable skyrmion logic gate. (a) A simulated logic gate
with a skyrmion at its centre. (b) Direction vectors derived from the
calculated basis vectors. (c) Logic table for current inputs and NAND and NOR
gate outputs. NAND and NOR where chosen as they are the universal logic gates
and can construct all others.

The SCLG can be preconfigured with a given logic gate configuration with each
current pairing giving a different vector to a distinct MTJ, see Fig. 2, and
each MTJ configured to provide a current response. We illustrate this in Fig.
2 (c), providing examples of NAND and NOR, the universal logic gates. This is
a very simple device, theoretically ideal for production en masse in a large
scale field programmable processor.

Figure 3: A field configurable skyrmion neural network node. Currents are
accumulated from $6\times 10^{11}$A m-2 to MTJ activation at between $13.1$
and $14.3\times 10^{11}$A m-2 over a period of $160$ ns at zero external
field. (a) The progress of the skyrmion to activation over the $160$ ns time
period. (b) A close up of the skyrmion entering the MTJ disc. (c) We plot the
distance from the device center to the MTJ disc over time. Note that the
skyrmion comes to approximate rest between later current increases prior to
activation inside the yellow highlighted region.

The SCNN, see Fig. 3, is designed in a similar fashion to the SCLG and
features four arms with MTJ detectors. This design requires additional
external control in the form of an accumulator for the input currents, see
Fig. 1 (e). The skyrmion is driven by the accumulated currents from many pre-
synaptic SCNN’s until sufficient current is present to drive the skyrmion
through the edge repulsion barrier that is the constriction at the end of the
arm prior to the MTJ disc. This simulates the leaky–integrate-fire (LIF)
spiking model of neurons [29, 30] . Only one arm of the device can be utilised
during any one activation sequence and each arm can provide a different
activation profile by tuning the externally applied field, see Fig. 1 (c).

Both devices feature MTJ’s at their centres in order to detect the base state.
In both cases the skyrmion can be driven, by an opposing applied pulse
current, back to that base state. In the case of the SCNN this must be done to
push the skyrmion from the MTJ activation site at the end of the arm as the
constriction will not allow free passage. Once the skyrmion is free from the
MTJ disc in the SCNN case and after activation in the SCLG case the skyrmion
will return to the center of each device with no further applied current as
this is the lowest energy state, this is however much slower than driving the
skyrmion to the center. In both cases the returning current is lower than the
activation current as the edge repulsion acts to drive the skyrmion toward the
lowest energy state.

##  IV Discussion

We have shown that skyrmions can be driven with fine control in 2D. This extra
dimension allows us to reconsider the design of many proposed skyrmionic
devices. We use this 2D control to demonstrate skyrmion based configurable
logic gates that could then be deployed in more complex information processing
circuits. We also showed a configurable neural network node device that could
have it’s firing activation functions altered and could simulate four separate
neurons at distinct times. These devices where chosen so that a complete next
generation computing architecture could be considered, drawing from classical
computing and neuron inspired brain-like computing.

One motivation behind this work was the authors attempt to drive skyrmions
down arbitrarily shaped tree like structures that could be used to represent
the dendritic arbors of a typical human neuron. In conclusion we hope that the
shift from 1D to 2D skyrmionic motion will enable a new perspective in the
application of skyrmions to modern and future data, logic and brain-like
applications.

##  V Appendix

We consider skyrmionic structures in thin film multilayers where the
interfacial DMI results in Néel type skyrmions. The simulations where
performed using the GPU-accelerated micromagnetic simulation program MuMax3
[27] with Landau-Lifshitz dynamics in the form

$\displaystyle\frac{\partial\mbox{\boldmath$n$}}{\partial
t}=\gamma\frac{1}{1+\alpha^{2}}\left(\mbox{\boldmath$n$}\times\mbox{\boldmath$B$}_{\mathrm{eff}}+\alpha\mbox{\boldmath$n$}\times\left(\mbox{\boldmath$n$}\times\mbox{\boldmath$B$}_{\mathrm{eff}}\right)\right),$
(6)

where $\gamma\approx 176$ rad(ns T)-1 is the electron gyromagnetic ratio,
$\alpha=0.3$ the dimensionless damping parameter,
$\mbox{\boldmath$B$}_{\mathrm{eff}}=\delta E/\delta\mbox{\boldmath$n$}$ the
effective field and
$\mbox{\boldmath$n$}(\mbox{\boldmath$x$})=\mbox{\boldmath$N$}(\mbox{\boldmath$x$})/M_{s}$
the magnetisation vector field normalised by the saturation magnetisation.

The simulations external geometry is a $512\times 512$ nm2 rectangle of $1$ nm
thickness, in order to represent a typical surface that could be fabricated
using lithographic processing. Cell size of $2\times 2\times 1$ nm3 has been
used. Material parameters are: saturation magnetisation $M_{s}=580$ kAm-1,
exchange $20$ pJm-1, interfacial DMI $3$ mJm-2 and uniaxial anisotropy along
the $+z$ direction $0.53$ MJm-3 [31].

The applied current is of Zhang-Li [16] (CIP) type, the non-adiabacity of spin
transfer torque constant is set to $\xi=0.05$, current flows in the
$\mbox{\boldmath$e$}_{x}$ and $\mbox{\boldmath$e$}_{y}$ directions at separate
times:

$\displaystyle\overrightarrow{\tau}_{ZL}$
$\displaystyle=\frac{1}{1+\alpha^{2}}\left(\left(1+\xi\alpha\right)\mbox{\boldmath$n$}\times\left(\mbox{\boldmath$n$}\times\left(\mbox{\boldmath$u$}\cdot\nabla\right)\right)\mbox{\boldmath$n$}\right.$
$\displaystyle+\left.\left(\xi-\alpha\right)\mbox{\boldmath$n$}\times\left(\mbox{\boldmath$u$}\cdot\nabla\right)\mbox{\boldmath$n$}\right)$
(7) $u$
$\displaystyle=\frac{\mu_{B}\mu_{0}}{2e\gamma_{0}B_{sat}\left(1+\xi^{2}\right)}\mbox{\boldmath$j$}$
(8)

where $j$ A m-2 is the current density, $\mu_{B}$ the Bohr magneton and
$B_{sat}$ the saturation magnetisation expressed in Tesla.

## Acknowledgements

This work was supported by the Engineering and Physical Sciences Research
Council (EPSRC) grant EP/M506473/1. The Titan V GPU used for parts of this
research was donated by the NVIDIA Corporation.

## References

  * Bogdanov and Hubert [1994] A. Bogdanov and A. Hubert, Thermodynamically stable magnetic vortex states in magnetic crystals, Journal of Magnetism and Magnetic Materials 138, 255 (1994).
  * Fert  _et al._ [2013] A. Fert, V. Cros, and J. Sampaio, Skyrmions on the track, Nature Nanotechnology 8, 152 (2013).
  * Zhang  _et al._ [2015] X. Zhang, M. Ezawa, and Y. Zhou, Magnetic skyrmion logic gates: Conversion, duplication and merging of skyrmions, Scientific Reports 5, 9400 (2015).
  * Chen  _et al._ [2020] R. Chen, Y. Li, V. F. Pavlidis, and C. Moutafis, Skyrmionic interconnect device, Physical Review Research 2, 043312 (2020).
  * Song  _et al._ [2020] K. M. Song, J.-S. Jeong, B. Pan, X. Zhang, J. Xia, S. Cha, T.-E. Park, K. Kim, S. Finizio, J. Raabe, J. Chang, Y. Zhou, W. Zhao, W. Kang, H. Ju, and S. Woo, Skyrmion-based artificial synapses for neuromorphic computing, Nature Electronics 3, 148 (2020).
  * Grollier  _et al._ [2020] J. Grollier, D. Querlioz, K. Y. Camsari, K. Everschor-Sitte, S. Fukami, and M. D. Stiles, Neuromorphic spintronics, Nature Electronics 3, 360 (2020).
  * Parkin  _et al._ [2008] S. S. P. Parkin, M. Hayashi, and L. Thomas, Magnetic Domain-Wall Racetrack Memory, Science 320, 190 (2008).
  * Woo  _et al._ [2018] S. Woo, K. M. Song, X. Zhang, Y. Zhou, M. Ezawa, X. Liu, S. Finizio, J. Raabe, N. J. Lee, S.-I. Kim, S.-Y. Park, Y. Kim, J.-Y. Kim, D. Lee, O. Lee, J. W. Choi, B.-C. Min, H. C. Koo, and J. Chang, Current-driven dynamics and inhibition of the skyrmion Hall effect of ferrimagnetic skyrmions in GdFeCo films, Nature Communications 9, 959 (2018).
  * Zhang  _et al._ [2016] X. Zhang, Y. Zhou, and M. Ezawa, Antiferromagnetic Skyrmion: Stability, Creation and Manipulation, Scientific Reports 6, 24795 (2016).
  * Jiang  _et al._ [2017] W. Jiang, X. Zhang, G. Yu, W. Zhang, X. Wang, M. Benjamin Jungfleisch, J. E. Pearson, X. Cheng, O. Heinonen, K. L. Wang, Y. Zhou, A. Hoffmann, and S. G. E. te Velthuis, Direct observation of the skyrmion Hall effect, Nature Physics 13, 162 (2017).
  * Gnoli  _et al._ [2021] L. Gnoli, F. Riente, M. Vacca, M. Ruo Roch, and M. Graziano, Skyrmion Logic-In-Memory Architecture for Maximum/Minimum Search, Electronics 10, 155 (2021).
  * Thiele [1973] A. A. Thiele, Steady-State Motion of Magnetic Domains, Physical Review Letters 30, 230 (1973).
  * Dzyaloshinsky [1958] I. Dzyaloshinsky, A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics, Journal of Physics and Chemistry of Solids 4, 241 (1958).
  * Moriya [1960] T. Moriya, Anisotropic Superexchange Interaction and Weak Ferromagnetism, Physical Review 120, 91 (1960).
  * Chen  _et al._ [2018] X. Chen, W. Kang, D. Zhu, X. Zhang, N. Lei, Y. Zhang, Y. Zhou, and W. Zhao, A compact skyrmionic leaky–integrate–fire spiking neuron device, Nanoscale 10, 6139 (2018).
  * Zhang and Li [2004] S. Zhang and Z. Li, Roles of Nonequilibrium Conduction Electrons on the Magnetization Dynamics of Ferromagnets, Physical Review Letters 93, 127204 (2004).
  * Gilbert [2004] T. Gilbert, A phenomenological theory of damping in ferromagnetic materials, IEEE Transactions on Magnetics 40, 3443 (2004).
  * Thiaville  _et al._ [2005] A. Thiaville, Y. Nakatani, J. Miltat, and Y. Suzuki, Micromagnetic understanding of current-driven domain wall motion in patterned nanowires, EPL (Europhysics Letters) 69, 990 (2005).
  * Everschor  _et al._ [2011] K. Everschor, M. Garst, R. A. Duine, and A. Rosch, Current-induced rotational torques in the skyrmion lattice phase of chiral magnets, Physical Review B 84, 064401 (2011).
  * Elías  _et al._ [2017] R. G. Elías, N. Vidal-Silva, and A. Manchon, Steady motion of skyrmions and domains walls under diffusive spin torques, Physical Review B 95, 104406 (2017).
  * Huber [1982] D. L. Huber, Dynamics of spin vortices in two-dimensional planar magnets, Physical Review B 26, 3758 (1982).
  * Iwasaki  _et al._ [2013] J. Iwasaki, M. Mochizuki, and N. Nagaosa, Universal current-velocity relation of skyrmion motion in chiral magnets, Nature Communications 4, 1463 (2013).
  * Tomasello  _et al._ [2014] R. Tomasello, E. Martinez, R. Zivieri, L. Torres, M. Carpentieri, and G. Finocchio, A strategy for the design of skyrmion racetrack memories, Scientific Reports 4, 6784 (2014).
  * Wang  _et al._ [2022] W. Wang, D. Song, W. Wei, P. Nan, S. Zhang, B. Ge, M. Tian, J. Zang, and H. Du, Electrical manipulation of skyrmions in a chiral magnet, Nature Communications 13, 1593 (2022).
  * Foster  _et al._ [2019] D. Foster, C. Kind, P. J. Ackerman, J.-S. B. Tai, M. R. Dennis, and I. I. Smalyukh, Two-dimensional skyrmion bags in liquid crystals and ferromagnets, Nature Physics 15, 655 (2019).
  * Kind and Foster [2021] C. Kind and D. Foster, Magnetic skyrmion binning, Physical Review B 103, L100413 (2021).
  * Vansteenkiste  _et al._ [2014] A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. Van Waeyenberge, The design and verification of MuMax3, AIP Advances 4, 107133 (2014).
  * Guang  _et al._ [2023] Y. Guang, L. Zhang, J. Zhang, Y. Wang, Y. Zhao, R. Tomasello, S. Zhang, B. He, J. Li, Y. Liu, J. Feng, H. Wei, M. Carpentieri, Z. Hou, J. Liu, Y. Peng, Z. Zeng, G. Finocchio, X. Zhang, J. M. D. Coey, X. Han, and G. Yu, Electrical Detection of Magnetic Skyrmions in a Magnetic Tunnel Junction, Advanced Electronic Materials 9, 2200570 (2023).
  * Abbott [1999] L. F. Abbott, Lapicque’s introduction of the integrate-and-fire model neuron (1907), Brain Research Bulletin 50, 303 (1999).
  * Hodgkin and Huxley [1990] A. L. Hodgkin and A. F. Huxley, A quantitative description of membrane current and its application to conduction and excitation in nerve, Bulletin of Mathematical Biology 52, 25 (1990).
  * Kind  _et al._ [2020] C. Kind, S. Friedemann, and D. Read, Existence and stability of skyrmion bags in thin magnetic films, Applied Physics Letters 116, 022413 (2020).


∎

11institutetext: Zoltán Prekopcsák 22institutetext: Budapest University of
Technology and Economics, Hungary
Magyar tudósok körútja 2, Budapest, H-1117 Hungary
22email<EMAIL_ADDRESS>Phone: +36 1 4633119, Fax: +36 1 4633107
33institutetext: Daniel Lemire 44institutetext: LICEF, Université du Québec à
Montréal (UQAM)
100 Sherbrooke West, Montreal, QC, H2X 3P2 Canada

# Time Series Classification by Class-Based Mahalanobis Distances

Zoltán Prekopcsák  Daniel Lemire

(Received: date / Accepted: date)

###### Abstract

To classify time series by nearest neighbors, we need to specify or learn one
or several distances. We consider variations of the Mahalanobis distances
which rely on the inverse covariance matrix of the data. Unfortunately – for
time series data – the covariance matrix has often low rank. To alleviate this
problem we can either use a pseudoinverse, covariance shrinking or limit the
matrix to its diagonal. We review these alternatives and benchmark them
against competitive methods such as the related Large Margin Nearest Neighbor
Classification (LMNN) and the Dynamic Time Warping (DTW) distance. As we
expected, we find that the DTW is superior, but the Mahalanobis distances are
computationally inexpensive in comparison. To get best results with
Mahalanobis distances, we recommend learning one distance per class using
either covariance shrinking or the diagonal approach.

###### Keywords:

Time-series classification Distance learning Nearest Neighbor Mahalanobis
distance

###### MSC:

62-07 62H30

††journal: Advances in Data Analysis and Classification

##  1 Introduction

Time series are sequences of values measured over time. Examples include
financial data, such as stock prices, or medical data, such as blood sugar
levels. Classifying time series is an important class of problems which is
applicable to music classification (Weihs et al, 2007), medical diagnostic
(Sternickel, 2002) or bioinformatics (Legrand et al, 2008).

Nearest Neighbor (NN) methods classify time series efficiently and accurately
(Ding et al, 2008). The 1-NN method is especially simple: we merely have to
find the nearest labeled instance.

We need to specify a distance: the Euclidean and Dynamic Time Warping (Sakoe
and Chiba, 1978a) distances are popular choices. However, we can also learn a
distance based on some training data (Yang and Jin, 2006; Weinberger and Saul,
2009). Given the training data set made of classes of time series instances,
we can either learn a single (global) distance function, or learn one distance
function per class (Csatári and Prekopcsák, 2010; Paredes and Vidal, 2000,
2006). That is, to compute the distance between a test element and an instance
of class $j$, we use a distance function specific to class $j$.

Because the Euclidean distance is popular for NN classification, it is
tempting to consider generalized ellipsoid distances (Ishikawa et al, 1998),
that is, distances of the form

$\displaystyle D(x,y)=(x-y)^{\top}M(x-y).$

When $M$ is a positive semi-definite matrix, the square root of this distance
is a pseudometric: it is symmetric, non-negative and it satisfies the triangle
inequality ($\sqrt{D(x,y)}+\sqrt{D(y,z)}\leq\sqrt{D(x,z)}$). Further, when $M$
is a positive definite matrix, then the square root becomes a metric because
$D(x,x)=0\Rightarrow x=0$. When the matrix $M$ is the identity matrix, we
recover the (squared) Euclidean distance. We get the Mahalanobis distance when
solving for the matrix $M$ minimizing the sum of distances $\sum_{x,y}D(x,y)$
(see § 3). We can require $M$ to be diagonal, thus defining the diagonal
Mahalanobis distance (Paredes and Vidal, 2006; Ishikawa et al, 1998). In the
more general case, solving for $M$ can be difficult, as it often involves
inverting a low-rank matrix. Perhaps partly due to these mathematical
difficulties, there has been no attempt to use the general Mahalanobis
distance to classify time series (to our knowledge). Thus, for the first time,
we apply the full-matrix Mahalanobis distance for time series classification.
To solve the mathematical difficulties, we use both an approach based on a
pseudoinverse and on covariance shrinkage. With both the diagonal and the
covariance shrinkage approaches, the square root of the distance $D$ is a
metric (under mild assumptions) whereas our approach based on a pseudoinverse
merely generates a pseudo-metric.

While we find that the full-matrix Mahalanobis distance is not competitive
when relying on a pseudoinverse, we get good results with covariance shrinkage
or using the diagonal Mahalanobis distance. Moreover, we find that the class-
based Mahalanobis distance is preferable to the global Mahalanobis distance.

##  2 Related Works

Several distance functions are used for time series classification, such as

  * •

Dynamic Time Warping (DTW) (Ratanamahatana and Keogh, 2004),

  * •

DISSIM (Frentzos et al, 2007),

  * •

Threshold Queries (Aßfalg et al, 2006),

  * •

Edit distances (Chen and Ng, 2004; Chen et al, 2005),

  * •

Longest Common Subsequences (LCSS) (Vlachos et al, 2002),

  * •

Swale (Morse and Patel, 2007),

  * •

SpADe (Chen et al, 2007),

  * •

and Cluster, Then Classify (CTC) (Keogh and Pazzani, 1998).

Ding et al (2008) presented an extensive comparison of these distance
functions and concluded that DTW is among the best measures and that the
accuracy of the Euclidean distance converges to DTW as the size of the
training set increases.

In a general Machine Learning setting, Paredes and Vidal (2000, 2006) compared
Euclidean distance with the conventional and class-based Mahalanobis
distances. One of our contribution is to validate these generic results on
time series: instead of tens of features, we have hundreds or even thousands
of values which makes the problem mathematically more challenging: the rank of
our covariance matrices are often tiny compared to their sizes.

More generally, distance metric learning has an extensive literature
(Wettschereck et al, 1997; Hastie and Tibshirani, 1996; Chai et al, 2010;
Short and Fukunaga, 1980). We refer the reader to Weinberger and Saul (2009)
for a review.

A conventional distance-learning approach is to find an optimal generalized
ellipsoid distance with respect to a specific loss function. The LMNN
algorithm proposed by Weinberger and Saul (2009) takes a different approach.
It seeks to force nearest neighbors to belong to the same class and it
separates instances from different classes by a large margin. LMNN can be
formulated as a semi-definite programming problem. They also propose a
modification which they call multiple metrics LMNN as it learns different
distances for each class.

There are many extensions and alternatives to NN classification. For example,
Jahromi et al (2009) use instance weights to improve classification.
Meanwhile, Zhan et al (2009) learn a distance per instance.

###  2.1 Dynamic Time Warping (DTW)

As already stated, The DTW is one of the most accurate distance function for
time-series classification. The DTW was invented to recognize spoken words
(Sakoe and Chiba, 1978b), but it is also used for various problems such as
handwriting recognition (Bahlmann, 2004; Niels and Vuurpijl, 2005), chromosome
classification (Legrand et al, 2007), networking (Ang et al, 2010) or shape
retrieval (Bartolini et al, 2005; Marzal et al, 2006).

For simplicity, consider two time series $x$ and $y$ of equal length ($n$).
Recall that the $l_{p}$ distance between two time series is

$\displaystyle\sqrt[p]{\sum_{i=1}^{n}|x_{i}-y_{i}|^{p}}$

if $p$ is a positive integer or $\max_{i=1,\ldots,n}|x_{i}-y_{i}|$ if
$p=\infty$. The $l_{1}$ distance is also called the Manhattan distance whereas
the $l_{2}$ distance is the Euclidean distance.

The intuition behind the DTW is that one could speak the same sentence by
speeding up and slowing down without changing the meaning of the sounds. The
DTW is a generalization of the $l_{p}$ distance which allows the data to be
realigned. To compute the DTW between $x$ and $y$, you must find a many-to-
many matching between the data points in $x$ and the data points in $y$. That
is each data point from one series must be matched with at least one data
point with the other series. One such matching is the trivial one, which maps
the first data point from $x$ to the first data point in $y$, the second data
point in $x$ to the second data point in $y$, and so on. Write the set matches
$(i,j)$ as $\Gamma$ so that the trivial matching is just
$\Gamma=\\{(1,1),(2,2),\ldots,(n,n)\\}$. The $l_{p}$ cost corresponding to a
matching is defined as

$\displaystyle\sqrt[p]{\sum_{(i,j)\in\Gamma}|x_{i}-y_{j}|^{p}}$

if $p$ is in an integer or $\max_{(i,j)\in\Gamma}|x_{i}-y_{j}|$ if $p=\infty$.
Typically, $p$ is either 1 or 2: for our purposes we assume $p=2$. For a given
$p$, the DTW is defined as the minimal cost over all allowed matchings
$\Gamma$. Typically, we require matchings to be monotonic: if both $(i,j)$ and
$(i+1,j^{\prime})$ are in $\Gamma$ then $j^{\prime}\geq j$, that is, we cannot
warp back in time. Moreover, some matches might be forbidden, maybe because
the data points are too far apart (Itakura, 1975; Sakoe and Chiba, 1978b). Yu
et al (2011) has proposed learning this warping constraint from the data.
Except when $p=\infty$, the DTW fails to satisfy the triangle inequality: the
DTW is not a metric distance.

The computational cost of the DTW is sometimes a challenge (Salvador and Chan,
2007). To alleviate this problem, several strategies have been proposed
including lower bounds and R*-tree indexes (Ratanamahatana and Keogh, 2005;
Lemire, 2009; Ouyang and Zhang, 2010).

Gaudin and Nicoloyannis (2006) proposed a weighted version of the DTW called
Adaptable Time Warping. Instead of computing
$\sum_{(i,j)\in\Gamma}|x_{i}-y_{j}|^{p}$, it computes
$\sum_{(i,j)\in\Gamma}M_{i,j}|x_{i}-y_{j}|^{p}$ where $M$ is some matrix.
Unfortunately, finding the optimal matrix $M$ can be a challenge. Jeong et al
(2010) investigated another form of weighted DTW where you seek the minimize

$\displaystyle\sqrt[p]{\sum_{(i,j)\in\Gamma}w_{|i-j|}|x_{i}-y_{j}|^{p}}$

where $w$ is some weight vector. Many other variations on the DTW distance
have been proposed, e.g., Chouakria and Nagabhushan (2007).

##  3 Mahalanobis distance

For completeness, we derive the Mahalanobis distance (Mahalanobis, 1936) as an
optimal form of generalized ellipsoid distance. We seek $M$ minimizing

$\displaystyle\sum_{x,y\in S}{(x-y)^{\top}M(x-y)}=\sum_{x,y\in
S}{\left(\sum_{k=1}^{n}\sum_{l=1}^{n}{(x_{k}-y_{k})m_{kl}(x_{l}-y_{l})}\right)}$

where $S$ is some class of time series. We add a regularization constraint on
the determinant ($\det(M)=1$).

We solve the minimization problem by the Lagrange’s multiplier method with the
Lagrangian

$\displaystyle L(M,\lambda)$ $\displaystyle=$ $\displaystyle\sum_{x,y\in
S}\left(\sum_{k=1}^{n}\sum_{l=1}^{n}{(x_{k}-y_{k})m_{k,l}(x_{l}-y_{l})}\right)-\lambda(\det
M-1).$

We want to compute the derivative of the Lagrangian, and of $\det(M)$, with
respect to $m_{k,l}$. By Laplace expansion, we have for all $l$ that

$\displaystyle\det(M)=\sum_{k=1}^{n}{(-1)^{k+l}m_{k,l}\det(M_{k,l})}=1$

where $M_{k,l}$ is the $(k,l)$ minor of $M$: an $(n-1)\times(n-1)$ matrix
obtained by deleting $k$-th row and $l$-th column of $M$. Thus, we have
${\partial\det(M)}/{\partial m_{k,l}}=(-1)^{k+l}\det(M_{k,l})$. Setting the
derivatives to zero, we get

$\displaystyle\frac{\partial L(M,\lambda)}{\partial m_{k,l}}=\sum_{x,y\in
S}{(x_{k}-y_{k})(x_{l}-y_{l})}-\lambda(-1)^{k+l}\det(M_{k,l})=0$

and therefore

$\displaystyle\det(M_{k,l})=\frac{\sum_{x,y\in
S}{(x_{k}-y_{k})(x_{l}-y_{l})}}{\lambda(-1)^{k+l}}.$

Because $\det(M)=1$ and using Cramer’s rule, the inverse matrix $M^{-1}$ can
be represented as

$\displaystyle
m_{k,l}^{-1}=\frac{(-1)^{k+l}\det(M_{k,l})}{\det(M)}=(-1)^{k+l}\det(M_{k,l}).$

Hence, we have

$\displaystyle m_{k,l}^{-1}=\frac{\sum_{x,y\in
S}{(x_{k}-y_{k})(x_{l}-y_{l})}}{\lambda}.$

Thus, we have that $M^{-1}\propto C$ where $C$ is the covariance matrix.
Because we require $M$ to be positive definite and to satisfy $\det(M)=1$, we
set $M=(\det(C))^{\frac{1}{n}}C^{-1}$ which produces the Mahalanobis
distance111The original Mahalanobis distance is defined with $M=C^{-1}$.

Our derivation assumes that the covariance matrix is (numerically) invertible.
This fails in practice. In § 4, we review some solutions.

##  4 Computing Mahalanobis distances for time series

As a rule of thumb, the covariance matrix becomes singular when the number of
instances is smaller or about the same as the number of attributes. This is a
common problem with time series: whereas individual time series might have
thousands of samples, there may only be a few labeled time series in each
class.

The most straight-forward solution is to limit the matrix $M$ to its
diagonal–thus producing a weighted Euclidean distance. Revisiting the
derivation of § 3 where we require $m_{k,l}=0$ for $k\neq l$, we get that the
inverse of the Mahalanobis matrix $M$ must be equal to the inverse of the
diagonal of the covariance matrix: $M^{-1}\propto\mathrm{diag}(C)$. As long as
the variance of each attribute in our training sets is different from zero – a
condition satisfied in practice in our experiments, the problem is well posed
and the result is a positive-definite matrix. In such a diagonal case, the
number of parameters to learn grows only linearly with the number of
attributes in the time series. In contrast, the number of elements in the full
covariance matrix grows quadratically. The speed of the computation of the
distances also depends on the number of non-zero elements in the Mahalanobis
matrix $M$.

Alas, the diagonal Mahalanobis distance fails to take into account the
information off the diagonal in the covariance matrix. See Figure 1 for the
covariance matrix of a class of time series. It is clear from the figure that
the covariance matrix has significant values off the diagonal. There are even
block-like patterns in the matrix corresponding to specific time intervals.

(a) Sample of time series (b) Sample covariance

Figure 1:  The Cylinder class from the CBF data set and its sample covariance.
Higher absolute values in the matrix are presented using darker colors.

Could it be that non-diagonal Mahalanobis distance could be superior or at
least competitive with the diagonal Mahalanobis distance? It is tempting to
use banded matrices, but the restriction of a positive definite matrix to a
band may fail to be positive definite. Block-diagonal matrices (Matton et al,
2010) can preserve positive definiteness, but learning which blocks to use in
the context of time series might be difficult. Instead, we propose two
approaches: one is based on the widely used Moore-Penrose pseudoinverse, and
the other is covariance shrinkage. See Figure 2 for the three different
covariance estimates of the same class.

(a) Sample covariance (b) Shrinked covariance (c) Diagonal covariance

Figure 2:  The covariance estimates of the Funnel class in the CBF data set.
Large absolute values are in darker colors. Both the shrinked and diagonal
covariances are positive definite whereas the sample covariance matrix is
singular.

The approach based on the pseudoinverse is based on the singular value
decomposition (SVD). We write the SVD as $C=U\Sigma V^{*}$ where $\Sigma$ is a
diagonal matrix with eigenvalues $\lambda_{1},\lambda_{2},\ldots$ and $U$ and
$V$ are unitary matrices. The Moore-Penrose pseudoinverse is given by
$V\Sigma^{+}U^{*}$ where $\Sigma^{+}$ is the diagonal matrix made of the
eigenvalues $1/\lambda_{1},1/\lambda_{2},\ldots$ with the convention that
$1/0=0$. The pseudo-determinant is the product of the non-zero eigenvalues of
$\Sigma$. We set $M$ equal to the pseudoinverse of the covariance
matrix—normalized so that it has a pseudo-determinant of one. This solution is
equivalent to projecting the time series data on the subspace corresponding to
the non-zero eigenvalues of $\Sigma$. That is, the matrix $M$ will be
singular.

Covariance shrinkage is an estimation method for problems with small number of
instances and large number of attributes (Stein, 1956). It has better
theoretical and practical properties for such data sets as the estimated
covariance matrix is guaranteed to be non-singular. Let
$x_{1},x_{2},\ldots,x_{n}$ be a set of time series. We write
$x_{i}=(x_{i1},x_{i2},\ldots,x_{ip})$ where $x_{ik}$ is the $k^{\mathrm{th}}$
value (attribute) of the time series $x_{i}$. One element of the (sample)
covariance matrix $S$ is

$\displaystyle
s_{ij}=\frac{1}{n-1}\sum_{k=1}^{n}(x_{ki}-\bar{x}_{i})(x_{kj}-\bar{x}_{j})$

where $\bar{x}_{i}=\frac{1}{n}\sum_{k=1}^{n}x_{ki}$. It is positive semi-
definite but can be singular. To fix the problem that $S$ is singular, we can
replace it with an estimation of the form

$\displaystyle C^{*}=\lambda T+(1-\lambda)S$

for some suitably chosen target $T$. If $T$ is a positive definite matrix and
$\lambda\in(0,1]$, we have that $\lambda T+(1-\lambda)S$ must be positive
definite. Moreover, the smallest eigenvalue of $\lambda T+(1-\lambda)S$ must
be at least as large as $\lambda$ times the smallest eigenvalue of $T$. We
have used the target recommended by Schäfer and Strimmer (2005) which is the
diagonal of the unrestricted covariance estimate, $T=\mathrm{diag}(C)$. It is
positive definite in our examples. For $\lambda$, we use the estimation
proposed by Schäfer and Strimmer (2005): it is computationally inexpensive.

Thus, finally, we consider six types of Mahalanobis distances: two localities
(global or class-based) and three estimators (pseudoinverse, shrinkage, or
diagonal).

##  5 Experiments

The main goal of our experiments is to evaluate Mahalanobis distances and the
class-based approach on time series. A secondary goal is to evaluate the LMNN
method.

We begin all tests with training data set made of classes of instances. When
applicable, distances are learned from this data set. We then attempt to
classify some test data using 1-NN and we report the classification error.

The code for the experiments is available online (Prekopcsák, 2011) with
instructions on how the results can be reproduced. For LMNN, we use the source
code provided by Weinberger and Saul (2008) for the experiments with default
parameters.

###  5.1 Data sets

We use the UCR time series classification benchmark (Keogh et al, 2006) for
our experiments as it includes diverse time series data sets from many
domains. It has predefined training-test splits for the experiments, so the
results can be compared across different papers. We removed the two data that
are not z-normalized by default (Beef and Coffee). Indeed, z-normalization
improves substantially the classification accuracy—irrespective of the chosen
distance. Thus, for fair results, we should z-normalize them, but this may
create confusion with previously reported numbers. We also removed the Wafer
data set as all distances classify it nearly perfectly. The remaining 17 data
sets were used for the comparison of different methods.

###  5.2 Best Mahalanobis distance for 1-NN accuracy

We compare the various Mahalanobis distances in Table 1. We have left out the
Moore-Penrose pseudoinverse, because its error rates were twice as high on
average compared to the other variants. What is immediately apparent is that
the class-based metrics give better classification results.

The diagonal Mahalanobis is somewhat better and they are also considerably
faster computationally, but the shrinkage estimate yields significantly better
results for several data sets (e.g. Adiac and Face (four)). Thus, out of the
six variations, we recommend the class-based shrinkage estimate and the class-
based diagonal Mahalanobis distance.

Table 1: Classification error for the various Mahalanobis distances. Data set | Shrinkage | Diagonal
---|---|---
global | class-based | global | class-based
50 words | 0.36 | 0.71 | 0.34 | 0.32
Adiac | 0.33 | 0.28 | 0.37 | 0.36
CBF | 0.52 | 0.04 | 0.16 | 0.05
ECG | 0.13 | 0.09 | 0.10 | 0.08
Fish | 0.33 | 0.15 | 0.19 | 0.18
Face (all) | 0.31 | 0.27 | 0.32 | 0.25
Face (four) | 0.45 | 0.10 | 0.16 | 0.17
Gun-Point | 0.06 | 0.10 | 0.10 | 0.11
Lighting-2 | 0.49 | 0.31 | 0.25 | 0.25
Lighting-7 | 0.59 | 0.32 | 0.36 | 0.23
OSU Leaf | 0.69 | 0.69 | 0.46 | 0.46
OliveOil | 0.17 | 0.17 | 0.17 | 0.13
Swedish Leaf | 0.25 | 0.14 | 0.21 | 0.18
Trace | 0.27 | 0.09 | 0.21 | 0.07
Two Patterns | 0.10 | 0.10 | 0.12 | 0.12
Synthetic Control | 0.23 | 0.10 | 0.13 | 0.09
Yoga | 0.24 | 0.21 | 0.17 | 0.17
# of best errors | 2 | 6 | 3 | 10

###  5.3 Comparing competitive distances

How do the class-based Mahalanobis distances fare compared to competitive
distances? Computationally, the diagonal Mahalanobis is inexpensive compared
to schemes such as the DTW or LMNN. Regarding the 1-NN classification error
rate, we give the results in Table 2. As expected (Ding et al, 2008), no
distance is better on all data sets. However, because the diagonal Mahalanobis
distance is closely related to the Euclidean distance, we compare their
classification accuracy. In two data sets, the Euclidean distance outperformed
the class-based Mahalanobis distance and only by small differences (0.09
versus 0.10-0.12). Meanwhile, the class-based diagonal Mahalanobis
outperformed the Euclidean distance 12 times, and sometimes by large margins
(0.07 versus 0.24 and 0.05 versus 0.15).

The LMNN is also competitive: its classification error is sometimes half that
of the Euclidean distance. However, the class-based LMNN gets the best result
among all methods only twice as opposed to five times for the global LMNN.
Moreover, the global LMNN significantly outperforms the class-based LMNN on
the Two Patterns data set (0.05 versus 0.24). For time series data sets, the
class-based LMNN is not an improvement over the global LMNN.

We have to note that DTW has the lowest error rates and provides best results
for half of the data sets, but it is much slower than Mahalanobis distances.

Table 2: Classification errors for some competitive schemes. We use class-based Mahalanobis distances. For the 50 words data set, the LMNN computation fails because it has a class with only one instance. Data set | Euclidean | DTW | C.-b. Mahalanobis | LMNN |
---|---|---|---|---|---
|  |  | shrink. | diag. |  | c.-b. |
50 words | 0.37 | 0.31 | 0.71 | 0.32 | — | — |
Adiac | 0.39 | 0.40 | 0.28 | 0.36 | 0.23 | 0.32 |
CBF | 0.15 | 0.00 | 0.04 | 0.05 | 0.15 | 0.15 |
ECG | 0.12 | 0.23 | 0.09 | 0.08 | 0.10 | 0.07 |
Fish | 0.22 | 0.17 | 0.15 | 0.18 | 0.13 | 0.14 |
Face (all) | 0.29 | 0.19 | 0.27 | 0.25 | 0.16 | 0.20 |
Face (four) | 0.22 | 0.17 | 0.10 | 0.17 | 0.16 | 0.16 |
Gun-Point | 0.09 | 0.09 | 0.10 | 0.11 | 0.05 | 0.09 |
Lighting-2 | 0.25 | 0.13 | 0.31 | 0.25 | 0.41 | 0.34 |
Lighting-7 | 0.42 | 0.27 | 0.32 | 0.23 | 0.51 | 0.48 |
OSU Leaf | 0.48 | 0.41 | 0.69 | 0.46 | 0.57 | 0.54 |
OliveOil | 0.13 | 0.13 | 0.17 | 0.13 | 0.13 | 0.13 |
Swedish Leaf | 0.21 | 0.21 | 0.14 | 0.18 | 0.21 | 0.19 |
Trace | 0.24 | 0.00 | 0.09 | 0.07 | 0.20 | 0.20 |
Two Patterns | 0.09 | 0.00 | 0.10 | 0.12 | 0.05 | 0.24 |
Synthetic Control | 0.12 | 0.01 | 0.10 | 0.09 | 0.03 | 0.09 |
Yoga | 0.17 | 0.16 | 0.21 | 0.17 | 0.18 | 0.18 |
# of best errors | 1 | 9 | 2 | 2 | 5 | 2 |

###  5.4 Effect of the number of instances per class

Whereas Table 2 shows that the Mahalanobis distances are far superior to the
Euclidean distance on some data sets, this result is linked to the number of
instances per class. For example, on the Wafer data set (which we removed),
there are many instances per class (500), and correspondingly, all distances
give a negligible classification error.

Thus, we considered three different synthetic time-series data sets with
varying numbers of instances per class: Cylinder-Bell-Funnel (CBF) (Saito,
1994), Control Charts (CC) (Pham and Chan, 1998) and Waveform (Breiman, 1998).
Test sets have 1 000 instances per class whereas training sets have between 10
to 1 000 instances. We repeated each test ten times, with different training
sets. Fig. 3 shows that whereas the class-based diagonal Mahalanobis is
superior to the Euclidean distance when there are few instances, this benefit
is less significant as the number of instances increases. Indeed, the
classification accuracy of the Euclidean distance grows closer to perfection
and it becomes more difficult for alternatives to be far superior.

Figure 3: Ratios of the 1-NN classification accuracies using the class-based
diagonal Mahalanobis and Euclidean distances

##  6 Conclusion

The Mahalanobis distances have received little attention for time series
classification and we are not surprised given their poor performance as a 1-NN
classifier when used in a straight-forward manner. However, by learning one
Mahalanobis distance per class we get a competitive classifier when using
either covariance shrinkage or a diagonal approach. Moreover, the diagonal
Mahalanobis distance is particularly appealing computationally: we only need
to compute the variances of the components. Meanwhile, we get good results
with the LMNN on time series data, though it is more expensive. The DTW is
superior, but computationally expensive.

###### Acknowledgements.

This work is supported by NSERC grant 261437.

## References

  * Ang et al (2010) Ang H, Gopalkrishnan V, Ng W, Hoi S (2010) On classifying drifting concepts in P2P networks. In: Balcázar J, Bonchi F, Gionis A, Sebag M (eds) Machine Learning and Knowledge Discovery in Databases, Lecture Notes in Computer Science, vol 6321, Springer Berlin / Heidelberg, pp 24–39
  * Aßfalg et al (2006) Aßfalg J, Kriegel H, Kröger P, Kunath P, Pryakhin A, Renz M (2006) Similarity search on time series based on threshold queries. EDBT 2006 pp 276–294
  * Bahlmann (2004) Bahlmann C (2004) The writer independent online handwriting recognition system frog on hand and cluster generative statistical Dynamic Time Warping. Writer 26(3):299–310
  * Bartolini et al (2005) Bartolini I, Ciaccia P, Patella M (2005) WARP: Accurate retrieval of shapes using phase of fourier descriptors and time warping distance. IEEE Transactions on Pattern Analysis and Machine Intelligence 27(1):142–147, DOI http://doi.ieeecomputersociety.org/10.1109/TPAMI.2005.21
  * Breiman (1998) Breiman L (1998) Classification and Regression Trees. Chapman & Hall/CRC
  * Chai et al (2010) Chai J, Liu H, Chen B, Bao Z (2010) Large margin nearest local mean classifier. Signal Process 90(1):236–248
  * Chen and Ng (2004) Chen L, Ng R (2004) On the marriage of $l_{p}$-norms and edit distance. In: VLDB’04, pp 1040–1049
  * Chen et al (2005) Chen L, Özsu M, Oria V (2005) Robust and fast similarity search for moving object trajectories. In: SIGMOD’05, ACM, pp 491–502
  * Chen et al (2007) Chen Y, Nascimento M, Ooi B, Tung A (2007) SpaDe: On shape-based pattern detection in streaming time series. In: ICDE’07, pp 786–795
  * Chouakria and Nagabhushan (2007) Chouakria A, Nagabhushan P (2007) Adaptive dissimilarity index for measuring time series proximity. Advances in Data Analysis and Classification 1:5–21
  * Csatári and Prekopcsák (2010) Csatári B, Prekopcsák Z (2010) Class-based attribute weighting for time series classification. In: POSTER 2010: Proceedings of the 14th International Student Conference on Electrical Engineering
  * Ding et al (2008) Ding H, Trajcevski G, Scheuermann P, Wang X, Keogh E (2008) Querying and mining of time series data: experimental comparison of representations and distance measures. In: VLDB’08, pp 1542–1552
  * Frentzos et al (2007) Frentzos E, Gratsias K, Theodoridis Y (2007) Index-based most similar trajectory search. In: ICDE 2007, pp 816–825
  * Gaudin and Nicoloyannis (2006) Gaudin R, Nicoloyannis N (2006) An adaptable time warping distance for time series learning. In: Proceedings of the 5th International Conference on Machine Learning and Applications, IEEE Computer Society, Washington, DC, USA, pp 213–218, DOI 10.1109/ICMLA.2006.12
  * Hastie and Tibshirani (1996) Hastie T, Tibshirani R (1996) Discriminant adaptive nearest neighbor classification. IEEE Transactions on Pattern Analysis and Machine Intelligence 18(6):607–616
  * Ishikawa et al (1998) Ishikawa Y, Subramanya R, Faloutsos C (1998) Mindreader: Querying databases through multiple examples. In: VLDB’98, pp 218–227
  * Itakura (1975) Itakura F (1975) Minimum prediction residual principle applied to speech recognition. IEEE Transactions on Acoustics, Speech, and Signal Processing 23(1):67–72
  * Jahromi et al (2009) Jahromi MZ, Parvinnia E, John R (2009) A method of learning weighted similarity function to improve the performance of nearest neighbor. Inf Sci 179(17):2964–2973
  * Jeong et al (2010) Jeong YS, Jeong MK, Omitaomu OA (2010) Weighted dynamic time warping for time series classification. Pattern Recognition In Press, Corrected Proof:–, DOI DOI:10.1016/j.patcog.2010.09.022
  * Keogh and Pazzani (1998) Keogh E, Pazzani M (1998) An enhanced representation of time series which allows fast and accurate classification, clustering and relevance feedback. In: Proceedings of the 4th International Conference of Knowledge Discovery and Data Mining, pp 239–241
  * Keogh et al (2006) Keogh E, Xi X, Wei L, Ratanamahatana CA (2006) The UCR time series classification/clustering homepage. http://www.cs.ucr.edu/~eamonn/time_series_data/ [last checked on 9/10/2010]
  * Legrand et al (2007) Legrand B, Chang CS, Ong SH, Neo SY, Palanisamy N (2007) Chromosome classification using dynamic time warping. Pattern Recognition Letters 29(3):215–222
  * Legrand et al (2008) Legrand B, Chang C, Ong S, Neo SY, Palanisamy N (2008) Chromosome classification using dynamic time warping. Pattern Recognition Letters 29(3):215 – 222
  * Lemire (2009) Lemire D (2009) Faster retrieval with a two-pass dynamic-time-warping lower bound. Pattern Recognition 42:2169–2180
  * Mahalanobis (1936) Mahalanobis PC (1936) On the generalised distance in statistics. In: Proceedings of the National Institute of Sciences of India, pp 49–55
  * Marzal et al (2006) Marzal A, Palazon V, Peris G (2006) Contour-based shape retrieval using Dynamic Time Warping. Lecture notes in Computer Science 4177:190
  * Matton et al (2010) Matton M, Compernolle DV, Cools R (2010) Minimum classification error training in example based speech and pattern recognition using sparse weight matrices. Journal of Computational and Applied Mathematics 234(4):1303–1311
  * Morse and Patel (2007) Morse MD, Patel JM (2007) An efficient and accurate method for evaluating time series similarity. In: SIGMOD ’07, pp 569–580
  * Niels and Vuurpijl (2005) Niels R, Vuurpijl L (2005) Using Dynamic Time Warping for intuitive handwriting recognition. In: IGS2005, pp 217–221
  * Ouyang and Zhang (2010) Ouyang Y, Zhang F (2010) Histogram distance for similarity search in large time series database. In: Proceedings of the 11th international conference on Intelligent data engineering and automated learning, Springer-Verlag, Berlin, Heidelberg, IDEAL’10, pp 170–177
  * Paredes and Vidal (2000) Paredes R, Vidal E (2000) A class-dependent weighted dissimilarity measure for nearest neighbor classification problems. Pattern Recognition Letters 21(12):1027–1036
  * Paredes and Vidal (2006) Paredes R, Vidal E (2006) Learning prototypes and distances: A prototype reduction technique based on nearest neighbor error minimization. Pattern Recognition 39(2):180–188
  * Pham and Chan (1998) Pham DT, Chan AB (1998) Control chart pattern recognition using a new type of self-organizing neural network. Proceedings of the Institution of Mechanical Engineers, Part I: Journal of Systems and Control Engineering 212(2):115–127
  * Prekopcsák (2011) Prekopcsák Z (2011) Matlab code for the experiments. http://github.com/Preko/Time-series-classification
  * Ratanamahatana and Keogh (2004) Ratanamahatana CA, Keogh E (2004) Making time-series classification more accurate using learned constraints. In: SDM’04, pp 11–22
  * Ratanamahatana and Keogh (2005) Ratanamahatana CA, Keogh E (2005) Three myths about Dynamic Time Warping data mining. In: SDM’05
  * Saito (1994) Saito N (1994) Local feature extraction and its applications using a library of bases. PhD thesis, Yale University, New Haven, CT, USA
  * Sakoe and Chiba (1978a) Sakoe H, Chiba S (1978a) Dynamic programming algorithm optimization for spoken word recognition. Acoustics, Speech and Signal Processing, IEEE Transactions on 26(1):43 – 49, DOI 10.1109/TASSP.1978.1163055
  * Sakoe and Chiba (1978b) Sakoe H, Chiba S (1978b) Dynamic programming algorithm optimization for spoken word recognition. IEEE Transactions on Acoustics, Speech, and Signal Processing 26(1):43–49
  * Salvador and Chan (2007) Salvador S, Chan P (2007) FastDTW: Toward accurate dynamic time warping in linear time and space. Intelligent Data Analysis 11(5):561–580
  * Schäfer and Strimmer (2005) Schäfer J, Strimmer K (2005) A shrinkage approach to large-scale covariance matrix estimation and implications for functional genomics. Statistical applications in genetics and molecular biology 4(1):32
  * Short and Fukunaga (1980) Short R, Fukunaga K (1980) A new nearest neighbor distance measure. In: Proceedings of the Fifth International Conference on Pattern Recognition, pp 81–86
  * Stein (1956) Stein C (1956) Inadmissibility of the usual estimator for the mean of a multivariate normal distribution. In: Proceedings of the Third Berkeley Symposium on Mathematical and Statistical Probability, pp 197––206
  * Sternickel (2002) Sternickel K (2002) Automatic pattern recognition in ECG time series. Computer Methods and Programs in Biomedicine 68(2):109 – 115
  * Vlachos et al (2002) Vlachos M, Gunopoulos D, Kollios G (2002) Discovering similar multidimensional trajectories. In: ICDE ’02, p 673
  * Weihs et al (2007) Weihs C, Ligges U, Mörchen F, Müllensiefen D (2007) Classification in music research. Advances in Data Analysis and Classification 1:255–291
  * Weinberger and Saul (2008) Weinberger K, Saul L (2008) Large margin nearest neighbor – matlab code. http://www.cse.wustl.edu/~kilian/Downloads/LMNN.html
  * Weinberger and Saul (2009) Weinberger K, Saul L (2009) Distance metric learning for large margin nearest neighbor classification. The Journal of Machine Learning Research 10:207–244
  * Wettschereck et al (1997) Wettschereck D, Aha DW, Mohri T (1997) A review and empirical evaluation of feature weighting methods for a class of lazy learning algorithms. Artif Intell Rev 11(1-5):273–314
  * Yang and Jin (2006) Yang L, Jin R (2006) Distance metric learning: A comprehensive survey. Tech. rep., Michigan State University
  * Yu et al (2011) Yu D, Yu X, Hu Q, Liu J, Wu A (2011) Dynamic time warping constraint learning for large margin nearest neighbor classification. Information Sciences 181(13):2787 – 2796, including Special Section on Databases and Software Engineering
  * Zhan et al (2009) Zhan DC, Li M, Li YF, Zhou ZH (2009) Learning instance specific distances using metric propagation. In: ICML ’09, pp 1225–1232