nips nips2007 nips2007-213 knowledge-graph by maker-knowledge-mining

213 nips-2007-Variational Inference for Diffusion Processes

Source: pdf

Author: Cédric Archambeau, Manfred Opper, Yuan Shen, Dan Cornford, John S. Shawe-taylor

Abstract: Diffusion processes are a family of continuous-time continuous-state stochastic processes that are in general only partially observed. The joint estimation of the forcing parameters and the system noise (volatility) in these dynamical systems is a crucial, but non-trivial task, especially when the system is nonlinear and multimodal. We propose a variational treatment of diffusion processes, which allows us to compute type II maximum likelihood estimates of the parameters by simple gradient techniques and which is computationally less demanding than most MCMC approaches. We also show how a cheap estimate of the posterior over the parameters can be constructed based on the variational free energy. 1

Reference: text

Summary: the most important sentenses genereted by tfidf model

sentIndex sentText sentNum sentScore

1 uk Abstract Diffusion processes are a family of continuous-time continuous-state stochastic processes that are in general only partially observed. [sent-16, score-0.218]

2 The joint estimation of the forcing parameters and the system noise (volatility) in these dynamical systems is a crucial, but non-trivial task, especially when the system is nonlinear and multimodal. [sent-17, score-0.541]

3 We propose a variational treatment of diffusion processes, which allows us to compute type II maximum likelihood estimates of the parameters by simple gradient techniques and which is computationally less demanding than most MCMC approaches. [sent-18, score-0.579]

4 We also show how a cheap estimate of the posterior over the parameters can be constructed based on the variational free energy. [sent-19, score-0.423]

5 1 Introduction Continuous-time diffusion processes, described by stochastic differential equations (SDEs), arise naturally in a range of applications from environmental modelling to mathematical ﬁnance [13]. [sent-20, score-0.432]

6 In this work we develop a novel variational approach to the problem of approximate inference in continuous-time diffusion processes, including a marginal likelihood (evidence) based inference technique for the forcing parameters. [sent-23, score-0.74]

7 In general, joint parameter and state inference using naive methods is complicated due to dependencies between state and system noise parameters. [sent-24, score-0.289]

8 We work in continuous time, computing distributions over sample paths1 , and discretise only in our posterior approximation, which has advantages over methods based on discretising the SDE directly [3]. [sent-25, score-0.154]

9 The approximate inference approach we describe is more computationally efﬁcient than competing Monte Carlo algorithms and could be further improved in speed by deﬁning a variety of sub-optimal approximations. [sent-26, score-0.086]

10 The approximation is also more accurate than existing Kalman smoothing methods applied to non-linear systems [4]. [sent-27, score-0.079]

11 1 A sample path is a continuous-time realisation of a stochatic process in a certain time interval. [sent-29, score-0.171]

12 Hence, a sample path is an inﬁnite dimensional object. [sent-30, score-0.102]

13 1 In Section 2 and 3, we introduce the formalism for a variational treatment of partially observed diffusion processes with measurement noise and we provide the tools to estimate the optimal variational posterior process [4]. [sent-31, score-1.146]

14 Section 4 deals with the estimation of the drift and the system noise parameter, as well as the estimation of the optimal initial conditions. [sent-32, score-0.397]

15 Finally, the approach is validated on a bi-stable nonlinear system in Section 5. [sent-33, score-0.16]

16 In this context, we also discuss how to construct an estimate of the posterior distribution over parameters based on the variational free energy. [sent-34, score-0.423]

17 2 Diffusion processes with measurement error Consider the continuous-time continuous-state stochastic process X = {Xt , t0 ≤ t ≤ tf }. [sent-35, score-0.561]

18 We assume this process is a d-dimensional diffusion process. [sent-36, score-0.301]

19 Its time evolution is described by the following SDE (to be interpreted as an Ito stochastic integral): dXt = fθ (t, Xt ) dt + Σ1/2 dWt , dWt ∼ N (0, dtI). [sent-37, score-0.198]

20 (1) The nonlinear vector function fθ deﬁnes the deterministic drift and the positive semi-deﬁnite matrix Σ ∈ Rd×d is the system noise covariance. [sent-38, score-0.411]

21 The diffusion is modelled by a d-dimensional Wiener process W = {Wt , t0 ≤ t ≤ tf } (see e. [sent-39, score-0.619]

22 However, it can be shown [13, 17, 6] that a range of state dependent stochastic forcings can be transformed into this form. [sent-45, score-0.084]

23 It is further assumed that only a small number of discrete-time latent states are observed and that the observations are subject to measurement error. [sent-46, score-0.127]

24 For n=1 n=1 n=1 simplicity, the measurement noise is modelled by a zero-mean multivariate Gaussian density,with covariance matrix R ∈ Rd×d . [sent-48, score-0.222]

25 3 Approximate inference for diffusion processes Our approximate inference scheme builds on [4] and is based on a variational inference approach (see for example [7]). [sent-49, score-0.702]

26 (3) As noted in Appendix A, this bound is ﬁnite if the approximate process is another diffusion process with a system noise covariance chosen to be identical to that of the prior process induced by (1). [sent-52, score-0.789]

27 However, since the variational distribution is restricted to have the same system noise covariance (see Appendix A) as the true posterior, the EM algorithm would leave this covariance completely unchanged in the M step and cannot be used for learning this crucial parameter. [sent-54, score-0.555]

28 Therefore, we adopt a different approach, which is based on a conjugate gradient method. [sent-55, score-0.102]

29 1 Optimal approximate posterior process We consider an approximate time-varying linear process with the same diffusion term, that is the same system noise covariance: dXt = g(t, Xt ) dt + Σ1/2 dWt , dWt ∼ N (0, dtI), (4) where g(t, x) = −A(t)x+b(t), with A(t) ∈ Rd×d and b(t) ∈ Rd . [sent-57, score-0.905]

30 In other words, the approximate posterior process q(X|Σ) is restricted to be a Gaussian process [4]. [sent-58, score-0.288]

31 The Gaussian marginal at time t is deﬁned as follows: q(Xt |Σ) = N (Xt |m(t), S(t)), 2 t0 ≤ t ≤ tf , (5) where m(t) ∈ Rd and S(t) ∈ Rd×d are respectively the marginal mean and the marginal covariance at time t. [sent-59, score-0.557]

32 In the rest of the paper, we denote q(Xt |Σ) by the shorthand notation qt . [sent-60, score-0.201]

33 For ﬁxed parameters θ and assuming that there is no observation at the initial time t0 , the optimal approximate posterior process q(X|Σ) is the one minimizing the variational free energy, which is given by (see Appendix A) tf FΣ (q, θ) = tf Esde (t) dt + t0 δ(t − tn ) dt + KL [q0 p0 ] . [sent-61, score-1.468]

34 The energy functions Esde (t) and Eobs (t) are deﬁned as follows: 1 Esde (t) = (fθ (t, Xt ) − g(t, Xt )) Σ−1 (fθ (t, Xt ) − g(t, Xt )) q , (7) t 2 1 d 1 Eobs (t) = (Yt − Xt ) R−1 (Yt − Xt ) qt + ln 2π + ln |R|. [sent-63, score-0.323]

35 (8) 2 2 2 where {Yt , t0 ≤ t ≤ tf } is the underlying continuous-time observable process. [sent-64, score-0.318]

36 2 Smoothing algorithm The variational parameters to optimise in order to ﬁnd the optimal Gaussian process approximation are A(t), b(t), m(t) and S(t). [sent-66, score-0.422]

37 For a linear SDE with additive system noise, it can be shown that the time evolution of the means and the covariances are described by a set of ordinary differential equations [13, 4]: ˙ m(t) = −A(t)m(t) + b(t), (9) ˙ S(t) = −A(t)S(t) − S(t)A (t) + Σ, (10) where ˙ denotes the time derivtive. [sent-67, score-0.343]

38 These equations provide us with consistency constraints for the marginal means and the marginal covariances along sample paths. [sent-68, score-0.266]

39 To enforce these constraints we formulate the Lagrangian tf Lθ,Σ = FΣ (q, θ) − ˙ λ (t) m(t) + A(t)m(t) − b(t) dt t0 tf − ˙ tr Ψ(t) S(t) + 2A(t)S(t) − Σ t0 d×d where λ(t) ∈ Rd and Ψ(t) ∈ R dt, (11) are time dependent Lagrange multipliers, with Ψ(t) symmetric. [sent-69, score-0.809]

40 First, taking the functional derivatives of Lθ,Σ with respect to A(t) and b(t) results in the following gradient functions: A Lθ,Σ (t) The gradients = A Esde (t) − λ(t)m (t) − 2Ψ(t)S(t), b Lθ,Σ (t) = b Esde (t) + λ(t). [sent-70, score-0.178]

41 (15) The optimal variational functions can be computed by means of a gradient descent technique, such as the conjugate gradient (see e. [sent-73, score-0.412]

42 Since m(t), S(t), λ(t) and Ψ(t) are dependent on these parameters, one needs also to take the corresponding implicit derivatives into account. [sent-77, score-0.106]

43 However, these implicit gradients vanish if the consistency constraints for the means (9) and the covariances (10), as well as the ones for the Lagrange multipliers (14-15), are satisﬁed. [sent-78, score-0.262]

44 One way to achieve this is to perform a forward propagation of the means and the covariances, followed by a backward propagation of the Lagrange multipliers, and then to take a gradient step. [sent-79, score-0.13]

45 The resulting algorithm for computing the optimal posterior q(X|Σ) over sample paths is detailed in Algorithm 1. [sent-80, score-0.29]

46 The estimation of the parameters related to the observable process are not discussed in this work, although it is a straightforward extension. [sent-83, score-0.104]

47 The smoothing algorithm described in the previous section computes the optimal posterior process by providing us with the stationary solution functions A(t) and b(t). [sent-84, score-0.251]

48 Therefore, when subsequently optimising the parameters we only need to compute their explicit derivatives; the implicit ones vanish since A Lθ,Σ = 0 and b Lθ,Σ = 0. [sent-85, score-0.098]

49 This leads to tf Lθ,Σ = FΣ (q, θ) − ˙ λ (t) A(t)m(t) − b(t) − λ (t)m(t) dt − λf mf + λ0 m0 t0 tf − ˙ tr Ψ(t) 2A(t)S(t) − Σ − Ψ(t)S(t) dt − tr {Ψf Sf } + tr {Ψ0 S0 } , (16) t0 At the ﬁnal time tf , there are no consistency constraints, that is λf and Ψf are both equal to zero. [sent-87, score-1.28]

50 1 Initial state The initial variational posterior q(x0 ) is chosen equal to N (x0 |m0 , S0 ) to ensure that the approximate process is a Gaussian one. [sent-89, score-0.436]

51 2 Drift The gradients for the drift function parameters θ f only depend on the total energy associated to the SDE. [sent-93, score-0.28]

52 Their general expression is given by tf θ f Lθ,Σ = θ f Esde (t) t0 4 dt, (18) where θ f Esde (t) = (fθ (t, Xt ) − g(t, Xt )) Σ−1 θ f fθ (t, Xt ) qt . [sent-94, score-0.491]

53 Note that the observations do play a role in this gradient as they enter through g(t, Xt ) and the expectation w. [sent-95, score-0.102]

54 3 System noise Estimating the system noise covariance (or volatility) is essential as the system noise, together with the drift function, determines the dynamics. [sent-100, score-0.657]

55 For example in a Bayesian MCMC approach, which alternates between sampling paths and parameters, the latent paths imputed between observations must have a system noise parameter which is arbitrarily close to its previous value in order to be accepted by a Metropolis sampler. [sent-102, score-0.529]

56 However, in our method, we can simply compute approximations to the marginal likelihood and its gradient directly. [sent-105, score-0.158]

57 In the next section, we will compare our results to a direct MCMC estimate of the marginal likelihood which is a time consuming method. [sent-106, score-0.093]

58 The gradient of (16) with respect to Σ is given by tf Σ Lθ,Σ = tf Σ Esde (t) dt + t0 where 5 Σ Esde (t) Ψ(t) dt, (19) t0 1 = − 2 Σ−1 (fθ (t, Xt ) − g(t, Xt )) (fθ (t, Xt ) − g(t, Xt )) qt Σ−1 . [sent-107, score-0.986]

59 Experimental validation on a bi-stable system In order to validate the approach, we consider the 1 dimensional double-well system: fθ (t, x) = 4x θ − x2 , θ > 0, (20) where fθ (t, x) is the drift function. [sent-108, score-0.249]

60 This dynamical system is highly nonlinear and its stationary distribution is multi-modal. [sent-109, score-0.207]

61 The system is driven by the system noise, which makes it occasionally ﬂip from one well to the other. [sent-111, score-0.236]

62 In the experiments, we set the drift parameter θ to 1, the system noise standard deviation σ to 0. [sent-112, score-0.369]

63 The time step for the variational approximation is set to ∆t = 0. [sent-115, score-0.257]

64 Figure 1(a) compares the variational solution to the outcomes of a hybrid MCMC simulation of the posterior process using the true parameter values. [sent-119, score-0.464]

65 At each step of the sampling process, an entire sample path is generated. [sent-121, score-0.102]

66 In order to keep the acceptance of new paths sufﬁciently high, the basic MCMC algorithm is combined with ideas from Molecular Dynamics, such that the MCMC sampler moves towards regions of high probability in the state space. [sent-122, score-0.108]

67 In particular, over 100, 000 sample paths were necessary to reach convergence in the case of the double-well system. [sent-124, score-0.147]

68 One can observe that the variational solution underestimates the uncertainty (smaller error bars). [sent-126, score-0.254]

69 Convergence of the smoothing algorithm was achieved in approximately 180 conjugate gradient steps, each one involving a forward and backward sweep. [sent-128, score-0.206]

70 The optimal parameters and the optimal initial conditions for the variational solution are given by ˆ σ = 0. [sent-129, score-0.308]

71 ˆ ˆ ˆ (21) Convergence of the outer optimization loop is typically reached after less then 10 conjugate gradient steps. [sent-134, score-0.102]

72 While the estimated value for the drift parameter is within 15% percent from its true value, 5 2 1. [sent-135, score-0.131]

73 The curves denote the mean paths and the shaded regions are the twostandard deviation noise tubes. [sent-159, score-0.228]

74 (b) Posterior of the system noise variance (diffusion term). [sent-160, score-0.238]

75 The plain curve and the dashed curve are respectively the approximations of the posterior shape based on the variational free energy and MCMC. [sent-161, score-0.44]

76 Furthermore, we have chosen a sample path which contains a transition between the two wells within a small time interval and is thus highly untypical with respect to the prior distribution. [sent-164, score-0.235]

77 This fact was experimentally assessed by estimating the parameters on a sample path without transition, in a time window of the same size. [sent-165, score-0.137]

78 ˆ Posterior distribution over the parameters Interestingly, minimizing the free energy Fσ2 for different values of σ provides us with much more information than a single point estimate for the parameters [14]. [sent-170, score-0.178]

79 Using a suitable prior over p(σ), we can approximate the posterior over the system noise variance via p(σ 2 |Y ) ∝ e−Fσ2 p(σ 2 ), (22) where we take e−Fσ2 (at its minimum) as an approximation to the marginal likelihood of the observations p(Y |σ 2 ). [sent-171, score-0.585]

80 A comparison with preliminary MCMC estimates for p(Y |σ 2 ) for θ = 1 and a set of system noise variances indicates that the shape of our approximation is a reasonable indicator of the shape of the posterior. [sent-173, score-0.278]

81 6 Conclusion We have presented a variational approach to the approximate inference of stochastic differential equations from a ﬁnite set of noisy observations. [sent-175, score-0.461]

82 So far, we have tested the method on a one dimensional bi-stable system only. [sent-176, score-0.118]

83 Comparison with a Monte Carlo approach suggests that our method can reproduce the posterior mean fairly well but underestimates the variance in the region of the transition. [sent-177, score-0.152]

84 Although our approach is based on a Gaussian approximation of the posterior process, we expect that one can improve on it and obtain non-Gaussian predictions at least for various marginal posterior distributions, including that of the latent variable Xt at a ﬁxed time t. [sent-180, score-0.371]

85 This should be possible by generalising our method for the computation of a non-Gaussian shaped probability density for the system noise parameter using the free energy. [sent-181, score-0.294]

86 6 We hope that the possibility of using simpler suboptimal parametrisations of the approximating Gaussian process will allow us to obtain a tractable inference method that scales well to higher dimensions. [sent-183, score-0.12]

87 A The Kullback-Leibler divergence interpreted as a path integral In this section, we show that the Kullback-Leibler divergence between the posterior process p(Xt |Y, θ, Σ) and its approximation q(X|Σ) can be interpreted as a path integral over time. [sent-185, score-0.428]

88 Consider the Euler-Muryama discrete approximation (see for example [13]) of the SDE (1) and its linear approximation (4): ∆xk = fk ∆t + Σ1/2 ∆wk , (23) 1/2 ∆xk = gk ∆t + Σ ∆wk , (24) where ∆xk ≡ xk+1 − xk and wk ∼ N (0, ∆tI). [sent-189, score-0.472]

89 The vectors fk and gk are shorthand notations for fθ (tk , xk ) and g(tk , xk ). [sent-190, score-0.586]

90 Hence, the joint distributions of discrete sample paths {xk }k≥0 for the true process and its approximation follow from the Markov property: p(x0 , . [sent-191, score-0.256]

91 , xK |Σ) = p(x0 ) N (xk+1 |xk + fk ∆t, Σ∆t), (25) N (xk+1 |xk + gk ∆t, Σ∆t), (26) k>0 q(x0 , . [sent-194, score-0.144]

92 Note thate we do not restrict the variational posterior to factorise over the latent states. [sent-198, score-0.37]

93 The distribution qt = q(Xt |Σ) is the marginal at time t for a given system noise covariance Σ. [sent-202, score-0.524]

94 It is important to realise that the KL between the induced prior process and its approximation is ﬁnite because the system noise covariances are chosen to be identical. [sent-203, score-0.442]

95 Parameter estimation in an intermediate complexity earth system model using an ensemble Kalman ﬁlter. [sent-229, score-0.118]

96 Expectation correction for smoothed inference in switching linear dynamical systems. [sent-247, score-0.098]

97 Exact and computationally efﬁcient likelihood-based estimation for discretely observed diffusion processes (with discussion). [sent-254, score-0.297]

98 A mean ﬁeld approximation in data assimilation for nonlinear dynamics. [sent-281, score-0.131]

99 Bayesian inference for nonlinear multivariate diffusion models observed with error. [sent-287, score-0.325]

100 On inference for partially observed non-linear diffusion models using the Metropolis-Hastings algorithm. [sent-321, score-0.314]

similar papers computed by tfidf model

tfidf for this paper:

wordName wordTfidf (topN-words)

[('esde', 0.469), ('tf', 0.318), ('xt', 0.247), ('diffusion', 0.232), ('variational', 0.217), ('xk', 0.207), ('qt', 0.173), ('eobs', 0.173), ('mcmc', 0.153), ('sde', 0.137), ('drift', 0.131), ('noise', 0.12), ('system', 0.118), ('posterior', 0.115), ('dt', 0.112), ('ak', 0.109), ('paths', 0.108), ('kl', 0.107), ('tk', 0.106), ('dwt', 0.099), ('fk', 0.073), ('gk', 0.071), ('multipliers', 0.069), ('process', 0.069), ('differential', 0.068), ('covariances', 0.068), ('processes', 0.065), ('gradient', 0.065), ('marginal', 0.063), ('hybrid', 0.063), ('path', 0.063), ('gradients', 0.062), ('lagrange', 0.061), ('forcing', 0.061), ('stochastic', 0.057), ('mk', 0.056), ('free', 0.056), ('kalman', 0.055), ('bk', 0.055), ('rd', 0.054), ('sk', 0.053), ('tn', 0.053), ('measurement', 0.052), ('carlo', 0.052), ('monte', 0.052), ('energy', 0.052), ('inference', 0.051), ('derivatives', 0.051), ('covariance', 0.05), ('assimilation', 0.049), ('dxt', 0.049), ('eyink', 0.049), ('ln', 0.049), ('dynamical', 0.047), ('wells', 0.043), ('dti', 0.043), ('volatility', 0.043), ('environmental', 0.042), ('nonlinear', 0.042), ('wk', 0.041), ('approximation', 0.04), ('roberts', 0.039), ('cornford', 0.039), ('archambeau', 0.039), ('aston', 0.039), ('physica', 0.039), ('integral', 0.039), ('smoothing', 0.039), ('appendix', 0.039), ('sample', 0.039), ('latent', 0.038), ('observations', 0.037), ('conjugate', 0.037), ('backward', 0.037), ('underestimates', 0.037), ('approximate', 0.035), ('parameters', 0.035), ('transition', 0.035), ('minimising', 0.035), ('vanish', 0.035), ('opper', 0.035), ('tr', 0.034), ('optimise', 0.033), ('em', 0.033), ('equations', 0.033), ('partially', 0.031), ('likelihood', 0.03), ('evolution', 0.029), ('shorthand', 0.028), ('optimal', 0.028), ('implicit', 0.028), ('interval', 0.028), ('forward', 0.028), ('gt', 0.027), ('particle', 0.027), ('ordinary', 0.027), ('dependent', 0.027), ('prior', 0.027), ('ft', 0.027), ('yt', 0.026)]

similar papers list:

simIndex simValue paperId paperTitle

same-paper 1 1.0 213 nips-2007-Variational Inference for Diffusion Processes

Author: Cédric Archambeau, Manfred Opper, Yuan Shen, Dan Cornford, John S. Shawe-taylor

2 0.20031244 17 nips-2007-A neural network implementing optimal state estimation based on dynamic spike train decoding

Author: Omer Bobrowski, Ron Meir, Shy Shoham, Yonina Eldar

Abstract: It is becoming increasingly evident that organisms acting in uncertain dynamical environments often employ exact or approximate Bayesian statistical calculations in order to continuously estimate the environmental state, integrate information from multiple sensory modalities, form predictions and choose actions. What is less clear is how these putative computations are implemented by cortical neural networks. An additional level of complexity is introduced because these networks observe the world through spike trains received from primary sensory afferents, rather than directly. A recent line of research has described mechanisms by which such computations can be implemented using a network of neurons whose activity directly represents a probability distribution across the possible “world states”. Much of this work, however, uses various approximations, which severely restrict the domain of applicability of these implementations. Here we make use of rigorous mathematical results from the theory of continuous time point process ﬁltering, and show how optimal real-time state estimation and prediction may be implemented in a general setting using linear neural networks. We demonstrate the applicability of the approach with several examples, and relate the required network properties to the statistical nature of the environment, thereby quantifying the compatibility of a given network with its environment. 1

3 0.18589564 214 nips-2007-Variational inference for Markov jump processes

Author: Manfred Opper, Guido Sanguinetti

Abstract: Markov jump processes play an important role in a large number of application domains. However, realistic systems are analytically intractable and they have traditionally been analysed using simulation based techniques, which do not provide a framework for statistical inference. We propose a mean ﬁeld approximation to perform posterior inference and parameter estimation. The approximation allows a practical solution to the inference problem, while still retaining a good degree of accuracy. We illustrate our approach on two biologically motivated systems.

4 0.18150143 146 nips-2007-On higher-order perceptron algorithms

Author: Claudio Gentile, Fabio Vitale, Cristian Brotto

Abstract: A new algorithm for on-line learning linear-threshold functions is proposed which efﬁciently combines second-order statistics about the data with the ”logarithmic behavior” of multiplicative/dual-norm algorithms. An initial theoretical analysis is provided suggesting that our algorithm might be viewed as a standard Perceptron algorithm operating on a transformed sequence of examples with improved margin properties. We also report on experiments carried out on datasets from diverse domains, with the goal of comparing to known Perceptron algorithms (ﬁrst-order, second-order, additive, multiplicative). Our learning procedure seems to generalize quite well, and converges faster than the corresponding multiplicative baseline algorithms. 1 Introduction and preliminaries The problem of on-line learning linear-threshold functions from labeled data is one which have spurred a substantial amount of research in Machine Learning. The relevance of this task from both the theoretical and the practical point of view is widely recognized: On the one hand, linear functions combine ﬂexiblity with analytical and computational tractability, on the other hand, online algorithms provide efﬁcient methods for processing massive amounts of data. Moreover, the widespread use of kernel methods in Machine Learning (e.g., [24]) have greatly improved the scope of this learning technology, thereby increasing even further the general attention towards the speciﬁc task of incremental learning (generalized) linear functions. Many models/algorithms have been proposed in the literature (stochastic, adversarial, noisy, etc.) : Any list of references would not do justice of the existing work on this subject. In this paper, we are interested in the problem of online learning linear-threshold functions from adversarially generated examples. We introduce a new family of algorithms, collectively called the Higher-order Perceptron algorithm (where ”higher” means here ”higher than one”, i.e., ”higher than ﬁrst-order” descent algorithms such as gradientdescent or standard Perceptron-like algorithms”). Contrary to other higher-order algorithms, such as the ridge regression-like algorithms considered in, e.g., [4, 7], Higher-order Perceptron has the ability to put together in a principled and ﬂexible manner second-order statistics about the data with the ”logarithmic behavior” of multiplicative/dual-norm algorithms (e.g., [18, 19, 6, 13, 15, 20]). Our algorithm exploits a simpliﬁed form of the inverse data matrix, lending itself to be easily combined with the dual norms machinery introduced by [13] (see also [12, 23]). As we will see, this has also computational advantages, allowing us to formulate an efﬁcient (subquadratic) implementation. Our contribution is twofold. First, we provide an initial theoretical analysis suggesting that our algorithm might be seen as a standard Perceptron algorithm [21] operating on a transformed sequence of examples with improved margin properties. The same analysis also suggests a simple (but principled) way of switching on the ﬂy between higher-order and ﬁrst-order updates. This is ∗ The authors gratefully acknowledge partial support by the PASCAL Network of Excellence under EC grant n. 506778. This publication only reﬂects the authors’ views. especially convenient when we deal with kernel functions, a major concern being the sparsity of the computed solution. The second contribution of this paper is an experimental investigation of our algorithm on artiﬁcial and real-world datasets from various domains: We compared Higher-order Perceptron to baseline Perceptron algorithms, like the Second-order Perceptron algorithm deﬁned in [7] and the standard (p-norm) Perceptron algorithm, as in [13, 12]. We found in our experiments that Higher-order Perceptron generalizes quite well. Among our experimental ﬁndings are the following: 1) Higher-order Perceptron is always outperforming the corresponding multiplicative (p-norm) baseline (thus the stored data matrix is always beneﬁcial in terms of convergence speed); 2) When dealing with Euclidean norms (p = 2), the comparison to Second-order Perceptron is less clear and depends on the speciﬁc task at hand. Learning protocol and notation. Our algorithm works in the well-known mistake bound model of on-line learning, as introduced in [18, 2], and further investigated by many authors (e.g., [19, 6, 13, 15, 7, 20, 23] and references therein). Prediction proceeds in a sequence of trials. In each trial t = 1, 2, . . . the prediction algorithm is given an instance vector in Rn (for simplicity, all vectors are normalized, i.e., ||xt || = 1, where || · || is the Euclidean norm unless otherwise speciﬁed), and then guesses the binary label yt ∈ {−1, 1} associated with xt . We denote the algorithm’s prediction by yt ∈ {−1, 1}. Then the true label yt is disclosed. In the case when yt = yt we say that the algorithm has made a prediction mistake. We call an example a pair (xt , yt ), and a sequence of examples S any sequence S = (x1 , y1 ), (x2 , y2 ), . . . , (xT , yT ). In this paper, we are competing against the class of linear-threshold predictors, parametrized by normal vectors u ∈ {v ∈ Rn : ||v|| = 1}. In this case, a common way of measuring the (relative) prediction performance of an algorithm A is to compare the total number of mistakes of A on S to some measure of the linear separability of S. One such measure (e.g., [24]) is the cumulative hinge-loss (or soft-margin) Dγ (u; S) of S w.r.t. a T linear classiﬁer u at a given margin value γ > 0: Dγ (u; S) = t=1 max{0, γ − yt u xt } (observe that Dγ (u; S) vanishes if and only if u separates S with margin at least γ. A mistake-driven algorithm A is one which updates its internal state only upon mistakes. One can therefore associate with the run of A on S a subsequence M = M(S, A) ⊆ {1, . . . , T } of mistaken trials. Now, the standard analysis of these algorithms allows us to restrict the behavior of the comparison class to mistaken trials only and, as a consequence, to reﬁne Dγ (u; S) so as to include only trials in M: Dγ (u; S) = t∈M max{0, γ − yt u xt }. This gives bounds on A’s performance relative to the best u over a sequence of examples produced (or, actually, selected) by A during its on-line functioning. Our analysis in Section 3 goes one step further: the number of mistakes of A on S is contrasted to the cumulative hinge loss of the best u on a transformed ˜ sequence S = ((˜ i1 , yi1 ), (˜ i2 , yi2 ), . . . , (˜ im , yim )), where each instance xik gets transformed x x x ˜ into xik through a mapping depending only on the past behavior of the algorithm (i.e., only on ˜ examples up to trial t = ik−1 ). As we will see in Section 3, this new sequence S tends to be ”more separable” than the original sequence, in the sense that if S is linearly separable with some margin, ˜ then the transformed sequence S is likely to be separable with a larger margin. 2 The Higher-order Perceptron algorithm The algorithm (described in Figure 1) takes as input a sequence of nonnegative parameters ρ1 , ρ2 , ..., and maintains a product matrix Bk (initialized to the identity matrix I) and a sum vector v k (initialized to 0). Both of them are indexed by k, a counter storing the current number of mistakes (plus one). Upon receiving the t-th normalized instance vector xt ∈ Rn , the algorithm computes its binary prediction value yt as the sign of the inner product between vector Bk−1 v k−1 and vector Bk−1 xt . If yt = yt then matrix Bk−1 is updates multiplicatively as Bk = Bk−1 (I − ρk xt xt ) while vector v k−1 is updated additively through the standard Perceptron rule v k = v k−1 + yt xt . The new matrix Bk and the new vector v k will be used in the next trial. If yt = yt no update is performed (hence the algorithm is mistake driven). Observe that ρk = 0 for any k makes this algorithm degenerate into the standard Perceptron algorithm [21]. Moreover, one can easily see that, in order to let this algorithm exploit the information collected in the matrix B (and let the algorithm’s ∞ behavior be substantially different from Perceptron’s) we need to ensure k=1 ρk = ∞. In the sequel, our standard choice will be ρk = c/k, with c ∈ (0, 1). See Sections 3 and 4. Implementing Higher-Order Perceptron can be done in many ways. Below, we quickly describe three of them, each one having its own merits. 1) Primal version. We store and update an n×n matrix Ak = Bk Bk and an n-dimensional column Parameters: ρ1 , ρ2 , ... ∈ [0, 1). Initialization: B0 = I; v 0 = 0; k = 1. Repeat for t = 1, 2, . . . , T : 1. Get instance xt ∈ Rn , ||xt || = 1; 2. Predict yt = SGN(wk−1 xt ) ∈ {−1, +1}, where wk−1 = Bk−1 Bk−1 v k−1 ; 3. Get label yt ∈ {−1, +1}; v k = v k−1 + yt xt 4. if yt = yt then: Bk k = Bk−1 (I − ρk xt xt ) ← k + 1. Figure 1: The Higher-order Perceptron algorithm (for p = 2). vector v k . Matrix Ak is updated as Ak = Ak−1 − ρAk−1 xx − ρxx Ak−1 + ρ2 (x Ak−1 x)xx , taking O(n2 ) operations, while v k is updated as in Figure 1. Computing the algorithm’s margin v Ax can then be carried out in time quadratic in the dimension n of the input space. 2) Dual version. This implementation allows us the use of kernel functions (e.g., [24]). Let us denote by Xk the n × k matrix whose columns are the n-dimensional instance vectors x1 , ..., xk where a mistake occurred so far, and y k be the k-dimensional column vector of the corresponding (k) labels. We store and update the k × k matrix Dk = [di,j ]k i,j=1 , the k × k diagonal matrix Hk = DIAG {hk }, (k) (k) hk = (h1 , ..., hk ) = Xk Xk y k , and the k-dimensional column vector g k = y k + Dk Hk 1k , being 1k a vector of k ones. If we interpret the primal matrix Ak above as Ak = (k) k I + i,j=1 di,j xi xj , it is not hard to show that the margin value wk−1 x is equal to g k−1 Xk−1 x, and can be computed through O(k) extra inner products. Now, on the k-th mistake, vector g can be updated with O(k 2 ) extra inner products by updating D and H in the following way. We let D0 and H0 be empty matrices. Then, given Dk−1 and Hk−1 = DIAG{hk−1 }, we have1 Dk = Dk−1 −ρk bk (k) , where bk = Dk−1 Xk−1 xk , and dk,k = ρ2 xk Xk−1 bk − 2ρk + ρ2 . On (k) k k −ρk bk dk,k the other hand, Hk = DIAG {hk−1 (k) (k) + yk Xk−1 xk , hk }, with hk = y k−1 Xk−1 xk + yk . Observe that on trials when ρk = 0 matrix Dk−1 is padded with a zero row and a zero column. (k) k This amounts to say that matrix Ak = I + i,j=1 di,j xi xj , is not updated, i.e., Ak = Ak−1 . A closer look at the above update mechanism allows us to conclude that the overall extra inner products needed to compute g k is actually quadratic only in the number of past mistaken trials having ρk > 0. This turns out to be especially important when using a sparse version of our algorithm which, on a mistaken trial, decides whether to update both B and v or just v (see Section 4). 3) Implicit primal version and the dual norms algorithm. This is based on the simple observation that for any vector z we can compute Bk z by unwrapping Bk as in Bk z = Bk−1 (I − ρxx )z = Bk−1 z , where vector z = (z − ρx x z) can be calculated in time O(n). Thus computing the margin v Bk−1 Bk−1 x actually takes O(nk). Maintaining this implicit representation for the product matrix B can be convenient when an efﬁcient dual version is likely to be unavailable, as is the case for the multiplicative (or, more generally, dual norms) extension of our algorithm. We recall that a multiplicative algorithm is useful when learning sparse target hyperplanes (e.g., [18, 15, 3, 12, 11, 20]). We obtain a dual norms algorithm by introducing a norm parameter p ≥ 2, and the associated gradient mapping2 g : θ ∈ Rn → θ ||θ||2 / 2 ∈ Rn . Then, in Figure 1, we p normalize instance vectors xt w.r.t. the p-norm, we deﬁne wk−1 = Bk−1 g(Bk−1 v k−1 ), and generalize the matrix update as Bk = Bk−1 (I − ρk xt g(xt ) ). As we will see, the resulting algorithm combines the multiplicative behavior of the p-norm algorithms with the ”second-order” information contained in the matrix Bk . One can easily see that the above-mentioned argument for computing the margin g(Bk−1 v k−1 ) Bk−1 x in time O(nk) still holds. 1 Observe that, by construction, Dk is a symmetric matrix. This mapping has also been used in [12, 11]. Recall that setting p = O(log n) yields an algorithm similar to Winnow [18]. Also, notice that p = 2 yields g = identity. 2 3 Analysis We express the performance of the Higher-order Perceptron algorithm in terms of the hinge-loss behavior of the best linear classiﬁer over the transformed sequence ˜ S = (B0 xt(1) , yt(1) ), (B1 xt(2) , yt(2) ), (B2 xt(3) , yt(3) ), . . . , (1) being t(k) the trial where the k-th mistake occurs, and Bk the k-th matrix produced by the algorithm. Observe that each feature vector xt(k) gets transformed by a matrix Bk depending on past examples ˜ only. This is relevant to the argument that S tends to have a larger margin than the original sequence (see the discussion at the end of this section). This neat ”on-line structure” does not seem to be shared by other competing higher-order algorithms, such as the ”ridge regression-like” algorithms considered, e.g., in [25, 4, 7, 23]. For the sake of simplicity, we state the theorem below only in the case p = 2. A more general statement holds when p ≥ 2. Theorem 1 Let the Higher-order Perceptron algorithm in Figure 1 be run on a sequence of examples S = (x1 , y1 ), (x2 , y2 ), . . . , (xT , yT ). Let the sequence of parameters ρk satisfy 0 ≤ ρk ≤ 1−c , where xt is the k-th mistaken instance vector, and c ∈ (0, 1]. Then the total number m 1+|v k−1 xt | of mistakes satisﬁes3 ˜ ˜ Dγ (u; Sc )) α2 α Dγ (u; Sc )) α2 m≤α + 2+ α + 2, (2) γ 2γ γ γ 4γ holding for any γ > 0 and any unit norm vector u ∈ Rn , where α = α(c) = (2 − c)/c. Proof. The analysis deliberately mimics the standard Perceptron convergence analysis [21]. We ﬁx an arbitrary sequence S = (x1 , y1 ), (x2 , y2 ), . . . , (xT , yT ) and let M ⊆ {1, 2, . . . , T } be the set of trials where the algorithm in Figure 1 made a mistake. Let t = t(k) be the trial where the k-th mistake occurred. We study the evolution of ||Bk v k ||2 over mistaken trials. Notice that the matrix Bk Bk is positive semideﬁnite for any k. We can write ||Bk v k ||2 = ||Bk−1 (I − ρk xt xt ) (v k−1 + yt xt ) ||2 (from the update rule v k = v k−1 + yt xt and Bk = Bk−1 (I − ρk xt xt ) ) = ||Bk−1 v k−1 + yt (1 − ρk yt v k−1 xt − ρk )Bk−1 xt ||2 2 = ||Bk−1 v k−1 || + 2 yt rk v k−1 Bk−1 Bk−1 xt + (using ||xt || = 1) 2 rk ||Bk−1 xt ||2 , where we set for brevity rk = 1 − ρk yt v k−1 xt − ρk . We proceed by upper and lower bounding the above chain of equalities. To this end, we need to ensure rk ≥ 0. Observe that yt v k−1 xt ≥ 0 implies rk ≥ 0 if and only if ρk ≤ 1/(1 + yt v k−1 xt ). On the other hand, if yt v k−1 xt < 0 then, in order for rk to be nonnegative, it sufﬁces to pick ρk ≤ 1. In both cases ρk ≤ (1 − c)/(1 + |v k−1 xt |) implies 2 rk ≥ c > 0, and also rk ≤ (1+ρk |v k−1 xt |−ρk )2 ≤ (2−c)2 . Now, using yt v k−1 Bk−1 Bk−1 xt ≤ 0 (combined with rk ≥ 0), we conclude that ||Bk v k ||2 − ||Bk−1 v k−1 ||2 ≤ (2 − c)2 ||Bk−1 xt ||2 = (2 − c)2 xt Ak−1 xt , where we set Ak = Bk Bk . A simple4 (and crude) upper bound on the last term follows by observing that ||xt || = 1 implies xt Ak−1 xt ≤ ||Ak−1 ||, the spectral norm (largest eigenvalue) of Ak−1 . Since a factor matrix of the form (I − ρ xx ) with ρ ≤ 1 and ||x|| = 1 has k−1 spectral norm one, we have xt Ak−1 xt ≤ ||Ak−1 || ≤ i=1 ||I − ρi xt(i) xt(i) ||2 ≤ 1. Therefore, summing over k = 1, . . . , m = |M| (or, equivalently, over t ∈ M) and using v 0 = 0 yields the upper bound ||Bm v m ||2 ≤ (2 − c)2 m. (3) To ﬁnd a lower bound of the left-hand side of (3), we ﬁrst pick any unit norm vector u ∈ Rn , and apply the standard Cauchy-Schwartz inequality: ||Bm v m || ≥ u Bm v m . Then, we observe that for a generic trial t = t(k) the update rule of our algorithm allows us to write u Bk v k − u Bk−1 v k−1 = rk yt u Bk−1 xt ≥ rk (γ − max{0, γ − yt u Bk−1 xt }), where the last inequality follows from rk ≥ 0 and holds for any margin value γ > 0. We sum 3 ˜ The subscript c in Sc emphasizes the dependence of the transformed sequence on the choice of c. Note that in the special case c = 1 we have ρk = 0 for any k and α = 1, thereby recovering the standard Perceptron bound for nonseparable sequences (see, e.g., [12]). 4 A slightly more reﬁned bound can be derived which depends on the trace of matrices I − Ak . Details will be given in the full version of this paper. the above over k = 1, . . . , m and exploit c ≤ rk ≤ 2 − c after rearranging terms. This gets ˜ ||Bm v m || ≥ u Bm v m ≥ c γ m − (2 − c)Dγ (u; Sc ). Combining with (3) and solving for m gives the claimed bound. From the above result one can see that our algorithm might be viewed as a standard Perceptron ˜ algorithm operating on the transformed sequence Sc in (1). We now give a qualitative argument, ˜ which is suggestive of the improved margin properties of Sc . Assume for simplicity that all examples (xt , yt ) in the original sequence are correctly classiﬁed by hyperplane u with the same margin γ = yt u xt > 0, where t = t(k). According to Theorem 1, the parameters ρ1 , ρ2 , . . . should be small positive numbers. Assume, again for simplicity, that all ρk are set to the same small enough k value ρ > 0. Then, up to ﬁrst order, matrix Bk = i=1 (I − ρ xt(i) xt(i) ) can be approximated as Bk I −ρ k i=1 xt(i) xt(i) . Then, to the extent that the above approximation holds, we can write:5 yt u Bk−1 xt = yt u I −ρ = yt u xt − ρ yt k−1 i=1 xt(i) xt(i) xt = yt u k−1 i=1 I −ρ k−1 i=1 yt(i) xt(i) yt(i) xt(i) xt yt(i) u xt(i) yt(i) xt(i) xt = γ − ρ γ yt v k−1 xt . Now, yt v k−1 xt is the margin of the (ﬁrst-order) Perceptron vector v k−1 over a mistaken trial for the Higher-order Perceptron vector wk−1 . Since the two vectors v k−1 and wk−1 are correlated (recall that v k−1 wk−1 = v k−1 Bk−1 Bk−1 v k−1 = ||Bk−1 v k−1 ||2 ≥ 0) the mistaken condition yt wk−1 xt ≤ 0 is more likely to imply yt v k−1 xt ≤ 0 than the opposite. This tends to yield a margin larger than the original margin γ. As we mentioned in Section 2, this is also advantageous from a computational standpoint, since in those cases the matrix update Bk−1 → Bk might be skipped (this is equivalent to setting ρk = 0), still Theorem 1 would hold. Though the above might be the starting point of a more thorough theoretical understanding of the margin properties of our algorithm, in this paper we prefer to stop early and leave any further investigation to collecting experimental evidence. 4 Experiments We tested the empirical performance of our algorithm by conducting a number of experiments on a collection of datasets, both artiﬁcial and real-world from diverse domains (Optical Character Recognition, text categorization, DNA microarrays). The main goal of these experiments was to compare Higher-order Perceptron (with both p = 2 and p > 2) to known Perceptron-like algorithms, such as ﬁrst-order [21] and second-order Perceptron [7], in terms of training accuracy (i.e., convergence speed) and test set accuracy. The results are contained in Tables 1, 2, 3, and in Figure 2. Task 1: DNA microarrays and artiﬁcial data. The goal here was to test the convergence properties of our algorithms on sparse target learning tasks. We ﬁrst tested on a couple of well-known DNA microarray datasets. For each dataset, we ﬁrst generated a number of random training/test splits (our random splits also included random permutations of the training set). The reported results are averaged over these random splits. The two DNA datasets are: i. The ER+/ER− dataset from [14]. Here the task is to analyze expression proﬁles of breast cancer and classify breast tumors according to ER (Estrogen Receptor) status. This dataset (which we call the “Breast” dataset) contains 58 expression proﬁles concerning 3389 genes. We randomly split 1000 times into a training set of size 47 and a test set of size 11. ii. The “Lymphoma” dataset [1]. Here the goal is to separate cancerous and normal tissues in a large B-Cell lymphoma problem. The dataset contains 96 expression proﬁles concerning 4026 genes. We randomly split the dataset into a training set of size 60 and a test set of size 36. Again, the random split was performed 1000 times. On both datasets, the tested algorithms have been run by cycling 5 times over the current training set. No kernel functions have been used. We also artiﬁcially generated two (moderately) sparse learning problems with margin γ ≥ 0.005 at labeling noise levels η = 0.0 (linearly separable) and η = 0.1, respectively. The datasets have been generated at random by ﬁrst generating two (normalized) target vectors u ∈ {−1, 0, +1}500 , where the ﬁrst 50 components are selected independently at random in {−1, +1} and the remaining 450 5 Again, a similar argument holds in the more general setting p ≥ 2. The reader should notice how important the dependence of Bk on the past is to this argument. components are 0. Then we set η = 0.0 for the ﬁrst target and η = 0.1 for the second one and, corresponding to each of the two settings, we randomly generated 1000 training examples and 1000 test examples. The instance vectors are chosen at random from [−1, +1]500 and then normalized. If u · xt ≥ γ then a +1 label is associated with xt . If u · xt ≤ −γ then a −1 label is associated with xt . The labels so obtained are ﬂipped with probability η. If |u · xt | < γ then xt is rejected and a new vector xt is drawn. We call the two datasets ”Artiﬁcial 0.0 ” and ”Artiﬁcial 0.1 ”. We tested our algorithms by training over an increasing number of epochs and checking the evolution of the corresponding test set accuracy. Again, no kernel functions have been used. Task 2: Text categorization. The text categorization datasets are derived from the ﬁrst 20,000 newswire stories in the Reuters Corpus Volume 1 (RCV1, [22]). A standard TF - IDF bag-of-words encoding was used to transform each news story into a normalized vector of real attributes. We built four binary classiﬁcation problems by “binarizing” consecutive news stories against the four target categories 70, 101, 4, and 59. These are the 2nd, 3rd, 4th, and 5th most frequent6 categories, respectively, within the ﬁrst 20,000 news stories of RCV1. We call these datasets RCV1x , where x = 70, 101, 4, 59. Each dataset was split into a training set of size 10,000 and a test set of the same size. All algorithms have been trained for a single epoch. We initially tried polynomial kernels, then realized that kernel functions did not signiﬁcantly alter our conclusions on this task. Thus the reported results refer to algorithms with no kernel functions. Task 3: Optical character recognition (OCR). We used two well-known OCR benchmarks: the USPS dataset and the MNIST dataset [16] and followed standard experimental setups, such as the one in [9], including the one-versus-rest scheme for reducing a multiclass problem to a set of binary tasks. We used for each algorithm the standard Gaussian and polynomial kernels, with parameters chosen via 5-fold cross validation on the training set across standard ranges. Again, all algorithms have been trained for a single epoch over the training set. The results in Table 3 only refer to the best parameter settings for each kernel. Algorithms. We implemented the standard Perceptron algorithm (with and without kernels), the Second-order Perceptron algorithm, as described in [7] (with and without kernels), and our Higherorder Perceptron algorithm. The implementation of the latter algorithm (for both p = 2 and p > 2) was ”implicit primal” when tested on the sparse learning tasks, and in dual variables for the other two tasks. When using Second-order Perceptron, we set its parameter a (see [7] for details) by testing on a generous range of values. For brevity, only the settings achieving the best results are reported. On the sparse learning tasks we tried Higher-order Perceptron with norm p = 2, 4, 7, 10, while on the other two tasks we set p = 2. In any case, for each value of p, we set7 ρk = c/k, with c = 0, 0.2, 0.4, 0.6, 0.8. Since c = 0 corresponds to a standard p-norm Perceptron algorithm [13, 12] we tried to emphasize the comparison c = 0 vs. c > 0. Finally, when using kernels on the OCR tasks, we also compared to a sparse dual version of Higher-order Perceptron. On a mistaken round t = t(k), this algorithm sets ρk = c/k if yt v k−1 xt ≥ 0, and ρk = 0 otherwise (thus, when yt v k−1 xt < 0 the matrix Bk−1 is not updated). For the sake of brevity, the standard Perceptron algorithm is called FO (”First Order”), the Second-order algorithm is denoted by SO (”Second Order”), while the Higher-order algorithm with norm parameter p and ρk = c/k is abbreviated as HOp (c). Thus, for instance, FO = HO2 (0). Results and conclusions. Our Higher-order Perceptron algorithm seems to deliver interesting results. In all our experiments HOp (c) with c > 0 outperforms HOp (0). On the other hand, the comparison HOp (c) vs. SO depends on the speciﬁc task. On the DNA datasets, HOp (c) with c > 0 is clearly superior in Breast. On Lymphoma, HOp (c) gets worse as p increases. This is a good indication that, in general, a multiplicative algorithm is not suitable for this dataset. In any case, HO2 turns out to be only slightly worse than SO. On the artiﬁcial datasets HOp (c) with c > 0 is always better than the corresponding p-norm Perceptron algorithm. On the text categorization tasks, HO2 tends to perform better than SO. On USPS, HO2 is superior to the other competitors, while on MNIST it performs similarly when combined with Gaussian kernels (though it turns out to be relatively sparser), while it is slightly inferior to SO when using polynomial kernels. The sparse version of HO2 cuts the matrix updates roughly by half, still maintaining a good performance. In all cases HO2 (either sparse or not) signiﬁcantly outperforms FO. In conclusion, the Higher-order Perceptron algorithm is an interesting tool for on-line binary clas6 7 We did not use the most frequent category because of its signiﬁcant overlap with the other ones. Notice that this setting fulﬁlls the condition on ρk stated in Theorem 1. Table 1: Training and test error on the two datasets ”Breast” and ”Lymphoma”. Training error is the average total number of updates over 5 training epochs, while test error is the average fraction of misclassiﬁed patterns in the test set, The results refer to the same training/test splits. For each algorithm, only the best setting is shown (best training and best test setting coincided in these experiments). Thus, for instance, HO2 differs from FO because of the c parameter. We emphasized the comparison HO7 (0) vs. HO7 (c) with best c among the tested values. According to Wilcoxon signed rank test, an error difference of 0.5% or larger might be considered signiﬁcant. In bold are the smallest ﬁgures achieved on each row of the table. FO TRAIN TEST TRAIN TEST LYMPHOMA HO 2 HO 4 HO 7 (0) HO 7 HO 10 SO 45.2 23.4% 22.1 11.8% 21.7 16.4% 19.6 10.0% 24.5 13.3% 18.9 10.0% 47.4 15.7% 23.0 11.5% 24.5 12.0% 20.0 11.5% 32.4 13.5 23.1 11.9% 29.6 15.0% 19.3 9.6% FO = HO 2(0.0) Training updates vs training epochs on Artificial 0.0 SO # of training updates 800 * HO 4(0.4) 600 HO 7(0.0) * 400 300 * * * * * SO 2400 HO 2(0.4) 700 500 HO 7 (0.4) * 2000 * 1200 400 2 3 5 10 15 20 * 1 * * 2 3 * Test error rates * * * (a = 0.2) HO 2(0.4) HO 4(0.4) * * * * HO 7(0.0) HO 7 (0.4) 14% Test error rates (minus 10%) FO = HO 2(0.0) SO 18% * HO 7(0.0) HO 7(0.4) 5 10 15 20 # of training epochs Test error rates vs training epochs on Artificial 0.0 22% * * # of training epochs 26% (a = 0.2) HO 2(0.4) HO 4(0.4) 1600 800 * 1 FO = HO 2(0.0) Training updates vs training epochs on Artificial 0.1 (a = 0.2) # of training updates B REAST FO = HO 2(0.0) Test error rates vs training epochs on Artificial 0.1 SO 26% 22% * * * * * * * * (a = 0.2) HO 2(0.4) HO 4(0.4) 18% HO 7(0.0) 14% HO 7 (0.4) 10% 10% 6% 6% 1 2 3 5 10 # of training epochs 15 20 1 2 3 5 10 15 20 # of training epochs Figure 2: Experiments on the two artiﬁcial datasets (Artiﬁcial0.0 , on the left, and Artiﬁcial0.1 , on the right). The plots give training and test behavior as a function of the number of training epochs. Notice that the test set in Artiﬁcial0.1 is affected by labelling noise of rate 10%. Hence, a visual comparison between the two plots at the bottom can only be made once we shift down the y-axis of the noisy plot by 10%. On the other hand, the two training plots (top) are not readily comparable. The reader might have difﬁculty telling apart the two kinds of algorithms HOp (0.0) and HOp (c) with c > 0. In practice, the latter turned out to be always slightly superior in performance to the former. siﬁcation, having the ability to combine multiplicative (or nonadditive) and second-order behavior into a single inference procedure. Like other algorithms, HOp can be extended (details omitted due to space limitations) in several ways through known worst-case learning technologies, such as large margin (e.g., [17, 11]), label-efﬁcient/active learning (e.g., [5, 8]), and bounded memory (e.g., [10]). References [1] A. Alizadeh, et al. (2000). Distinct types of diffuse large b-cell lymphoma identiﬁed by gene expression proﬁling. Nature, 403, 503–511. [2] D. Angluin (1988). Queries and concept learning. Machine Learning, 2(4), 319–342. [3] P. Auer & M.K. Warmuth (1998). Tracking the best disjunction. Machine Learning, 32(2), 127–150. [4] K.S. Azoury & M.K. Warmuth (2001). Relative loss bounds for on-line density estimation with the exponential familiy of distributions. Machine Learning, 43(3), 211–246. [5] A. Bordes, S. Ertekin, J. Weston, & L. Bottou (2005). Fast kernel classiﬁers with on-line and active learning. JMLR, 6, 1579–1619. [6] N. Cesa-Bianchi, Y. Freund, D. Haussler, D.P. Helmbold, R.E. Schapire, & M.K. Warmuth (1997). How to use expert advice. J. ACM, 44(3), 427–485. Table 2: Experimental results on the four binary classiﬁcation tasks derived from RCV1. ”Train” denotes the number of training corrections, while ”Test” gives the fraction of misclassiﬁed patterns in the test set. Only the results corresponding to the best test set accuracy are shown. In bold are the smallest ﬁgures achieved for each of the 8 combinations of dataset (RCV1x , x = 70, 101, 4, 59) and phase (training or test). FO TRAIN TEST 993 673 803 767 7.20% 6.39% 6.14% 6.45% RCV170 RCV1101 RCV14 RCV159 HO 2 TRAIN TEST 941 665 783 762 6.83% 5.81% 5.94% 6.04% SO TRAIN TEST 880 677 819 760 6.95% 5.48% 6.05% 6.84% Table 3: Experimental results on the OCR tasks. ”Train” denotes the total number of training corrections, summed over the 10 categories, while ”Test” denotes the fraction of misclassiﬁed patterns in the test set. Only the results corresponding to the best test set accuracy are shown. For the sparse version of HO2 we also reported (in parentheses) the number of matrix updates during training. In bold are the smallest ﬁgures achieved for each of the 8 combinations of dataset (USPS or MNIST), kernel type (Gaussian or Polynomial), and phase (training or test). FO TRAIN U SPS M NIST G AUSS P OLY G AUSS P OLY TEST 1385 1609 5834 8148 6.53% 7.37% 2.10% 3.04% HO 2 TRAIN TEST 945 1090 5351 6404 4.76% 5.71% 1.79% 2.27% Sparse HO2 SO TRAIN TEST TRAIN TEST 965 (440) 1081 (551) 5363 (2596) 6476 (3311) 5.13% 5.52% 1.81% 2.28% 1003 1054 5684 6440 5.05% 5.53% 1.82% 2.03% [7] N. Cesa-Bianchi, A. Conconi & C. Gentile (2005). A second-order perceptron algorithm. SIAM Journal of Computing, 34(3), 640–668. [8] N. Cesa-Bianchi, C. Gentile, & L. Zaniboni (2006). Worst-case analysis of selective sampling for linearthreshold algorithms. JMLR, 7, 1205–1230. [9] C. Cortes & V. Vapnik (1995). Support-vector networks. Machine Learning, 20(3), 273–297. [10] O. Dekel, S. Shalev-Shwartz, & Y. Singer (2006). The Forgetron: a kernel-based Perceptron on a ﬁxed budget. NIPS 18, MIT Press, pp. 259–266. [11] C. Gentile (2001). A new approximate maximal margin classiﬁcation algorithm. JMLR, 2, 213–242. [12] C. Gentile (2003). The Robustness of the p-norm Algorithms. Machine Learning, 53(3), pp. 265–299. [13] A.J. Grove, N. Littlestone & D. Schuurmans (2001). General convergence results for linear discriminant updates. Machine Learning Journal, 43(3), 173–210. [14] S. Gruvberger, et al. (2001). Estrogen receptor status in breast cancer is associated with remarkably distinct gene expression patterns. Cancer Res., 61, 5979–5984. [15] J. Kivinen, M.K. Warmuth, & P. Auer (1997). The perceptron algorithm vs. winnow: linear vs. logarithmic mistake bounds when few input variables are relevant. Artiﬁcial Intelligence, 97, 325–343. [16] Y. Le Cun, et al. (1995). Comparison of learning algorithms for handwritten digit recognition. ICANN 1995, pp. 53–60. [17] Y. Li & P. Long (2002). The relaxed online maximum margin algorithm. Machine Learning, 46(1-3), 361–387. [18] N. Littlestone (1988). Learning quickly when irrelevant attributes abound: a new linear-threshold algorithm. Machine Learning, 2(4), 285–318. [19] N. Littlestone & M.K. Warmuth (1994). The weighted majority algorithm. Information and Computation, 108(2), 212–261. [20] P. Long & X. Wu (2004). Mistake bounds for maximum entropy discrimination. NIPS 2004. [21] A.B.J. Novikov (1962). On convergence proofs on perceptrons. Proc. of the Symposium on the Mathematical Theory of Automata, vol. XII, pp. 615–622. [22] Reuters: 2000. http://about.reuters.com/researchandstandards/corpus/. [23] S. Shalev-Shwartz & Y. Singer (2006). Online Learning Meets Optimization in the Dual. COLT 2006, pp. 423–437. [24] B. Schoelkopf & A. Smola (2002). Learning with kernels. MIT Press. [25] Vovk, V. (2001). Competitive on-line statistics. International Statistical Review, 69, 213-248.

5 0.15798862 48 nips-2007-Collective Inference on Markov Models for Modeling Bird Migration

Author: M.a. S. Elmohamed, Dexter Kozen, Daniel R. Sheldon

Abstract: We investigate a family of inference problems on Markov models, where many sample paths are drawn from a Markov chain and partial information is revealed to an observer who attempts to reconstruct the sample paths. We present algorithms and hardness results for several variants of this problem which arise by revealing different information to the observer and imposing different requirements for the reconstruction of sample paths. Our algorithms are analogous to the classical Viterbi algorithm for Hidden Markov Models, which ﬁnds the single most probable sample path given a sequence of observations. Our work is motivated by an important application in ecology: inferring bird migration paths from a large database of observations. 1

6 0.1485358 34 nips-2007-Bayesian Policy Learning with Trans-Dimensional MCMC

7 0.13084976 21 nips-2007-Adaptive Online Gradient Descent

8 0.11673644 206 nips-2007-Topmoumoute Online Natural Gradient Algorithm

9 0.11164635 47 nips-2007-Collapsed Variational Inference for HDP

10 0.10502413 203 nips-2007-The rat as particle filter

11 0.10148083 145 nips-2007-On Sparsity and Overcompleteness in Image Models

12 0.099139884 59 nips-2007-Continuous Time Particle Filtering for fMRI

13 0.094071202 125 nips-2007-Markov Chain Monte Carlo with People

14 0.091028407 91 nips-2007-Fitted Q-iteration in continuous action-space MDPs

15 0.089454323 87 nips-2007-Fast Variational Inference for Large-scale Internet Diagnosis

16 0.087181732 148 nips-2007-Online Linear Regression and Its Application to Model-Based Reinforcement Learning

17 0.083650783 127 nips-2007-Measuring Neural Synchrony by Message Passing

18 0.083120376 163 nips-2007-Receding Horizon Differential Dynamic Programming

19 0.081543542 140 nips-2007-Neural characterization in partially observed populations of spiking neurons

20 0.080685869 170 nips-2007-Robust Regression with Twinned Gaussian Processes

similar papers computed by lsi model

lsi for this paper:

topicId topicWeight

[(0, -0.245), (1, -0.066), (2, 0.083), (3, -0.025), (4, -0.001), (5, -0.048), (6, -0.089), (7, -0.11), (8, -0.197), (9, -0.239), (10, 0.138), (11, 0.124), (12, -0.086), (13, 0.15), (14, 0.03), (15, 0.155), (16, -0.15), (17, 0.041), (18, -0.046), (19, -0.079), (20, 0.166), (21, -0.007), (22, -0.115), (23, 0.079), (24, -0.004), (25, -0.022), (26, 0.013), (27, 0.082), (28, -0.034), (29, -0.08), (30, -0.054), (31, 0.024), (32, 0.119), (33, -0.069), (34, 0.056), (35, 0.068), (36, 0.089), (37, 0.056), (38, 0.05), (39, -0.092), (40, -0.111), (41, -0.022), (42, -0.097), (43, 0.001), (44, 0.102), (45, 0.026), (46, -0.045), (47, -0.032), (48, -0.04), (49, 0.064)]

similar papers list:

simIndex simValue paperId paperTitle

same-paper 1 0.96800292 213 nips-2007-Variational Inference for Diffusion Processes

Author: Cédric Archambeau, Manfred Opper, Yuan Shen, Dan Cornford, John S. Shawe-taylor

2 0.77454448 214 nips-2007-Variational inference for Markov jump processes

Author: Manfred Opper, Guido Sanguinetti

3 0.62976241 17 nips-2007-A neural network implementing optimal state estimation based on dynamic spike train decoding

Author: Omer Bobrowski, Ron Meir, Shy Shoham, Yonina Eldar

4 0.61929625 48 nips-2007-Collective Inference on Markov Models for Modeling Bird Migration

Author: M.a. S. Elmohamed, Dexter Kozen, Daniel R. Sheldon

5 0.59584075 146 nips-2007-On higher-order perceptron algorithms

Author: Claudio Gentile, Fabio Vitale, Cristian Brotto

6 0.5264765 59 nips-2007-Continuous Time Particle Filtering for fMRI

7 0.52283853 87 nips-2007-Fast Variational Inference for Large-scale Internet Diagnosis

8 0.49791861 206 nips-2007-Topmoumoute Online Natural Gradient Algorithm

9 0.46367663 203 nips-2007-The rat as particle filter

10 0.44934052 34 nips-2007-Bayesian Policy Learning with Trans-Dimensional MCMC

11 0.44713601 47 nips-2007-Collapsed Variational Inference for HDP

12 0.41653916 167 nips-2007-Regulator Discovery from Gene Expression Time Series of Malaria Parasites: a Hierachical Approach

13 0.40020308 21 nips-2007-Adaptive Online Gradient Descent

14 0.39748916 91 nips-2007-Fitted Q-iteration in continuous action-space MDPs

15 0.38859981 68 nips-2007-Discovering Weakly-Interacting Factors in a Complex Stochastic Process

16 0.38800225 28 nips-2007-Augmented Functional Time Series Representation and Forecasting with Gaussian Processes

17 0.37850872 31 nips-2007-Bayesian Agglomerative Clustering with Coalescents

18 0.37242836 170 nips-2007-Robust Regression with Twinned Gaussian Processes

19 0.34310824 163 nips-2007-Receding Horizon Differential Dynamic Programming

20 0.34037861 125 nips-2007-Markov Chain Monte Carlo with People

similar papers computed by lda model

lda for this paper:

topicId topicWeight

[(5, 0.032), (13, 0.031), (16, 0.06), (18, 0.386), (21, 0.064), (31, 0.02), (34, 0.044), (35, 0.02), (47, 0.087), (49, 0.011), (83, 0.082), (85, 0.021), (87, 0.015), (90, 0.051)]

similar papers list:

simIndex simValue paperId paperTitle

1 0.92673445 203 nips-2007-The rat as particle filter

Author: Aaron C. Courville, Nathaniel D. Daw

Abstract: Although theorists have interpreted classical conditioning as a laboratory model of Bayesian belief updating, a recent reanalysis showed that the key features that theoretical models capture about learning are artifacts of averaging over subjects. Rather than learning smoothly to asymptote (reﬂecting, according to Bayesian models, the gradual tradeoff from prior to posterior as data accumulate), subjects learn suddenly and their predictions ﬂuctuate perpetually. We suggest that abrupt and unstable learning can be modeled by assuming subjects are conducting inference using sequential Monte Carlo sampling with a small number of samples — one, in our simulations. Ensemble behavior resembles exact Bayesian models since, as in particle ﬁlters, it averages over many samples. Further, the model is capable of exhibiting sophisticated behaviors like retrospective revaluation at the ensemble level, even given minimally sophisticated individuals that do not track uncertainty in their beliefs over trials. 1

2 0.88455522 3 nips-2007-A Bayesian Model of Conditioned Perception

Author: Alan Stocker, Eero P. Simoncelli

Abstract: unkown-abstract

same-paper 3 0.8467049 213 nips-2007-Variational Inference for Diffusion Processes

Author: Cédric Archambeau, Manfred Opper, Yuan Shen, Dan Cornford, John S. Shawe-taylor

4 0.53657132 125 nips-2007-Markov Chain Monte Carlo with People

Author: Adam Sanborn, Thomas L. Griffiths

Abstract: Many formal models of cognition implicitly use subjective probability distributions to capture the assumptions of human learners. Most applications of these models determine these distributions indirectly. We propose a method for directly determining the assumptions of human learners by sampling from subjective probability distributions. Using a correspondence between a model of human choice and Markov chain Monte Carlo (MCMC), we describe a method for sampling from the distributions over objects that people associate with different categories. In our task, subjects choose whether to accept or reject a proposed change to an object. The task is constructed so that these decisions follow an MCMC acceptance rule, deﬁning a Markov chain for which the stationary distribution is the category distribution. We test this procedure for both artiﬁcial categories acquired in the laboratory, and natural categories acquired from experience. 1

5 0.53279215 34 nips-2007-Bayesian Policy Learning with Trans-Dimensional MCMC

Author: Matthew Hoffman, Arnaud Doucet, Nando D. Freitas, Ajay Jasra

Abstract: A recently proposed formulation of the stochastic planning and control problem as one of parameter estimation for suitable artiﬁcial statistical models has led to the adoption of inference algorithms for this notoriously hard problem. At the algorithmic level, the focus has been on developing Expectation-Maximization (EM) algorithms. In this paper, we begin by making the crucial observation that the stochastic control problem can be reinterpreted as one of trans-dimensional inference. With this new interpretation, we are able to propose a novel reversible jump Markov chain Monte Carlo (MCMC) algorithm that is more efﬁcient than its EM counterparts. Moreover, it enables us to implement full Bayesian policy search, without the need for gradients and with one single Markov chain. The new approach involves sampling directly from a distribution that is proportional to the reward and, consequently, performs better than classic simulations methods in situations where the reward is a rare event.

6 0.49057621 47 nips-2007-Collapsed Variational Inference for HDP

7 0.47889355 153 nips-2007-People Tracking with the Laplacian Eigenmaps Latent Variable Model

8 0.47005525 51 nips-2007-Comparing Bayesian models for multisensory cue combination without mandatory integration

9 0.46846005 93 nips-2007-GRIFT: A graphical model for inferring visual classification features from human data

10 0.46157143 74 nips-2007-EEG-Based Brain-Computer Interaction: Improved Accuracy by Automatic Single-Trial Error Detection

11 0.45821261 214 nips-2007-Variational inference for Markov jump processes

12 0.45416462 202 nips-2007-The discriminant center-surround hypothesis for bottom-up saliency

13 0.45187527 87 nips-2007-Fast Variational Inference for Large-scale Internet Diagnosis

14 0.43993998 100 nips-2007-Hippocampal Contributions to Control: The Third Way

15 0.43960705 155 nips-2007-Predicting human gaze using low-level saliency combined with face detection

16 0.43911332 140 nips-2007-Neural characterization in partially observed populations of spiking neurons

17 0.43670753 104 nips-2007-Inferring Neural Firing Rates from Spike Trains Using Gaussian Processes

18 0.43222666 122 nips-2007-Locality and low-dimensions in the prediction of natural experience from fMRI

19 0.42738882 154 nips-2007-Predicting Brain States from fMRI Data: Incremental Functional Principal Component Regression

20 0.42569235 48 nips-2007-Collective Inference on Markov Models for Modeling Bird Migration