Dynamic system classifier

0.47 MB·

by  Daniel Pumpe

Dynamic system classifier
Daniel Pumpe∗1,2, Maksim Greiner1,2, Ewald Mu¨ller1,3 and Torsten A. Enßlin1,2
1 Max-Planck-Institut fu¨r Astrophysik, Karl-Schwarzschild-Str. 1, D-85748 Garching, Germany
2 Ludwigs-Maximilians-Universit¨at Mu¨nchen, Geschwister-Scholl-Platz 1, D-80539 Mu¨nchen, Germany
3 Technische-Universit¨at Mu¨nchen, Arcisstr. 21, D-80333 Mu¨nchen, Germany
(Dated: July, 21st 2016)
Stochastic differential equationsdescribe well manyphysical, biological and sociological systems,
despite the simplification often made in their derivation. Here the usage of simple stochastic dif-
ferential equations to characterize and classify complex dynamical systems is proposed within a
6 Bayesianframework. Tothisend,wedevelopadynamicsystem classifier (DSC).TheDSCfirstab-
1 stractstrainingdataofasystemintermsoftimedependentcoefficientsofthedescriptivestochastic
0 differential equation. Thereby the DSC identifies unique correlation structures within the training
2 data. For definiteness we restrict the presentation of DSC to oscillation processes with a time de-
l pendent frequency ω(t) and damping factor γ(t). Although real systems might be more complex,
u thissimpleoscillatorcapturesmanycharacteristicfeatures. Theωandγ timelinesrepresenttheab-
J
stract system characterization and permit the construction of efficient signal classifiers. Numerical

6 experimentsshow that such classifiers perform well even in the low signal-to-noise regime.
2
Keywords: DataAnalysis–SignalInference–BayesianMethod–ModelSelection–PatternRecog-

] nition
n
a
-
a I. INTRODUCTION ploring parts of the parameter space, showed that one
t canclassifytheirgravitationalwavesignaturesintothree
a
d A. Basic idea different main types [4]. This should permit to discrim-
. inate among these types, if from the few available nu-
s
c merical samples an abstract description can be learned,
i A classificationproblem for complex dynamic systems which can be confronted to real, noisy data. A way to
s
might look as follows: For a number of different system
y abstract training samples appropriately, permitting one
h classes a few training samples of the evolution of some to recognisegeneralcharacteristicsof time lines and fea-
p of the variables for each class are obtained, either by turesofasystemclass,isthereforedesirable. Tothisend
[ observation or by numerical simulation. Now, different
we propose the usage of simplified stochastic dynamical

observedsystems shouldbe classifiedwith respectto the
2 systems with time dependent coefficients to abstract the
reference classes.
v time evolution of system classes. Once the time depen-
1 In seismology, for example, numerical simulations of dent coefficients are obtained, they serve as a reference
0 earthquakesarepossible,butsotime consumingthatfor foraninexpensivedata-systemclasscomparison. Wede-
9 agivenobservedseismiceventamatchingsimulationcan velopboth steps, the ’model training’, i.e. the coefficient
7
rarely be found [1, 2]. This demands for robust methods determination and the ’model selection’, i.e. the classi-
0
. toidentifycharacteristicpatternswithinalimitedsample fication of an observation into the previous learned sys-
1 ofsimulatedwaveformswhichthencanbeusedtoclassify tem classes. This is done within informationfield theory
0
seismicobservationswithrespecttothelearnedpatterns. (IFT)[5],aBayesianinference frameworkforcontinuous
6
In astrophysics, for instance, the field equations of quantities.
1
: general relativity predict the existence of periodic space As periodicity plays a remarkable role in the descrip-
v
time perturbations, so called gravitationalwaves. These
i tionofdynamicalsystemsweusethestochasticoscillator
X events are generated in violent cosmic events, like the
equation
r merging of neutron stars and block holes [3], or by vari-
a oussupernovaexplosions. Unfortunately,itisoftencom- d2x(t) dx(t)
+γ(t) +ω2(t)x(t)=ξ(t) (1)
putationally not feasible to simulate these processes for
dt2 dt
all potential initial conditions. A complete catalogue
of possible gravitational wave forms is therefore out of as our baseline model. Here the time dependent fre-
reach. Identification and classification of gravitational quency ω(t) and damping factor γ(t) characterize each
wave signals need to abstract the few available signal system class and the external force ξ(t) is regarded as
templates in an appropriate fashion, to be able to rec- a random quantity which generates the variance within
ognize similar, but not exactly known waveforms in real eachclass. ω(t)andγ(t)havetobelearnedfromatrain-
data. Numerical simulations of supernova explosions ex- ingdatasetofthissystemclass. Thissystemclassisthen
identified by its ω(t) and γ(t), while we assume that the
exact realization of ξ(t) is unknown and differs from re-
alizationtorealizationwithinthatclass. Inessence,ω(t)
∗[email protected] and γ(t) serve as a fingerprint of the system while ξ(t)
2
plays the role of the random seed generating individual tionalwavedataanalysis. Therethe qualityandnumber
realization. The non-stationarity of Eq. (1), induced by of the training samples is often not sufficient to train a
the time dependent parameters γ(t) and ω(t) allows us PCA or NN. A PCA based classifier,if trained with suf-
to also model anharmonic system classes which do not ficiently many templates might be too rigid with respect
show a clear oscillating behavior. to systematic errorsinthe training dataset. PCAbased
The adopted stochastic model, here the oscillator in classifierscomparelinearcombinationsofrigidtemplates
Eq. (1), is intended to be a simplified description of a to the data, while DSC might turn out to be more flex-
much more complex system in reality. We intend to find ible as it classifies the signal covariance structure. NN
empirical parameters, here ω(t) and γ(t), which char- might outperform both, PCA and DSC. However, they
acterize the system class well, without them necessarily do not provide easy insight into what features of a sig-
being unique, perfectly optimal or even physically inter- nal triggers the classification and typically need bigger
pretable. For this reason, we will be pragmatic when training samples. In the following we will show that the
determining the parameters. The strategy presented in DSC can properly distinguish between different system
this work is first to analyze system classes in terms of classes using just a few training samples for each class.
the class’time-wise varying frequency ω(t) and damping A comparison of DSC to NN and PCA is left for future
factor γ(t). To do so, we apply a hierarchical Bayesian worksas itwill stronglydependonthe use caseathand.
ansatz within the framework of IFT as this also allows
us to simultaneously reconstruct the two point correla-
tion function of ω(t) andγ(t). For eachanalyzedsystem C. Structure of this work
class, ω(t) and γ(t) serve as abstract classification sig-
natures to which observational data can be confronted. This paper is structured as follows: In Sec. II we in-
Consequently, we introduce a Bayesian model compari- troduce some basic notation for signal inference in IFT.
son approach using these signatures to state the proba- In Sec. III we establish the training algorithm to infer
bility that a givenobservationis explained by one of the the time dependent fields ω and γ from training data.
learned system classes. In Sec. IV the inferred fields are used to construct our
model selection algorithm. After deriving the whole dy-
namicsystemclassifier(DSC)algorithmweapplyitto a
B. Previous work realistic test scenario in Sec. V. We conclude in Sec. VI.
In the last decades, extensive research has been made
II. INFERENCE OF FIELDS
withinthefieldofestimatingparametersofdeterministic
and stochastic systems, respectively [6–13]. Most atten-
tion has been drawn to local and global nonlinear opti- A. Basic notation
mizationmethods. Theuseofstochasticprocesseswithin
the global optimization approaches is called Bayesian First of all, we have to establish some notations and
global optimization or the random function approach. In basic assumptions. To do so, we primarily follow the no-
a more generic setting general parameter estimation has tationusedinIFT.Inthispaperwewillsupposethatwe
been performed extensively by maximum likelihood esti- are analyzing a discrete set of data d = (d ,...,d )T ∈
1 r
mations [14, 15]. In financial mathematics the parame- Rr, which may depend on the underlying signal s : S →
terestimationforstochasticmodelshasbeenextensively R. By S we denote the continuous subspace of a u-
studiedinafrequentistmaximumlikelihoodaswellasin dimensional Euclidian space on which the signal is de-
a Bayesian framework [16]. fined.
In the past, model selection has primarily been per- From the principal of minimal updating [27] and the
formed by likelihood ratio tests [17]. Due to the enor- principalofmaximumentropy [28]itfollowsthataGaus-
mous increase of computational power, Bayesian meth- sian is the proper probability distribution for a quantity
ods have been coming into use more and more [18]. which is only characterized by its first and second mo-
The principal component analysis (PCA), emerging mentum, as it often occurs in physical experiments. By
from [19] has successfully been used to reduce data to
be represented by a linear superposition of a few uncor- G (φ,Φ)= 1 exp −1φ†Φ−1φ (2)
related principle components [20, 21]. As a PCA might |2πΦ|1/2 2
(cid:18) (cid:19)
performpoorlyincasethe superpositionprinciple isvio-
lated, neural networks (NN) for pattern recognition and we denote a multivariate Gaussian probability distribu-
model selection are used more extensively [22]. Despite tionfunction ofa continuous field φ. |Φ| denotes the de-
the success of NN in the domain of pattern recognition terminantandφ†thetransposedandcomplexconjugated
[23] there is still an ongoing discussion how to properly φ. The covarianceΦ=hφφ†iG(φ,Φ) canbe regardedas a
adjust NN to their desired task [24–26]. We believe that function of two arguments Φ(x,y)=hφ(x)φ(y)†i or as a
the dynamic system classifier (DSC) can perform better matrix-likeoperatorΦ =hφ φ†i,whereweintroduced
xy x y
than a PCA or NN, especially in the field of gravita- the index notations for functions ψ = ψ(x). Under the
x
3
assumptions of a statistical stationary or homogeneous βt γt
process over a u-dimensional space
s
Φ =C(x−y) (3)
xy
one canshowthat the covarianceΦ becomes diagonalin
R(s)
Fourier space,
Φ =(2π)uδ(k−q)P (k). (4)
kq φ
x
t
We use k to denotefrequenciesinthe Fourierconvention
f = dteiktf and f = dkf e−ikt. In Eq. (4) we
k t t 2π k
have also introduced the power spectrum P(k) which is
R R Figure 1: Hierarchical Bayes model for the key quantities.
identical to the Fourier transform of C(x−y). In order
Thetwofieldsβ andγ formtogetherasignals,definingthe
to apply any operator to a field we have to specify the t t
response R(s), with which the system reacts to the driving
scalar product, which we take as
white noise ξ. The application of R(s) on white Gaussian
noise solves Eq. (1) and consequently yields to its solution
φ†ψ = dxφ(x)ψ(x) ∀φ,ψ :S →R. (5) x=R(s)ξ. Fromsufficienttrainingdatavectorsxaplausible
ZS signal should be inferred in thelearning phase.
Stationary Gaussian processes with zero mean are fully
determined by their power spectrum. Here we aim at
given the observed d, we use the maximum a posteriori
characterizingnon-stationaryprocesses,forwhichEq.(3)
ansatz (MAP). The outcome of the MAP- Ansatz states
and (4) do not hold. We constructa non-stationarypro-
the most probable field configuration
cess x(t) by Eq. (1), in which a fixed time evolution of
ω(t) and γ(t) imprints non-stationary correlation struc-
m=argmax{P(s|d)}. (7)
tures onto x(t) for a set of noise ξ(t) realizations. The
s
noiserealizationsξ(t)aswellasthecharacteristicsofω(t)
and γ(t) of a model class are themselves assumed to be In the following sections we will now discuss the likeli-
realizations of stationary processes. Nevertheless, x(t) is hood and the prior of the evolution of a stochastic sys-
non-stationary for fixed ω(t) and γ(t). tem, which is described by Eq. (1).
B. Signal inference III. MODEL TRAINING
Informative physical experiments provide data from A. The likelihood of a stochastic differential
which an unknown signal of interest can be inferred. equation
Since there might be infinitely many possible signalfield
configurations leading to the same data set, there is no
As outlined, we use an oscillator model with evolving
exact solution to this inference problem. Consequently,
frequencyω anddampingfactorγ tocharacterizeaclass
we have to use probabilistic data analysis methods to
of systems. To ensure strict positive definiteness of the
obtainthe mostplausiblesignalfieldincluding itsuncer- time dependent frequency we parametrize it as ω2(t) =
tainty. ω2eβ(t) in Eq. (1). Here we choose time units such that
0
Given a set of data d, the posterior probability dis- ω = 1. Consequently, the actual fluctuations of the
0
tribution P(s|d) describes how probable the signal s is frequency are characterizedby β(t).
given the observed data d. This posterior can be calcu-
Figure 1 shows a Bayesian network to infer the time
lated according to Bayes’ theorem,
dependent parameters β and γ from the solution x of
t t t
a stochastic differential equation as in Eq. (1). The in-
P(d|s)P(s)
P(s|d)= , (6) ferred vector s = (βt,γt) will serve as the characteristic
P(d)
signature of a system class in section IV. To construct
this Bayesian inference network we note that the linear-
which is the quotient of the product of the likelihood
ity of Eq. (1) implies
P(d|s)andthesignalprior P(s)dividedbytheevidence
P(d). Thelikelihooddescribeshowlikelyitistomeasure x=Rξ, (8)
the observed data set d given a signal field s. It covers
all processes that are relevant for the measurement of d.
and in scalar product notation implying
The prior characterizesall a-prioriknowledge and there-
fore does not depend on the data itself. As we are inter-
estedto findthe mostplausible signalfieldconfiguration xt = dt′Rtt′ξt′, (9)
Z
4
where we have introduced the response operator Rtt′ = αβ, qβ σβ
R(t,t′). The force ξ to the oscillator is assumed to be
t
white Gaussian noise, i.e. P(ξ|Ξ) = G (ξ,Ξ), with a
given diagonal and constant covariance Ξtt′ = δtt′Ξ˜ and
the amplitude Ξ˜. Given that the response operator R τβ
and its inverse R−1 exist, Eq. (9) can be rewritten as
dt′R−tt′1xt =ξt, (10) βt
Z
which allows us to identify the functional form of the re-
Figure 2: The spectral parameters α andq definetogether
sponseoperatorRtt′ bycomparingEq.(10)withEq.(1). withthesmoothnessenforcingparamβeterσ βthepriorforthe
Hence Rtt′ is a reformulation of the differentials1 in shapeofthespectralparametersτ (k). Henβcethecorrelation
β
Eq. (1) and is defined as structureof β is described by τ (k).
t β
−1
R(tst)’ :=δt(t2′)−γ(t)δt(t1′)+eβ(t)δtt′. (11)
(cid:16) (cid:17) Gaussian distributions,
As already used in Eq. (11) we will from now on refer to
the responseoperatorasR(s) toindicate thatitdepends P(βt|Ω)=G(βt,Ω), (13)
onthesignals,whichischaracteristicforasystemclass. P(γ |Γ)=G(γ ,Γ). (14)
t t
With Eq. (9) the likelihood P(x|s) becomes
The covariances Ω and Γ describe the strength of the
P(x|s)= Dξδ(x−R(s)ξ)G(ξ,Ξ) temporal correlations of β and γ, respectively and thus
the smoothness of their fluctuations. A convenient
Z
=G x,R(s)†ΞR(s) , (12) parametrization of the covariances can be found, if our
knowledge of the parameters does a priorinot single out
(cid:16) (cid:17) any instance; i.e. correlations only depend on time in-
where wehavemarginalizedoverallpossible realizations
tervals. This is equivalent to assume β and γ to be
t t
of the driving force ξt. We see that xt is modeled as statistically stationary. Under this assumption, Γ and
a Gaussian random field with a temporarily structured Ω are diagonal in the Fourier space representation, see
covariance X = R(s)†ΞR(s). Thus, X is not of the Eq. (4).
from given by Eq. (3) and (4), respectively and specified
To keep notations simple, we will only show the cal-
a non-stationary process, which is characterized by s. culations for Ω in full detail since the ones for Γ can be
performed analogously. Later we will see that the priors
for β and γ only differ by hyperprior-parameters.
t t
Under the assumption of statistical stationarity, we
choose the following ansatz for the covariance:
B. Prior assumptions
Ω= eτβ(k)Ω (15)
k
AsthelikelihoodinEq.(12)onlydescribeshowβ and
t k
X
γ are transformed into x we need to model our prior
t t
knowledgeaboutthe time evolutionofthese parameters. Hereτ(k)arespectralparameters,determiningthepower
To do so we briefly outline a hierarchical Bayesian prior spectrum Pβ(k) and Ωk are projection operators onto
ansatz,whichleadstothe socalledcritical filter [29,30]. spectralbands,withapproximatelyidenticalpowerspec-
Figure 2 shows this hierarchicalparametermodel, which trumvalues. ByΩ−k1werefertothepseudoinverseofthe
will be introduced in the following. We assume a priori projection operators. Since the covariances Γ and Ω are
that β as well as the γ parameters obey multivariate unknown, one has to introduce another prior for them,
t t
i.e. a prior for properly describing the spectral param-
eters τ(k), for each β and γ. As the power spectrum’s
shape might span over several orders of magnitude, this
1 ThedefiningpropertyofaDiracdeltadistributionδ(x)is requiresalogarithmicallyuniformpriorfor eachelement
of the power spectrum and a uniform prior P for each
∞ ∞ un
Z−∞f(t)δ(t−a)dt≡Z−∞f(t)δtadt=f(a). spectral parameter τk, respectively. In accordance with
[29,31]wethereforeinitiallyassumeinverse-Gammadis-
Itsderivativesaregivenby tributions for the individual elements
Z f(t)δ(n)(t)dt≡−Z ∂f∂(tt)δ(n−1)(t)dt. P (eτk|α , q )= qkαk−1 e−(αkτk+qke−τk) (16)
un k k
Γ(α −1)
k
5
and hence C. Parameter reconstruction
qαk−1
Pun(τk|αk, qk)= k Now we are seeking for the most probable parame-
Γ(α −1)
k ter configurations of β and γ given the training data
t t
×e−(αkτk+qke−τk) deτk (17) xt. Due to the complexity of the posterior, given by
(cid:12) dτk (cid:12) Eq. (21), it is virtually impossible to solve this prob-
(cid:12) (cid:12) lem analytically. Numerical searches through the pa-
where α and q denote shape and scale(cid:12)param(cid:12) eters for
k k (cid:12) (cid:12) rameter space using Markov Chain Monte Carlo meth-
the spectral hyperpriors, and Γ the Gamma function. If
ods [32, 33] are possible, as similar problems have been
α → 1 and q → 0∀k, the inverse Gamma distribu-
k k solved by advanced sampling techniques [34–37]. Given
tion becomes asymptotically flat on a logarithmic scale.
thattheemployedstochasticpartialdifferentialequation
In practice, q ≥ 0 provide lower limits for the power
k is just an empirical guess we pragmatically use here the
spectra,whichleadto amorestableinference algorithm.
numerically and computationally cheaper MAP-ansatz.
Note thatthe variationsofα andq with k canbe used
k k It is not necessary to find the best parameters, a good
to model specific spectral prior knowledge. However, in
classification will usually suffice. Rather than maximiz-
the absence of such knowledge, these will get the same
ing the posterior it is thereby convenient to define the
values α =α, q =q =const..
k k negative logarithm of the posterior P(s,τ ,τ |x) as the
β γ
information Hamiltonian
H (s,τ ,τ |x )=−logP(s,τ ,τ |x )
β γ t β γ t
Until now we have only addressed each individual ele-
ment of the power spectrum separately, but empirically =log det R(s)
wfluecktunaotwiotnhsatonmnaneiyghpboowreinrgspFeocutrraiedroscnaoltese.xhIitbiitsstthreorneg- + 1(cid:16)x† Rh(s)−1i(cid:17)†Ξ−1R(s)−1x
2 t t
fore natural for the spectrum to request some sort of
1 (cid:16) (cid:17) 1
smoothness. To enforce this behavior, we further incor- + log(det[Ω])+ β†Ω−1β
2 2 t t (22)
porate a spectral smoothness prior P [29, 30]. This
sm 1 1
spectral smoothness prior is based on the second loga- + log(det[Γ])+ γ†Γ−1γ
rithmic derivative of the spectral parameters τ. Up to a 2 2 t t
1
normalization constant Psm can be written as + (αi−1)†τ +qi†e−τi + 2τi†Tτi,
1 i∈{Xβt,γt}
Psm(τ|σ)∝exp − τ†Tτ (18) +H ,
2 0
(cid:18) (cid:19)
where we have absorbed all terms constant in β , γ , τ ,
with t t β
and τ into H .
γ 0
2
1 ∂2τ By this reformulation the MAP solution is now seek-
τ†Tτ = d(logk) k , (19)
σ2 ∂(logk)2! ing for the minimum of Eq.(22). This minimum may be
Z foundby taking the derivativeofthe informationHamil-
tonianwithrespecttoβ ,γ ,τ ,andτ ,respectivelyand
whichispunishingdeviationsfromanypower-lawbehav- t t γ β
equating them with zero. This yields four implicit equa-
ior of the power spectrum. The strength of the punish-
tions. The minimum we are seeking for may be found
ment is encoded by σ. In total, the resulting prior for
by an iterative downhill algorithm such as the steepest
the spectral parameters τ is given by the product of the
descent. To better understand the MAP solution we fo-
priors discussed above
cusontheresultingfilteringformulasofthisansatz. The
onesforthefrequencyanddampingfactorevolutionread
P(τ|α, q,σ)∝P (τ|σ) P (τ |α , q ) , (20)
sm un k k k
Yk ∂H =Tr R(s)−1 ∂R(s)
with its three given quantities αk,qk and σ. ∂βt (cid:12) " ∂βt #
(cid:12)β=βmin
(cid:12)
(cid:12)(cid:12) + 1x† ∂R(s)−1 †Ξ−1R(s)−1x
With this hierarchical Bayesian model we are able to 2 t ∂βt ! t
(23)
state the posterior: + 1x† R(s)−1 †Ξ−1∂R(s)−1x
t
P(s,τ ,τ |x ,α ,q ,σ )= 2 ∂β
β γ t ω,γ ω,γ ω,γ t
(cid:16) (cid:17)
P(x |s) +Ω−1β
t t
P(i|τ )P(τ |α , q , σ ) (21)
i i i i i
P(xt) i∈Y{β,γ} =! 0
6
and where R(s), including s = (β ,γ ), serves as an abstract
t t
operatortoclassify,identify anddistinguishbetweendif-
∂H =Tr R(s)−1 ∂R(s) ferent physicalsystems s1,s2,.... Hence each s acquires
∂γt (cid:12) " ∂γt # its R(s) from training data x(t) of its specific physical
(cid:12)(cid:12)γt=γmin system si, according to the previously described algo-
(cid:12)(cid:12) + 12x†t ∂R∂(γs)−1 †Ξ−1R(s)−1xt roifthpmos.siTbloesstiagtnealtsheSprgoibvaenbiltihtye oofbsaermveoddedlastiainwethaegaseint
t !
(24) use Bayes’ theorem
+ 1x† R(s)−1 †Ξ−1∂R(s)−1x P(d|S,s )P(s |S)
2 (cid:16) (cid:17) ∂γt t P(si|S,d)= Pi(d) i . (29)
+Γ−1γ
t
The involved likelihood turns out to be
!
=0.
P(d|s )= DxP(d|x)P(x|s )
While the filter formula for the power spectra of β is i i
t
Z
∂H =! 0→eτβ = qβ + 12 Tr βtβt†Ω−k1 k =Z Dx G(d−Robsx,N)
∂τβ (αβ −1)+ 21(cid:16)TrhΩkΩ−k1 ik(cid:17)+Tτβ ×G(x,R(s)†ΞR(s))
(25)
(cid:0) (cid:2) (cid:3)(cid:1) 1 1
∝ exp j†Dj (30)
and the one for γ is |D| 2
t (cid:18) (cid:19)
∂∂Hτγ =! 0→eτγ = (αγ −qγ1)++2112(cid:16)TTrrhγΓtγktΓ†Γ−k−k11ik(cid:17)k+Tτγ . with j =pR(s)†R†obsN−1d (31)
(26)
(cid:0) (cid:2) (cid:3)(cid:1) and
The filtering formulas in Eq. (25) and (26) have previ-
ouslybeenderivedRefs.[30,38]. Duetotheconstruction D−1 =R(s)†R†obsN−1RobsR(s)+Ξ−1. (32)
ofthe hierarchicalBayesianparametermodelthe covari-
With this equation one is able to calculate the model
ance structures of β and γ get also inferred from x .
t t t
The spectral shapes of Ω and Γ are only constrained by posterior,Eq.(29),andtostatethemostpropablemodel
s .
Eq. (20). i
Figure 3 shows an overview of the suggested hierar-
Intotal,themodellearningphaseleadstoarestriction
chical Bayesian decision algorithm. Given the hyper pa-
of possible signals to the set
rameters,the algorithmfirst learnsthe frequency β and
t
damping factor γ evolution from each training data set.
S ={s ,s ,...s }, (27) t
1 2 n Thelogarithmicpowerspectrumτ forβ aswellasforγ
t t
canbe regardedas a set of nuisance parametersthat get
which is finite in case of a finite number of considered
reconstructedfromthedatatoproperlyinferthe param-
system classes.
etersofinterestβ andγ . Afters =(β ,γ )waslearned
t t i t t
foraclassi,itservesasanabstractcharacteristicforthis
model. Withtheknowledgeofs thealgorithmisableto
i
statehowprobablethe previouslylearnedmodeli would
have caused the observed data d. This then serves as a
IV. MODEL SELECTION
proxy probability for the system classification. I
Sofarweonlyfacedtheinverseproblemtoreconstruct
timedependentparameters,suchasthefrequencyβ and V. DYNAMIC SYSTEM CLASSIFIER
t
the damping factor γ including their their power spec- ALGORITHM
t
tra from an oscillator driven by a stochastic force. From
now on we will extend our model to a measurement sce- Inferring time dependent fields, such as a time-wise
nario, involving a measurement response R and ad- varying frequency β and damping factor γ , from a
obs t t
ditive Gaussian measurement noise n ←֓ G (n,N), with stochastically driven oscillator is a non-trivial task. The
related covariance N. Consequently the data model is reliability of the dynamic system classifier (DSC) algo-
now given by rithm strongly depends on the successful and proper re-
construction of β and γ as they serve as classifiers. To
t t
d=R x+n=R R(s)ξ+n. (28) showtheprincipalcapabilitiesofoursuggestedalgorithm
obs obs
7
α , q σ α , q σ mately the complexity of this inverse problem goes back
β β β γ γ γ
to the generally highly non-linear entanglement between
the two parameters β and γ . To overcome this prob-
t t
lem we strongly advise to analyze as many realizations
τβ τγ (x=x1,x2,x3,...xl) of the same system as possible. In
section VA2 we will discuss in more detail how many
data realizations are necessary for an appropriate re-
construction of the parameters. The training part of
βt γt the DSC-algorithm is based on an iterative optimization
s scheme, where certain parts of the problem get alterna-
tivelyoptimizedinsteadofthewholeproblemsimultane-
s1 s2 s3 ously. To some degree the optimization results are sen-
sitive to the starting values due to the non-convexity of
the considered Hamiltonian. However, remaining degen-
R(s)
eraciesbetweenβ andγ after exploiting sufficient large
training t t
trainingdatasets areirrelevant,asthesedo notstrongly
discriminate between the members of the training set of
x
t a given system class.
Based on our experience with variations of the DSC-
classification algorithm we propose the following scheme:
d 1. Initialize the algorithm with naive values, such as
β =γ =0 and τ =const.∀k.
t t k
Figure 3: Complete Bayesian network and algorithm. The 2. Infer β and γ via an iterative downhill algo-
t t
solid arrows show graphical stepsby theprobabilistic depen- rithm, such as steepest descent using the informa-
dencies. The top four boxes indicate the hyper parameters
tion Hamiltonian Eq. (22), as well as its gradient
α ,q ,andσ . Belowthemfollow thespectralparam-
β/γ β/γ β/ω Eq. (23) and Eq. (24). A more sophisticated mini-
etersτ andthetupleofthefrequencyβ andthedamping
β/γ t mizationscheme,suchas non-linearconjugate gra-
factorγ ,whichformtogetherthesignals. Eachsystemclass
t
dient, is conceivable to speed up the algorithmbut
i is thereby characterized by its s . Beneath that comes the
i
response operator R(s), the training data x , and finally the it would require the full Hessian of Eq. (22). Mul-
t
observed data, according to Eq. 28. The dashed lines on the tiple test runs have shown that it is sufficient to
left hand side of the figure display the workflow of the al- evaluate a simplified Hamiltonian
gorithm to learn the time evolution of β and γ , including
t t
their power spectra from the training data set. On the right H (s|x1,x2,...xl)∝
hand side the dotted lines display the workflow of the algo- 1 −1 † −1
rithm to state which of the previously learned training data x† R(s) Ξ−1 R(s) x
set explains theobserved datad best. 2 i i
i=X1,... (cid:18)h i (cid:19) (cid:16) (cid:17)
1 1
+ log(det[Ω])+ β†Ω−1β
2 2 t t
wewillfirstdiscussthealgorithmtoinfertime-dependent
1 1
parameters of a dynamic system driven by a stochastic + log(det[Γ])+ γ†Γ−1γ ,
force as described in Sec. III. Subsequently, we use the 2 2 t t
inferredparameterstotesttheperformanceofthemodel
and its corresponding gradient. The simplified
selection algorithm as described in Sec. IV.
Hamiltonianneglectsinparticularlog det R(s)
as it appears in Eq. (22). In contrast(cid:16)to thhe diaig(cid:17)-
A. Model learning algorithm onal covariance matrixes Ω and Γ the evaluation
of the determinant of R(s) is computationally time
1. Numerical application consuming due to its complex structure. Numeri-
cal experiments with and without detR(s) did not
show a significant importance of this term.
The information Hamiltonian, Eq. (22), is a scalar
quantity defined over the configurationspace of possible
3. Use Eq. (25) and Eq. (26), respectively, to update
model parameter evolutions. In addition to the param- the priors Ω and Γ.
eters β and γ , the two spectral parameters τ and τ
t t β γ
also need to get inferred from a single system trajectory 4. Repeat step 2 and 3 until convergence. This iter-
x(t). ative scheme will take a few cycles until the algo-
Hence, the algorithm faces an underdetermined in- rithm has reached its desired convergence level.
verse problem, which is also reflected in the possibility The spaces of possible parameter configurations of β
t
of local-minima of the non-convex Hamiltonian. Ulti- and γ are huge. Consequently, it seems impossible to
t
8
(a)
8000
1.0 Orginalβ
x
1 x
1,2
6000 x
2 0.5 x
1,2,3
x (cid:3)
4000 x34 22ω/ω0 0.0 xx1,2,3,4
n 1,2,3,4,5,
m] 2000 xx56 βl(cid:2)t −0.5 x1,2,3,4,5,6
[
xt() 0 −1.0
−2000 s 0.5
ual 0.0
d
−4000 Resi−−10..05
0 5 10 15 20
−6000 time[2π/ω0]
0 5 10 15 20
(b)
time[2π/ω0]
1.0
Figure 4: Six wave realizations according to Eq. (9) for the 0.5
same model class, i.e. all x share the same γ and β . The
gaimvepnlifcyoinngstafonrtcecoξvawraiasndcreawΞ˜nt=fr5o.mP(ξ|Ξ)=Gt(ξ,Ξ),twitha 2/ω(cid:3)0 0.0 Ox1r,2ginalγ
1
(cid:2) x
γt −0.5 1,2,3
x
1,2,3,4
judge whether the algorithm has converged into the de- x1,2,3,4,5
−1.0
sired global minimum or some local minimum. It might x1,2,3,4,5,6
alsohappenthatthereconstructedfieldsdisplayfeatures
1.0
s
whichareoriginallycausedbytheexcitingforceξandnot al 0.5
u
by the frequency and damping factor itself. These prob- esid−00..50
lems can be reduced by the above demonstrated joint R
−1.0
analysis of multiple data realizations, as we see in the 0 5 10 15 20
time[2π/ω0]
following, where we discuss the numerical tests of the
optimization scheme.
Figure5: Panel(a) showstheoriginal β aswell asitsrecon-
t
structionfordifferentnumbersofdatarealizationsx ranging
t
2. Numerical tests from two(x ) up tosix (x ). Theresiduals between
1,2 1,2,3,4,5,6
theoriginalβ andtheirreconstructionsarealsoshown. Panel
t
(b)shows thesamefor γ and itsreconstructions. Oneneeds
To test the performance of the DSC algorithm we ap- t
at least three x to get a proper reconstruction of the fields.
t
plied it to simulated but realistic training data sets (see
Otherwise, the reconstruction shows too many features im-
Fig.4). Theinitialconditions,x(0)anddx(0)/dtwereset
printed bythe drivingwhite noise.
to random values, drawn from a zero centered Gaussian
distribution with variance one. This data might repre-
sent physicalsystems whose frequency anddamping fac- rameter evolution, β and γ , are shown. Due to the
t t
tor are changing overtime. For example, in astrophysics ahead mentioned degeneracy of the inverse problem we
one could expect such a behavior from a gravitational expectedthatafewwaverealizationsareneededinorder
wave burst caused by a supernova [4, 39]. to get reliable reconstructionsof the parameters. This is
In the following tests we used a regular grid with 104 clearly visible in the outcome shown in Fig. 5, as the
pixels and the signal inference library Nifty [40] to im- residuals between original and reconstructed parameters
plement the algorithm. Fig. 4 shows six realizations of are largest if one only uses two timelines. Consequently,
the same simulated system class. This means that the one needs multiple wave realizations froma system class
wave realizations, calculated according to Eq. (9) with to get a proper reconstructionand to prevent overfitting
Ξ˜ = 5, share the same β and γ . Note that the waves of the classe’s parameters β and γ . However,these two
t t t t
displayed in Fig. 4 are not just rescaled versions of the parameters are not intended to describe the exact evo-
samewavetemplate. ForthedescribedMAPreconstruc- lution of the system. They only serve as a abstract and
tion we used α = 1, q = 10−30, and σ = 2 for both descriptive fingerprint of a system class, which is con-
γ and β . In Fig. 5 the inferred parameters, β and structed on the basis of Eq. (1). In case the wave real-
t t rec
γ including their residuals, as well as the original pa- izations x(t) of a system class do not provide sufficient
rec
9
(a)
100
10−2 700
10−4
350
10−6
m]
1100−−108 xt()[ 0 x1
10−12 −350 x2
k() 10−14 OriginalP(k) −700 x3
P 10−16
x
10−18 1,2
x (cid:3)
10−20 1,2,3 2ω0 0.5
10−22 x1,2,3,4 /1(cid:2)
10−24 x1,2,3,4,5 ,γt 0.0
(cid:3)
10−26 x1,2,3,4,5,6 2/ω0 −0.5
10−2180−1 100 101 2ωln(cid:2) −1.0 γrec
|k| βt βrec
0 5 10 15 20
Figure 6: The original power spectrum of βt as well as its time[2π/ω0]
reconstruction for different numbers of data realizations x
t (b)
used in its reconstruction.
40
information, i.e. have only very small amplitudes, the
inference problem becomes more and more degenerated.
m] 0
dInuctehipsrceacsiseetlyhethreecoonrisgtirnuaclteodneβsreucsaedndinγroecurdotenstottoregperno-- xt()[ −40 x1
erate the simulated waves. Nonetheless this degeneracy x2
does not destroy the performance of the DSC-algorithm −80 x3
because what counts for the model decision algorithm is
the ability of the two parameters to represent the co- 0.7
(cid:3)
variancestructure for a model classandnot whether the 2ω0
/
parameters are as those generating the timelines of the 1(cid:2) 0.2
γt
model classes. In Figure 6 the inferred power spectrum ,
ofβ is showntogetherwith the originalone. The power 2ω(cid:3)0 −0.3
t /
spectrum at small |k|, which corresponds to large-scale 2ω
correlations,iswellreconstructed. Incontrast,thepower ln(cid:2) −0.8 γrec
spectrumatlarge|k| is underestimated, whichmay have βt −1.3 βrec
various reasons. Small short term variations in γt are 0 5 10 15 20
nearly indistinguishable from random noise variations in time[2π/ω0]
ξ. Tobetter reconstructthe small-scalecorrelationfunc-
tionmanymorerealizationsofthesystemclasswouldbe Figure 7: In panel (a) the training data set of s2 with its
needed. reconstructed parameters, ωrec and γrec is shown. Panel (b)
In summary we conclude that the presented algo- shows the same for s3.
rithm can reproduce time dependent parameters from a
stochastic differential equation. In the next section we
willusetheinferredparameters,βrec andγrec todiscrim- Fig. 4, by s2 and s3 we refer to the training data sets
inate between different models si. shown in Fig. 7. It is a trivial task to distinguish be-
tween the models s , if the means of β and γ differ
1,2,3 t t
byseveralordersofmagnitude. Toavoidsuchtrivialsit-
B. Performance of model selection algorithm uations,allβ ’sandγ ’sweredrawnfromthesamepower
t t
spectrum
To show that β and γ can indeed serve as an ab-
t t
42
stractfingerprintofclasseswetrainedouralgorithmwith P(k)= . (33)
three different training sets, generated from three differ- (1+|k|)12
ent models, specified by the parameters s ,s , and s .
1 2 3
Eachofthemhaddifferentβ andγ ,andforeachmodel In a second step we generated one noisy data set for
t t
we used three wave realizations using the s to train the each model according to Eq. (28), using new random re-
i
classificator. By s we refer to the training data set in alization for ξ and n. The underlying model for each
1
10
(a)
data set is referred to as the correct model j from now 2000
on. The additive noise n is white and Gaussian, i.e
n ←֓ G (n,σnoise). σnoise was tuned to a specific signal- 1500
to-noise ratio (SNR) which we define as
1000
σ2
x(t)
SNR= . (34)
σn2oise ]m 500
[
d
Thereby σ refers to the variance of the wave realiza-
x(t)
tions over the displayed period. For the observational 0
response R , see Eq. (28) we assume the unity opera-
obs
tor Robs(t,t′) = δtt′. Fig. 8 shows three different SNR −500
scenariosforoneourwavesgeneratedfroms . Todemon-
1
strate the performance of the model selection algorithm
−1000
we choose P (s |S)= 1 for all s , i.e. we did not prefer 0 5 10 15 20
any model. Thiis allows3us to onliy look at the likelihood time[2π/ω0]
of each model instead of its posterior. (b)
In Table I we give the differences of the log-likelihood 6000
∆ =surprise of model i−surprise of correct model j
ij 4000
=H(d|s )−H(d|s )
i j
=−logP(d|s )+logP(d|s ) (35)
i j 2000
with i,j ∈ {1,2,3}, for various SNR. Note that the in-
]m 0
formationHamiltonian, the negative log probability, can [
d
also be considered as the amount of surprise. The larger
the value ∆ is, the less plausible an assumed class i is −2000
ij
compared to the correct class j. It would be as more
surprising from the perspective of the data that the as- −4000
sumedmodeliiscorrectcomparedtothecorrectj. Upto
a SNR= 0.01 all ∆ ≥ 0 and all ∆ = 0, which means
ij ii −6000
that all datasets are correctly classified. If the SNR is 0 5 10 15 20
worse than 0.01 the algorithm starts to give misleading time[2π/ω0]
classification results, however, only on the 1σ level, and (c)
therefore not with convincing significance. As one intu-
60000
itively expects the algorithm performs better in case of
high SNR, because the absolute differences between the
40000
likelihoods are the largest within this regime.
20000
VI. CONCLUSION
]m 0
[
d
We have established the dynamic system classifier
(DSC) algorithm for model selection between dynamic −20000
systems. The algorithm consists of two steps. First, it
analyzes training data from system classes to construct
−40000
abstractclassifyinginformationforeachmodelclass. For
distinguishing oscillating systems, a natural basis is the
−60000
systems’ time dependent frequency and damping factor 0 5 10 15 20
evolution. In the second step the algorithm confronts time[2π/ω0]
data with the previously learned models and states the
probability, which of the learned models s explains the Figure 8: Panel (a), (b), and (c) show data realizations,
databest. Withthesecapabilities,DSCisapowerfultool Eq. (28) of s1. The SNR was tuned to 10, 0.1 and 0.001,
to analyze stochastic and dynamically evolving systems. respectively.
Itcanabstractasetofsampletimelinesintocharacteris-
ticcoefficientswhichencodeanon-stationarycorrelation
structure of the signals of the class.
The theoretical foundations of the first step of DSC is