Particle filter-based Gaussian process optimisation for parameter - - PowerPoint PPT Presentation

particle filter based gaussian process optimisation for
SMART_READER_LITE
LIVE PREVIEW

Particle filter-based Gaussian process optimisation for parameter - - PowerPoint PPT Presentation

Particle filter-based Gaussian process optimisation for parameter inference IFAC World Congress 2014, Cape Town, South Africa, August 28, 2014. Johan Dahlin johan.dahlin@liu.se Division of Automatic Control, Linkping University, Sweden.


slide-1
SLIDE 1

Particle filter-based Gaussian process optimisation for parameter inference

IFAC World Congress 2014, Cape Town, South Africa, August 28, 2014.

Johan Dahlin

johan.dahlin@liu.se

Division of Automatic Control, Linköping University, Sweden.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-2
SLIDE 2

This is collaborative work with

  • Dr. Fredrik Lindsten (University of Cambridge, United Kingdom)

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-3
SLIDE 3

Summary

Aim

Efficient likelihood parameter inference in nonlinear SSMs.

Methods

Gaussian process optimisation. Sequential Monte Carlo methods.

Contributions

Decreased computational cost compared with popular methods. Interesting method for solving other costly optimisation problems.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-4
SLIDE 4

Example: Modelling volatility in OMXS30 returns

900 1000 1100 1200 1300

Month Daily closing prices jan 12 mar 12 maj 12 jul 12 sep 12 nov 12 jan 13 mar 13 maj 13 jul 13 sep 13 nov 13 jan 14

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-5
SLIDE 5

Example: Modelling volatility in OMXS30 returns

−4 −2 2 4

Month Daily log−returns jan 12 mar 12 maj 12 jul 12 sep 12 nov 12 jan 13 mar 13 maj 13 jul 13 sep 13 nov 13 jan 14

−3 −2 −1 1 2 3 −4 −2 2 4

Theoretical Quantiles Sample Quantiles Daily log−returns Density

−4 −2 2 4 0.0 0.1 0.2 0.3 0.4

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-6
SLIDE 6

Example: Modelling volatility in OMXS30 returns

xt+1|xt ∼ N

  • xt+1; φxt, σ2

, yt|xt ∼ N

  • yt; 0, β2 exp(xt)
  • .

Task: Estimate θML argmax

θ∈Θ

ℓ(θ), with θ = {φ, σ, β}.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-7
SLIDE 7

Overview of the algorithm

(i) Given iterate θk, estimate the log-likelihood ℓk ≈ ℓ(θk). (ii) Given {θj, ℓj}k

j=0, create a surrogate cost function of ℓ(θ).

(iii) Select a new θk+1 using the surrogate cost function.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-8
SLIDE 8

Overview of the algorithm

(i) Given iterate θk, estimate the log-likelihood ℓk ≈ ℓ(θk).

Estimated using a particle filter.

(ii) Given {θj, ℓj}k

j=0, create a surrogate cost function of ℓ(θ).

The predictive distribution of a Gaussian process.

(iii) Select a new θk+1 using the surrogate cost function.

Acquisition function, a heuristic based on the predictive distribution.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-9
SLIDE 9

Particle filtering: overview

Resampling Propagation Weighting

Given the particle system

  • x(i)

1:T , w(i) 1:T

N

i=1,

the filtering density can be approximated by

  • pθ(❞x1:T |y1:T ) =

N

  • i=1
  • w(i)

T

N

k=1 w(k) T

  • δx(i)

1:T (❞x1:T ).

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-10
SLIDE 10

Particle filtering: log-likelihood estimator

Estimator:

  • ℓ(θ) =

T

  • t=1

log

  • N
  • i=1

w(i)

t

  • − T log N.

Statistical properties (CLT): √ N

  • ℓ(θ) −

ℓ(θ) + σ2

l

2N

  • d

− → N

  • 0, σ2

l

  • .

Error in the log−likelihood estimate Density −1.0 −0.5 0.0 0.5 1.0 0.0 0.5 1.0 1.5 2.0 −0.5 0.0 0.5 Error in the log−likelihood estimate −3 −2 −1 1 2 3 −0.5 0.0 0.5 Theoretical Quantiles Sample Quantiles

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-11
SLIDE 11

Gaussian process regression: overview

We assume a priori that ℓ(θ) ∼ GP

  • m
  • θ
  • , κ
  • θ, θ′

, which gives the posterior predictive distribution ℓ(θ⋆)|Dk ∼ N

  • µ(θ⋆|Dk), σ2(θ⋆|Dk) + σ2

l

  • , with Dk = {θj,

ℓi}k

j=1.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-12
SLIDE 12

Gaussian process regression: toy example

0.2 0.4 0.6 0.8 1.0 1.2

  • 400
  • 390
  • 380
  • 370
  • 360
  • 350
  • 340
  • 330

θ surrogate function

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-13
SLIDE 13

Gaussian process regression: toy example

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-14
SLIDE 14

Acquisition rule for selecting sampling points

Consider, the 95% upper confidence bound as the acquisition rule θk+1 = argmax

θ⋆∈Θ

  • µ(θ⋆|Dk) + 1.96
  • σ2(θ⋆|Dk)
  • ,

to determine the next iterate θk+1.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-15
SLIDE 15

Overview of the algorithm

(i) Given iterate θk, estimate the log-likelihood ℓk ≈ ℓ(θk).

Estimated using a particle filter.

(ii) Given {θj, ℓj}k

j=0, create a surrogate cost function of ℓ(θ).

The predictive distribution of a Gaussian process.

(iii) Select a new θk+1 using the surrogate cost function.

Acquisition function, a heuristic based on the predictive distribution.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-16
SLIDE 16

A toy example of the algorithm in action

0.2 0.4 0.6 0.8 1.0 1.2

  • 400
  • 380
  • 360
  • 340
  • 320

θ surrogate function

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-17
SLIDE 17

A toy example of the algorithm in action

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-18
SLIDE 18

Example: Modelling volatility in OMXS30 returns

50 100 150 0.0 0.4 0.8 1.2

iteration estimate from GPO

50 100 150 0.0 0.4 0.8 1.2

iteration estimate from SPSA

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-19
SLIDE 19

Example: Modelling volatility in OMXS30 returns

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-20
SLIDE 20

Example: Modelling volatility in OMXS30 returns

0.90 0.92 0.94 0.96 0.98 1.00 5 10 15 φ Density 0.00 0.05 0.10 0.15 0.20 0.25 5 10 15 σv Density 0.6 0.7 0.8 0.9 1 2 3 4 5 6 β Density

Estimator

  • φ
  • σ
  • β

Maximum likelihood (GPO) 0.98 0.11 0.93 Bayesian posterior mode 0.98 0.12 0.88 Bayesian posterior mean 0.97 0.14 0.93

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-21
SLIDE 21

Example: Modelling volatility in OMXS30 returns

−4 −2 2 4

Month Daily log−returns jan 12 mar 12 maj 12 jul 12 sep 12 nov 12 jan 13 mar 13 maj 13 jul 13 sep 13 nov 13 jan 14

−0.5 0.0 0.5 1.0

Month Estimated volatility jan 12 mar 12 maj 12 jul 12 sep 12 nov 12 jan 13 mar 13 maj 13 jul 13 sep 13 nov 13 jan 14

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-22
SLIDE 22

Conclusions

Methods

Particle filtering for log-likelihood estimation. CLT for the log-likelihood and Gaussian process modelling. Acquisition rules.

Contributions

Decreased computational cost compared with popular methods. Only makes use of cheap zero-order information.

Future work

Bias compensation of log-likelihood estimate. Approximate Bayesian computations. (New paper!) Input design. (New paper!)

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-23
SLIDE 23

Thank you for your attention!

Questions, comments and suggestions are most welcome.

The paper and code to replicate the results within it are found at: http://work.johandahlin.com/.

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-24
SLIDE 24

(bootstrap) Particle filtering

Resampling Propagation Weighting

  • Resampling: P(a(i)

t

= j) = w(j)

t−1 and set

x(i)

t−1 = xa(i)

t

t−1.

  • Propagation: x(i)

t

∼ Rθ

  • xt
  • x(i)

t−1

  • = fθ(xt|

x(i)

t−1).

  • Weighting: w(i)

t

= Wθ

  • x(i)

t ,

x(i)

t−1

  • = gθ(yt|xt).

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET

slide-25
SLIDE 25

Likelihood estimation using the APF

The likelihood for an SSM can be decomposed by L(θ) = pθ(y1:T ) = pθ(y1)

T

  • t=2

pθ(yt|y1:t−1), where the one-step ahead predictor can be computed by pθ(yt|y1:t−1) =

  • fθ(xt|xt−1)gθ(yt|xt)pθ(xt−1|y1:t−1) ❞xt

=

  • Wθ(xt|xt−1)Rθ(xt|xt−1)pθ(xt−1|y1:t−1) ❞xt.

pθ(yt|y1:t−1) ≈ 1 N

N

  • i=1
  • Wθ(xt|xt−1)δx(i)

t ,

x(i)

t−1 ❞xt = 1

N

N

  • i=1

w(i)

t .

AUTOMATIC CONTROL REGLERTEKNIK LINKÖPINGS UNIVERSITET