Competitive Denoising Methods with Similar Performances

Overview

B We show that "weakly" denoised versions of the input image can serve to compute a performant patch-based aggregated estimator. B We evaluate the performance of each patch estimator to compute the Exponentially Weighted Aggregation (EWA) (Leung & Barron, 2006) (Dalayan & Tsybakov, 2008) (Salmon & Le Pennec, 2009). B The aggregation method is flexible enough to combine any standard denoising algorithms and has an interpretation with Gibbs distribution. B PEWA is based on a MCMC sampling and is able to produce results that are comparable to the current state-of-the-art. PEWA: A statistical aggregation method which combines denoised image patches, generalizes the NL-means and produces state-of-the-art results Similar ideas:

Image Patch Model and Estimator

B Notations (patches):

f (x): n-dimensional unknown patch at location

x 2 X ⇢ R 2 v(x): n-dimensional noisy patch (Gaussian noise) v(x) = f (x) + "(x) with "(x) ⇠ N (0, 2 I n⇥n ) b f (x): n-dimensional patch estimator B Empirical statistic (measure): detection of deviation from noise

R( b f (x)) = kv(x) b f (x)k 2 n n 2
with this choice, we have

E[R( b f (x))] = E[kf (x) b f (x)k 2 n ] (L 2 risk)

Aggregation by Exponential Weights

B Assume a set {f (x), 2 ⇤ = {1, • • • , M }} of pre-computed estimators.

We consider an aggregate that is the weighted average of estimators with some data-dependent weights:

b f (x) = M X =1 w (x) f (x) such that w (x) 0 and M X =1 w (x) = 1.
B We associate two probability measures w(x) = {w (x)} and ⇡(x) = {⇡ (x)}, 2 ⇤ and we define the Kullback-Leibler divergence as [START_REF] Rigollet | [END_REF]:

D KL (w(x), ⇡(x)) = M X =1 w (x) log ✓ w (x) ⇡ (x) ◆ .
The role of the distribution ⇡ is to put a prior weight on the estimators in the set.

M X =1 b (x)w (x)
)

where ↵, > 0, b (x)w (x) = 0 and : R ! R + .

B From the Karush-Kuhn-Tucker conditions, we get

b w (x) = exp( (R(f (x)))/ ) ⇡ (x) P M 0 =1 exp( (R(f 0 (x)))/ ) ⇡ 0 (x)
, and can be interpreted as a "temperature" parameter.

PEWA: Patch-based EWA estimator

B Assume {u 1 , • • • , u L } "
weakly" denoised versions of v (Gaussian, Wiener, DCT, Wavelet, Median, Bilateral ...).

B An estimator f (x) is a n-dimensional patch (denoted u l (y)) taken in u `, `2 {1, • • • , L} at any location y 2 X .
B Our estimator is of the following form:

b f (x) = 1 Z(x) L X `=1 X y2X e |R(u `(y))|/ ⇡ `(y) u `(y)
where Z(x) is a normalization constant and (z) = |z| favors estimators with small deviations from noise.

Patch-based EWA Estimator:

Gibbs Model and "Neighborhood" Prior

The patch-based EWA (PEWA) estimator is written in terms of Gibbs distributions as ( = 4 2 , Leung, 2006):

b f PEWA (x) = 1 Z(x) L X `=1 X y2X e E(u `(y)) u `(y) Z(x) = L X `0=1 X y 0 2X e E(u `0 (y 0 )) , E(u `(y)) = |kv(x) u `(y)k 2 n n 2 | 4 2 + kx yk 2 2 2⌧ 2 .
B Prior: favors patches located in the spatial neighborhood of x.

B Monte-Carlo sampling: method to approximately compute PEWA when the number of patch estimators is large.

Patch-based EWA Estimator: "Data-driven" Monte-Carlo Sampling:

Generalization of Non-Local means

PEWA is equivalent to NL-means if we choose L = 1, u `=1 = v, (z) = z and ⌧ ! 1 (⇡ flat prior): b f NLM (x) = 1 Z(x) X y2X e E(v(y)) v(y), Z(x) = X y 0 2X e E(v(y 0 )) E(v(y)) = kv(x) v(y)k 2 n n 2 + kx yk 2 2 2⌧ 2 ⇡ kv(x) v(y)k 2 n + cte

Computational Issues

Assume a random process (F m (x)) m 0 and an initial noisy patch F 0 (x) = v(x). The "data-driven" MCMC procedure is based on the Metropolis-Hastings algorithm: Draw a patch by considering a two-stage drawing procedure:

1. draw uniformly a value `in the set {1, 2, • • • , L}.

draw a pixel y

= y c + , y 2 X , with ⇠ N (0, I 2⇥2 ⌧ 2 ) and y c is the position of the current patch. At the initialization y c = x. Define F m+1 (x) = ⇢ u `(y) if ↵ ⇠ U [0, 1]  e (E(u `(y)) E(F m (x))) F m (x) otherwise.
"Data-driven" Monte-Carlo Sampling (contd')

Computational Issues

If we assume the Markov chain is ergodic, homogeneous, irreductible, reversible and stationary, for any F 0 (x), we have almost surely

lim T !+1 1 T T b T X m=T b F m (x) ⇡ b f PEWA (x)
B MCMC algorithm (one patch):

-T ⇡ 1000 is the maximum number of samples.

-The first T b = 250 samples are discarded (burn-in phase).

B Global image denoising:

-Patch overlapping: multiple estimates b f PEWA (x) at a given pixel x.

-Uniform averaging: "fusion" of n independent Markov chains at each pixel.

Denoising: Additive White Gaussian Noise on Natural Images

Set of 25 images. Top left: images from the BM3D website (cs.tut.fi/˜foi/GCFBM3D/); Bottom left: images from IPOL (ipol.im); Right: images from the Berkeley segmentation database (eecs.berkeley.edu/Research/Projects/CS/ vision/bsds/).

Experimental Results: Artificially Noisy Data

A two-step procedure with the parameters ⌧ = 7, n = 7 ⇥ 7 and L = 4: B 1 st iteration: estimation using the noisy image v and 3 denoised images u l (DCT shrinkage thresholds: {1.25; 1.50; 1.75} ⇥ ) [START_REF] Yu | Image Processing On Line[END_REF]. B 2 nd iteration: estimation as before using the 1 st PEWA estimator considered as an additional denoised image (improvement in the range of 0.2 to 0.5 dB).

Experimental Results:

Artificially Noisy "Lena" Image

Experimental Results:

Artificially Noisy "Cameraman" Image Figure 2: "Barabra" image corrupted with white Gaussian noise ( = 20). Left: MCMC-based PEWA (1000 samples) applied to 7 ⇥ 7 non-overlapping patches. Right: "exact" PEWA applied to 7 ⇥ 7 non-overlapping patches and inspection of all image patches in the image (L ⇥ |X | patches). BM3D [Dabov, 2007] (PSNR = 29.19) NL-Bayes [Lebrun, 2013] (PSNR = 29.22)

Cameraman (z) = exp(z) (z) = z 2 (z) = z (z) = |z| (z) = H 4 (z) (z) = p z (z) = (z) + ( = 20) = 0.

Experimental Results

Artificially Noisy Data noisy (PSNR = 24.61) PEWA (PSNR = 29.25)

BM3D [Dabov, 2007] (PSNR = 29.19) NL-Bayes [Lebrun, 2013] (PSNR = 29.22) BM3D [Dabov, 2007] (PSNR = 29.36)

NL-Bayes [Lebrun, 2013] (PSNR = 29.48)

Experimental Results

Artificially Noisy Data noisy (PSNR = 20.18) PEWA (PSNR = 29.49) BM3D [Dabov, 2007] (PSNR = 29.36)

NL-Bayes [Lebrun, 2013] (PSNR = 29.48) The PSNR values of the 3 images denoised with DCT-based transform (Yu, 2011) and combined with PEWA are 25.77, 24.26 and 22.85.

  Non-local means (practice): selection of patches in a fixed-size search window (21 ⇥ 21 pixels).

Figure 1 :

 1 Figure 1: "Lena" image corrupted with white Gaussian noise ( = 20). Left: MCMC-based PEWA (1000 samples) applied to 7 ⇥ 7 non-overlapping patches. Right: "exact" PEWA applied to 7 ⇥ 7 non-overlapping patches and inspection of all image patches in the image (L ⇥ |X | patches).

Figure 3 :

 3 Figure 3: "Valldemossa" image corrupted with white Gaussian noise ( = 15). The PSNR values of the 3 images denoised with DCT-based transform (Yu, 2011) and combined with PEWA are 27.78, 27.04 and 26.26.

Figure 4 :

 4 Figure 4: "Valldemossa" image corrupted with white Gaussian noise ( = 15).The PSNR values of the 3 images denoised with DCT-based transform[START_REF] Yu | Image Processing On Line[END_REF] and combined with PEWA are 27.78, 27.04 and 26.26.

Figure 5 :

 5 Figure 5: Castle image corrupted with white Gaussian noise ( = 25).The PSNR values of the 3 images denoised with DCT-based transform(Yu, 2011) and combined with PEWA are 25.77, 24.26 and 22.85.

Figure 6 :

 6 Figure 6: Castle image corrupted with white Gaussian noise ( = 25).The PSNR values of the 3 images denoised with DCT-based transform(Yu, 2011) and combined with PEWA are 25.77, 24.26 and 22.85.

Table 1 :

 1 Average of denoising results over 25 tested images for several values of (white Gaussian noise). The experiments with NL-Bayes, S-PLE, NL-means and DCT have been performed using the implementations of IPOL (www.ipol.im).

	PEWA 1 st	iteration	= 5 38.27	= 10 34.39	= 15 32.26	= 20 30.76	= 25 29.62	= 50 26.00	= 100 22.35
	PEWA 2 nd	iteration	38.54	34.75	32.67	31.26	30.15	26.95	23.76
	BM3D [Dabov, 2007]	38.64	34.78	32.68	31.25	30.19	26.97	24.08
	NL-Bayes [Lebrun, 2013]	38.60	34.75	32.48	31.22	30.12	26.90	23.65
	S-PLE [Wang, 2013]	38.17	34.38	32.35	30.67	29.77	26.46	23.21
	NL-means [Buades, 2005]	37.44	33.35	31.00	30.16	28.96	25.53	22.29
	DCT [Yu, 2011]	37.81	33.57	31.87	29.95	28.97	25.91	23.08

Table 2 :

 2 PSNR values are averaged over 3 di↵erent noise realizations.

		Experimental Results
	Image	Peppers	House	Lena	Barbara
		(256 ⇥ 256)	(256 ⇥ 256)	(512 ⇥ 512)	(512 ⇥ 512)

Table 3 :

 3 Comparison of several versions of PEWA ((W)iener), (D)CT, Basic) and BM3D on a few standard images corrupted with white Gaussian noise.

Denoising of Real Old Pictures

Gaussian noise and spatially-varying variance:

v(x) = f (x) + "(x), "(x) ⇠ N (0, 2 (x))

Empiricial statistic:

Denoising of Real Old Pictures

Gaussian noise and spatially-varying variance:

Empiricial statistic: Poisson-Gaussian Noise in Fluorescence Microscopy

. B g 0 is gain of the overall electronic system. B @(x) is the number of photo-electrons at pixel x assumed to be Poisson distributed with unknown mean ✓(x). B ✏(x) ⇠ N (0, 2 ✏ ) is a white Gaussian noise and represents "dark current".

Denoising in Fluorescence Microscopy

(single micro-patterned cell)

R( b

Summary

A statistical aggregation method which combines denoised image patches and generalizes the NL-means method:


B Maury

Congested transport at microscopic and macroscopic scales

Keywords: 2010 Mathematics Subject Classification. 34G25, 35R70, 70E55 Continuity equation, crowd motion, optimal transportation, Wasserstein space, granular flows, non-smooth analysis

This note addresses mathematical issues raised by congestion constraints in transport equations that arise in the modeling of crowd motion or more general active entities. We address in particular the differences between the microscopic and the macroscopic settings.

Introduction

We are interested here in the role played by congestion constraints in some evolution processes of active entities, and in the way it modifies the mathematical properties of the associated equations. By congestion we mean the following: at the microscopic level, entities are represented by moving geometrical objects (like discs, or spheres), and congestion simply means that two entities may not occupy the same space at the same time; at the macroscopic level, the population is represented by a density that cannot increase above a prescribed maximal value. Most part of the approach presented here has been developped in the context of human crowd motions, but it may as well apply to other entities like swimming cells (we refer to [START_REF] Sokolov | Concentration dependence of the collective dynamics of swimming bacteria[END_REF] for experiments on bacteria), insects ( [START_REF] Parisi | Faster-is-slower effect in escaping ants revisited: Ants do not behave like humans[END_REF]), or larger animals like sheeps ( [START_REF] Garcimartín | Flow and clogging of a sheep herd passing through a bottleneck[END_REF]). Strategies to account for congestion can be classified according to criteria that are standard in particle physics or granular flow modeling: Soft Sphere / Hard Sphere approaches or, similarly, Molecular Dynamics / Contact Dynamics. In the first (soft) approach, congestion is treated in a smooth way, typically by adding a repulsion force that acts whenever entities become close to each other. If properly tuned, this force prevents a full overlapping of entities (see e.g. [START_REF] Hong | Molecular dynamics simulations of hard sphere granular particles[END_REF], or [START_REF] Helbing | Social force model for pedestrian dynamics[END_REF] for an application of this strategy in the context of crowd motion modeling). In the macroscopic setting it may take different forms. The one that corresponds to short-range repulsion forces consists in adding some sort of pressure within the population, seen here as a compressible fluid (see [START_REF] Degond | Congestion in a macroscopic model of self-driven particles modeling gregariousness[END_REF]). Depending on the context and on the underlying physical reality, other choices are possible and justified. It may for example consist in directly inhibiting the mechanism that is identified as the cause of concentration (see [START_REF] Dolak | The Keller-Segel model with logistic sensitivity function and small diffusivity[END_REF], where the chemotactic velocity is assumed to decrease where local density increases), reducing the desired velocity in crowded areas, in the context of crowd motion ( [START_REF] Hughes | A continuum theory for the flow of pedestrian[END_REF][START_REF] Hughes | The Flow of Human Crowds[END_REF]), adding nonlinear diffusion terms (see e.g. [START_REF] Alexander | Quasi-static evolution and congested crowd transport[END_REF]), or incorporating correction terms that nonlocally depend on the crowd distribution, in order to deviate the trajectory of people to avoid overcrowded areas ( [START_REF] Piccoli | Transport equation with nonlocal velocity in Wasserstein spaces: convergence of numerical schemes[END_REF]).

Let us also mention higher order models that incorporate inertial or delay effects, leading to an evolution equation for velocities, see [START_REF] Bellomo | On the modelling crowd dynamics from scaling to hyperbolic macroscopic models[END_REF] or [START_REF] Twarogowska | Macroscopic modeling and simulations of room evacuation[END_REF].

The present text focuses on the alternative approach, that is accounting for congestion constraints in a hard way. At the microscopic level, it consists in strictly forbidding the overlapping of entities (identified as rigid objects), whereas at the macroscopic level, a set a feasible densities will be defined by prescribing a maximal value for the density. We shall investigate here how these constraints are likely to affect the evolution process that is considered, in the framework of very crude modeling assumptions: entity tendencies correspond to a pure transport process (by the velocity that each of them would like to have if it were alone) that is prevented and modified by congestion.

The present paper puts in perspective different works that have been carried out by the author and various collaborators, it aims in particular at identifying the similarities between the microscopic and macroscopic settings, and more importantly to highlight their deep discrepancies. It focuses on the hard congestion problem, but most remarks that we make here to distinguish the two scales of description are applicabled to soft congestion. Sections 2 and 3 present in parallel the models at both scales, and highlight their formal analogies, while Section 4 asserts their differences, and details how those differences at the mathematical level affect the very behavior of the associated evolution models.

Congested transport at both scales

2.1. Microscopic model. The first problem we shall consider consists in finding a trajectory t -→ q(t) in a euclidean space, starting from a given initial value, verifying dq dt = P Cq U (q), (

where U (q) is a "spontaneous" velocity field, and C q is a set of feasible velocities. The instantaneous projection on C q is intended to force q to remain in a given set K of feasible configurations.

The crowd motion microscopic model (firstly introduced in [START_REF] Maury | A discrete contact model for crowd motion[END_REF]) fits in this framework: individual are identified with rigid discs of radius r > 0 evolving in a two-dimensional space, so that the position vector q = (q 1 , . . . , q N ) ∈ R 2N is subject to remain in

K = q ∈ R 2N , D ij (q) = |q j -q i | -2r ≥ 0 ∀i < j . (2.2)
Feasible velocities are such that, when there is a contact (D ij = 0), the distance may not be further reduced. Thus, the set of feasible velocities is straightforwardly defined as

C q = v ∈ R 2N , D ij (q) = 0 ⇒ G ij • v ≥ 0 , (2.3)
where

G ij = ∇D ij . Given a vector of desired velocities U = (U 1 , . . . , U N ) ∈ R 2N
(where U i is the velocity that individual i would like to have), the model takes the form (2.1). In the present case, the projection on C q can be put in a saddle-point form u -

i∼j p ij G ij = U, -G ij • u ≤ 0 ∀i ∼ j, p ≥ 0, p ij (G ij • u) = 0.
(2.4)

where each p ij accounts for the constraints between i and j, and can be seen as an interaction repulsive force that acts against overlapping. Introducing the matrix B, each line of which expresses a constraint pertaining to a couple in contact, we obtain u

+ B ⋆ p = U, Bu ≤ 0, p ≥ 0, p • Bu = 0.
(2.5)

It formally takes the form of a unilateral Darcy problem (see (2.8) below, which corresponds to a standard Darcy problem in a PDE setting). The fact that B ⋆ can be seen as a discrete gradient is clear from the one dimensional case, with N persons in a row (like in Fig. 2a). Applied to a pressure vector, it computes the (opposite of the) resultant force on i, that is p i,i+1 -p i-1,i . We shall see in Section 4.2 that the interpretation of B ⋆ as a gradient is more tedious in the twodimensional setting, and that the richness of the model actually lies in this very discrepancy between B ⋆ and a proper discrete gradient operator, like the ones that are obtained by consistent space discretization.

2.2. Macroscopic model. The very same principles can be applied at the macroscopic level, as initially proposed in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF]. The population is described by a measure density ρ supported in a domain Ω, that aims at being transported according to a given velocity field U (desired velocity), yet subject to remain in a set of feasible densities

K = ρ ∈ L 1 (Ω) , 0 ≤ ρ ≤ 1 a.e. , Ω ρ(x) dx = 1 . (2.6)
A velocity is considered feasible if it does not lead to a violation of the constraint above, therefore it consists in velocities that do not overconcentrate the density in zones that are already saturated. Unformally said, it requires ∇ • u ≥ 0 wherever ρ = 1. This can be formalized in a dual way :

C ρ = v ∈ L 2 ρ (Ω) , Ω v • ∇p ≤ 0 ∀p ∈ H 1 ρ , p ≥ 0 a.e.
where the set H 1 ρ of pressure test functions is defined as

H 1 ρ = p ∈ H 1 (Ω) , p(1 -ρ) = 0 a.e. .
The macroscopic problem (macroscopic counterpart of (2.1)) simply writes

∂ρ ∂t + ∇ • (ρu) = 0 u = P Cq U, (2.7)
where the continuity equation is meant in a weak sense (see e.g. [START_REF] Piccoli | Time-evolving measures and macroscopic modeling of pedestrian flow[END_REF]), and the projection on C q corresponds to the L 2 norm. The formal similarity with the microscopic approach is underlined by the saddle-point formulation of the projection problem. Let us define the essential saturated zone as the largest open set ω ⊂ Ω such that ρ(x) = 1 for a.e. x ∈ ω. The projection problem takes the form:

Find (u, p) ∈ L 2 (ω) × H 1 0 (ω) such that u + ∇p = U in ω, -∇ • u ≤ 0 in ω, p ≥ 0 in ω, ω u • ∇p = 0, (2.8)
that is the macroscopic counterpart of system (2.5). It takes the form of a Darcy problem, that is commonly used to describe the flow of an incompressible fluid in a porous medium (see e.g. [START_REF] Bear | Dynamics of Fluids in Porous Media[END_REF]). The latter problem provides the effective velocity in ω only. The overall velocity u is obtained by extending u |ω with U (the actual velocity is obviously the desired one outside the saturated zone).

Proposition 2.1. Problem (2.8) has a unique solution (u, p) ∈ L 2 (ω) 2 × H 1 0 (ω).
Proof. This is a straightforward consequence on the fact that the operator

B : v ∈ L 2 (ω) 2 -→ -∇ • u ∈ H -1 (ω) , Bv, p = v • ∇p,
is surjective, thanks to Poincaré inequality in H 1 0 (Ω), so that B ⋆ is one-to-one and has closed range. As a consequence, the polar cone to the set of feasible velocities C q can be written

C • q = w ∈ L 2 (ω) 2 , Ω w • v ≤ 0 ∀v ∈ C q = ∇p , p ∈ H 1 0 (ω) , p ≥ 0 ,
hence the existence and uniqueness of a saddle-point (u, p).

Whether this macroscopic model could be properly obtained through a micromacro limit from the first will be addressed, and in some way answered, in a negative way in Section 4.

Well posedness of the time evolution problems

3.1. Microscopic model. The evolution problem (2.1) consists in finding a con-

tinuous path t → q(t) ∈ R 2N dq dt + ∂I K (q) ∋ U (q) a.e. in [0, T ] , q(0) = q 0 , (3.1)
where I K (q) = 0 as soon as q ∈ K (equal to +∞ otherwise), and ∂I K is the Fréchet subdifferential, defined for q ∈ K by

∂I K (q) = {w , I K (q) + h • w ≤ I K (q + h) + o(h) ∀h} = {w , h • w ≤ o(h) ∀h s.t. q + h ∈ K} . (3.2) Proposition 3.1.
Let the set of feasible configuration K be defined by (2.2), and let the desired velocity q → U (q) be defined as a Lipschtitz function of q. Equation (3.1) admits a unique continuous solution over [0, T ].

Proof. The proof mainly relies on a time-discretization scheme, inspired by the socalled catchin-up algorithm initially introduced in [START_REF] Moreau | Evolution problem associated with a moving convex set in a Hilbert space[END_REF] to build discrete solutions of sweeping processes. Let τ > 0 be a time step, the sequence (q n τ ) is built according to qn+1 = q n + τ U (q n ),

q n+1 = P K qn+1 , (3.3)
Although the set K of feasible configurations is not convex, it can be shown to be prox-regular, which essentially means that projection on it is well-defined in a neighborhood of its boundary (see [START_REF] Maury | A discrete contact model for crowd motion[END_REF][START_REF] Maury | Handling Congestion in Crowd Motion Modeling[END_REF] and the next section). Thus, for τ sufficiently small, the projection of qn+1 on K is well defined. The collection of discrete solutions can be interpolated, which forms a piecewise affine path t -→ q τ . This approximate solution uniformly converges to a solution of (3.1) (see [START_REF] Maury | A discrete contact model for crowd motion[END_REF] for details).

Macroscopic model. The macroscopic problem (2.7) has the general form

∂ρ ∂t + ∇ • (ρu(ρ)) = 0
where u(ρ) is a velocity field, the values of which depend on the density ρ overall the domain, and not only at the considered point. The mapping ρ → u(ρ) is far from being smooth, and the native regularity of the velocity itself is simply L 2 , so that standard theory is not applicable. An alternative approach was proposed in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF][START_REF] Maury | Handling Congestion in Crowd Motion Modeling[END_REF], based on optimal transportation. It consists in extending Moreau's approach, that was initially dedicated to Hilbert spaces, to measures. The scheme reads as follows:

ρn+1 = (Id + τ U ) ρ n , ρ n+1 = P K ρn+1 , (3.4)
where the projection is meant in the square Wasserstein sense, i.e. according to the distance defined by

W 2 (µ, ν) 2 = arg min Λµ,ν Ω |T (x) -y| 2 dµ(x),
where Λ µ,ν is the set of all those measurable maps that push µ forward to ν. The minimum above is attained whenever µ does not charge zero-measure sets (which holds in the situation we consider here), see [START_REF] Villani | Topics in Optimal Transportation[END_REF] or [START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF].

Proposition 3.2.
Let Ω be a bounded convex domain, ρ 0 ∈ K an initial density, and let U be a Lipschitz velocity field, such that U points inward the domain on ∂Ω. We consider t ∈ [0, T ] -→ ρ τ (t) ∈ K (defined by (2.6)) the path obtained by geodesic interpolation between the discrete densities obtained by application of the scheme (3.4), and by u τ the associated velocity. Then the couple (ρ τ , u τ ) converges to a solution (ρ, u) of Problem (2.7).

Proof. We refer to [START_REF] Maury | Handling Congestion in Crowd Motion Modeling[END_REF] for a detailed proof of this convergence. Let us simply address here a particular issue raised by the projection step (second line of (3.4)), that enlight the fact that this macroscopic model is not a straight extension of the microscopic model. Existence of a minimizer comes from standard compactness arguments. As for uniqueness, K is convex in the geodesical (or Mc Cann) sense: for any two measures ρ 0 and ρ 1 in K, and T an optimal map from ρ 0 to ρ 1 , then

ρ t = ((1 -t)Id + tT ) ♯ ρ 0 ∈ K ∀t ∈ [0, 1]. (3.5)
Yet, as pointed out in [START_REF] Ambrosio | Gradient Flows in Metric Spaces and in the Space of Probability Measures[END_REF], the functional ρ → W 2 (ρ, ρ) 2 is not geodesically convex. Uniqueness can nevertheless be obtained by using the fact that the latter functional is convex along generalized geodesics. In the present situation, such geodesics are defined as follows: consider ρ 0 and ρ 1 in K, ρ the measure that is to be projected, and T 0 and T 1 optimal transport maps from ρ to ρ 0 and ρ 1 , respectively. The generalized geodesic is obtained by interpolating the transport maps, i.e.

ρ s = ((1 -s)T 0 + sT 1 ) ♯ ρ.
It can be checked that K is still convex with respect to this new definition of geodesics (see [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF]), and uniqueness follows.

Micro-macro discrepancies

The microscopic and the macroscopic models express the very same principle at different levels of description: the actual velocity is the closest (in a least square sense) to the desired velocity, among all feasible velocities. They futhermore present formal analogies, up to the saddle-point formulations of the minimization problem that defines the instantaneous velocities (Systems (2.5) and (2.8)). Yet, those analogies covers deep discrepancies that are described in this section.

4.1. Geometrical properties of the Wasserstein space. This section addresses some micro-macro issues linked with the geometrical properties of the sets of feasible configurations: K defined by (2.2) for the microscopic situation, and its macroscopic counterpart K defined by (2.6). Although they obviously play and equivalent role, they have very distinct properties from a geometrical standpoint. At the microscopic level, K is the set of hard discs configurations without overlapping. It is obviously non convex, and the euclidean projection on this closed set is not uniquely defined outside a thin neighborhood. The notion of ηprox regularity quantifies this property (see e.g. [START_REF] Thibault | Sweeping process with regular and nonregular sets[END_REF]): K is said to be uniformly η-prox regular whenever ∀q ∈ ∂K , ∀v ∈ ∂I k (q) with |v| = 1 , q = P K (q + ηv),

where ∂I K is defined by (3.2). It can be proven (see [START_REF] Maury | A discrete contact model for crowd motion[END_REF]) that K (defined by (2.2)) is indeed η-prox regular, with a constant η that goes to 0 when N → +∞ together with r → 0. This degeneracy suggests a singular behavior in the micro-macro limit. Indeed, the macroscopic set of feasible density K does not present the characteristic of a prox-regular set. Let us first remark that it has empty interior (for the Wasserstein metric), so that any density ρ that never saturates the constraints belongs to ∂ K, while it is obviously the projection on K of itself only. Even when one considers densities that saturate the constraint on their support, uniform η-prox regularity is immediately ruled out. Consider for example a density that is the indicator function of a collection of non-overlapping discs of common radius 1/ √ N :

ρ N = 1 π N i=1 ½ B(xi,1/ √ N ) .
Among all those measures that projects on K at ρ N , the most remote is a sum of Dirac masses supported by the collection of disc centers, and the corresponding distance scales like 1/ √ N , which goes to 0 as N goes to +∞. In the Hilbertian setting, the projection on such a set (that rules out any uniform prox-regularity) would not be properly defined, even in a close neighborhood of the set. Yet, as mentioned in the proof of Proposition 3.2, the distance from a measure ρ to K is always attained at a unique point in K, no matter how far ρ is from K. The apparent contradiction between the properties of P K and P K are actually due to the fact that the Wasserstein setting corresponds to a Lagrangian description of particles up to permutations. It leads to a puzzling property, that we will illustrate by a simple example: in some cases, considering a density ρ ∈ K, and a vector field that belong to the subdifferential (in the wasserstein sense) of I K (such a velocity field can be seen as directed along the "outward normal direction" of K), the curve ρ t = (Id + tv) ♯ ρ starts indeed moving straight away from K, then turns back toward ρ, and re-enters K at ρ. Let us make it straight that (Id + tv) ♯ ρ does not represent a solution to the transport equation with the Eulerian velocity field v and initial density ρ, but rather the one-shot push-forward (image measure) The example is the following: in the one-dimensional setting, one considers two individuals identified to rigid segments of length one, initially in contact, we consider a velocity field (1, -1) that belongs to ∂K (see (5.1) below), and we consider the microscopic path q t = q + tv. For t ∈ [0, 1/2[, the projection of q t on K is q 0 = q. For t = 1/2 the two segments coincide, and the projection is not unique (the distance is attained for q and for the configuration obtained by interchanging 1 and 2, that are distinct since the description is fully Lagrangian). For larger t, the projection is well defined, but it is not the starting point q. The right-hand side of the figure represents the Wasserstein version of this situation, with the very same velocity field. Segments are represented in the same way, since "particles" are indistinguishable. For t ∈ [0, 1], the projection of ρ t = (Id + tv) ♯ ρ is well-defined, with a unique projection that is the initial configuration. The path ρ t , that follows a velocity field that point outward of K, actually moves away from ρ in a first phase, but then returns back to its very initial position, unlike in the microscopic situation.

B. Maury

Discrete counterparts of the differential operators.

We describe here the discrepancies between B, B ⋆ , and BB ⋆ , discrete analogs of gradient, divergence, and Laplacian operators, respectively, and their macroscopic counterparts.

Discrete divergence. As already mentioned, the matrix B that expresses the non overlapping constraints (see Equation (2.5)) is the discrete counterpart of the opposite of the divergence operator in the macroscopic Darcy equation (2.8). Let us denote by N c the number of contacts. Each row of B ∈ M Nc,2N (R) corresponds to a contact between two discs i and j:

(0, . . . , 0, e ij , 0, . . . , 0, -e ij , 0, . . . , 0) ∈ R 2N , where e ij = (q j -q i )/|q j -q i | is the unit vector between centers. In the onedimensional setting, Bu ≤ 0 prevents inter-center distances to decrease, thus it is a straight counterpart of divergence non-negativity. In the two-dimensional setting, the constraint highly depends on the disc arrangement. In the cartesian case a b c d e (Fig. 2b), Bu ≤ 0 imposes the same kind of constraint along x and y directions, independently, which is much stronger that imposing a local non-negative divergence like in the macroscopic situation. The constraints are even stronger for the triangular lattice (Fig. 2c): the monotone character of the velocity is imposed in three different directions. In the general case (unstructured cluster like in (Fig. 2d)), the contraints are imposed along directions that depend on the local arrangement on discs. Those examples show that the discrete constraint Bu ≤ 0 is much stronger than its macroscopic counterpart -∇ • u ≤ 0.

Discrete gradient. Such a discrepancy can be formulated in a dual way by considering the matrix B ⋆ defined in Section 2, whose role is to transfer Lagrange multipliers (interaction repulsive forces) onto velocity corrections (see Equation (2.5)). This matrix plays the role of the gradient in the macroscopic Darcy equation (2.8). This analogy is fully relevant in the one-dimensional setting. In particular, if one considers a cluster of aligned discs (see Fig. 2a), a uniform pressure fields induces no effect on the discs (except at both ends of the row). This remains true in the two-dimensional setting in very particular, structured situations, like the ones represented in Fig. 2b and Fig. 2c (cartesian and triangular lattice), but this property is lost for clusters with no symmetries (like in Fig. 2d). In the particular configuration represented in Fig. 2e, a uniform pressure field induces a centrifugal force field. In case all discs tend to reach the center of the circle, this allows for a static state to exist, with a exact balance between the desired velocity field and -B ⋆ p. This example can be interpreted as a static jam, and we shall see than such jams have no equivalent in the macroscopic setting. The fact that discs generically have more than 4 neighbors in congested situations induces a mathematical degeneracy that has deep consequences from the modeling standpoint. Consider a large population of N discs arranged according to the triangular lattice (Fig. 2c). The number of primal degrees of freedom (centers of discs) is 2N , whereas the number N c of contacts is of the order 3N , so that B ⋆ = M 2N,Nc (R) is highly singular: the pressure is not defined in a unique way. Let us note though that, in spite of this degeneracy in terms of pressure uniqueness, the solution set is bounded, as expressed by Proposition 4.2. The core of this property relies on the following lemma, which asserts that non trivial degenerate pressure fields (i.e. such that B ⋆ p = 0) always contain pressures with opposite signs. Proof. Let p be such that B ⋆ p = 0, with p ≥ 0. Consider a connected sub-cluster of q (a connected maximal subset of discs in contact), like the one represented in Fig. 2d. We denote by J the corresponding set of indices. Consider now the convex-hull of centers (q i ) i∈J , and pick one q i that is an extremal point of the convex hull. By Hahn Banach's Theorem, q i can be separated from the convex hull of the remaing centers, by an hyperplane normal to some vector w. Now write the force balance for q i along direction w:

j∼i p ij e ji • w = 0.
By construction e ji •w > 0 for all j in contact with i, so that all pressures involving i are equal to zero. Applying the same approach recursively makes it possible to eliminate extremal points one after the other, thus the pressure is identically 0.

Proposition 4.2. Let q ∈ K and U ∈ R 2N be given, u = P Cq U . Then the solution set for pressure, i.e. the set Λ of all those pressure fields p ∈ R Nc + such that (u, p) is a solution to (2.5), is bounded.

Proof. The solution set writes

Λ = p ∈ R Nc , B ⋆ p = U -u , p ≥ 0 .
If a sequence (p n ) of pressures in Λ goes to infinity, then, up to a subsequence,

p n / |p n | converges to p ∈ ker B ⋆ ∩ R Nc + , that is necessarily 0.
Discrete Laplacian. The previous remarks explain a difference in the behaviour of the two models, that is highly significant if one considers evacuation processes.

In the macroscopic setting, the congested zone upstream the exit has the typical shape represented in Fig. 4 (left). The congested zone is denoted by ω, and the desired velocity field points toward the exit, and verifies ∇ • U < 0. The pressure field is then solution of a Poisson problem

-∆p = -∇ • U > 0, (4.2)
with homogeneous Dirichlet boundary conditions on the upper free boundary, and on the exit, and homogeneous neuman condition on the wall (no normal velocity). By the maximum principle, the pressure is nonnegative in ω, so that

u • n = U • n -∂p/∂n ≥ U • n,
where n is the outward normal on the exit boundary. As a consequence, the crowd exits the room faster than it would if there were non congestion. This feature contradicts experimental evidence, in particular the so-called capacity drop phenomenon: beyond a certain threshold, congestion tends to decrease the evacuation speed (see e.g. [START_REF] Cepolina | A Pedestrian Movement Model That Takes into Account the Capacity Drop Phenomenon in the Motion of Crowd, Cellular Automata[END_REF]). This model could be corrected by imposing a flux reduction in case of upstream congestion, in the spirit of what is proposed in [START_REF] Andreianov | Riemann problems with non-local point constraints and capacity drop[END_REF], to force the capacity drop phenomenon, but such a ingredient does not really model, neither explains, the considered phenomenon. On the contrary, the microscopic model natively reproduces this behavior, because of the aforementioned nonstandard character of the underlying discrete operators. Indeed, there exist static configurations with all individuals pointing toward a door that is significantly wider than the widest individual, where the desired velocity is exactly balanced by the action of a positive pressure field. Such a static jam is represented in Fig. 4 (right). In such a situation, the pressure field p solves (2.5) with u = 0. The desired velocity field is concentrating, i.e. it tends to reduce the distances, thus by definition BU > 0. The pressure is then solution to a discrete Poisson-like problem, discrete counterpart of (4.2):

BB ⋆ p = BU, (4.3)
with a positive right hand side. In some structured situations, BB ⋆ is a standard discrete Laplacian. For a 1d row of individuals for example (Fig. 2a), it is the one dimensional discrete Laplace operator with Dirichlet boundary conditions at the ends. The Dirichlet character of boundary conditions expresses the fact that individuals at the end of the row are free to move away from the cluster (i.e. the outside pressure is zero). In the cartesian case (Fig. 2b), this matrix encodes two 1d Laplace operator (along x and y axes). Note the difference with the macroscopic situation: the discrete pressure field does not correspond to a scalar pressure fields, but to a collection of two scalar pressure fields, one for the horizontal direction, one for the vertical one. When the configuration is unstructured, like in (Fig. 2d), the situation is more complex. The operator BB ⋆ is defined on the networks of contact points, that is dual to the network of disc centers (see Fig. 3 (right)). Considering a pressure field p = (p kℓ ), where (k, ℓ) runs over active contacts, the vector BB ⋆ p is defined on the dual network, and the value that corresponds to the contact between i and j is

(k,ℓ)∼(i,j) p kℓ G ij • G kℓ ,
where (k, ℓ) ∼ (i, j) means that the two contacts share a common disc. This matrix shares some properties with matrices associated with a resistive network, or matrices resulting from the space discretization of the Laplace operator by a Finite Difference / Finite Element method, in particular it is symmetric and its diagonal elements are positive. Yet, some extradiagonal coefficients are positive, like in the situation represented in Fig. 5. The element of BB ⋆ that encodes the interaction between (i, j) and (i, k) is e ij • e ik > 0. The resistive networks associated with this matrix therefore admits negative resistances, so that the maximum principle does no longer hold, ruling out the arguments that we used in the macroscopic situation, to assert that individuals escape faster than they would if they were alone. Indeed, the microscopic model reproduces static jams, like the one represented in Fig. 4. In this situation, all individuals tend to reach the exit, but positive pressures create upstream correction terms that exactly balance the desired velocities. One may qualify this phenomenon as some sort of inverse surface tension: when the border of the cluster presents a locally concave shape, a positive pressure field between discs induces a force that points inward the cluster.

Gradient flow structure of the evolution process

Both models may be formulated in a gradient flow framework, in case of evacuation processes. As for the microscopic one, it is natural to define the desired velocity of an individual centered at q i as -∇D(q i ), where D(•) is the geodesic distance to the exit. The latter quantity can be considered as an individual dissatisfaction.

The global dissatisfaction can then be defined as the sum of individual ones, e.g.

Φ(q) = N i=1 D(q i ).
Since the set K is pro-regular (i.e. it verifies (4.1) for some η > 0, and the subdifferential of I K is well defined by (3.2)), it can be checked that the evolution equation (2.1) can be written (see [START_REF] Maury | A discrete contact model for crowd motion[END_REF])

dq dt = -∇Φ(q) -∂I K (q) = -∂ (Φ + I K ) (q),
that is a gradient flow associated to the dissatisfaction functional Φ + I K .

In the macroscopic setting, the Wasserstein framework makes it possible to exhibit a similar structure based on a macroscopic dissatisfaction function canonically defined as

Φ(ρ) = Ω D(x)ρ(x) dx.
Its subdifferential is defined as the set of all those fields w ∈ L 2 ρ (Ω) such that

Φ(ρ) + Ω w • v ρ(x) dx ≤ Φ(T ♯ ρ) + o( T -Id ), (5.1)
for all transport maps T . As detailed in [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF], the evolution takes the expected form of a gradient flow in the Wasserstein space: Again, the formal analogy hides deep differences in terms of behavior. In case of an evacuation like in Fig. 4 (left), the functional Φ + I K is geodesically convex, and there is no density distribution such that 0 ∈ -∂ Φ + I K , i.e. no equilibrium point that would correspond to a static jam. On the contrary, at the microscopic level, a situation like in Fig. 4 (right) correspond to a critical point of the nonsmooth functional Φ + I K . The latter dissatisfaction functional is not convex, and it admits local minima that correspond to static jams. Because of the non smooth character of the functional, stability of such jams cannot be asserted by standard tools, but the generalized Hessian that is introduced in [START_REF] Faure | Crowd motion from the granular standpoint[END_REF] makes it possible to verify that some of those jams that spontaneously appear actually correspond to local minima, and are therefore stable.

Conclusion

The previous considerations show that some features of the microscopic model are not captured by the macroscopic one. This discrepancy between the two models is not a contradiction: the macroscopic model has not been obtained as a rigorous limit of the microscopic one, but rather by simply expressing the same principles at the macroscopic level : the actual velocity field is the closest to the desired velocity field in the set of feasible fields. At the microscopic level, feasibility calls for monotonicity constraints in various directions (e.g. 3 directions for the triangular lattice, as detailed in Section 4.2), whereas, at the microscopic level, the divergence free constraint is scalar. The constraints are much stronger for a hard-sphere collection than for a continuous, deformable medium, solely based on local density. A continuous model fully respectful of microscopic features would necessarily rely on extra quantities reflecting the local arrangement of grains, as illustrated by the several discrepancies listed in the previous section (we also refer to [START_REF] Maury | Handling Congestion in Crowd Motion Modeling[END_REF], Section 5, for a further description of the obstacles that prevent from elaborating a rigorous macroscopic limit of a hard sphere system in highly packed situations). Up to our knowledge, there does not exist any micro-macro framework that allows to recover at the macroscopic scale the features of the microscopic hard-sphere model, in highly congested situations.

Figure 1 .

 1 Figure 1. Crossing of entities, Lagrangian and Wasserstein viewpoints

Figure 2 .

 2 Figure 2. Typical disc arrangments

Lemma 4 . 1 .Figure 3 .

 413 Figure 3. Primal and dual networks

  ω

Figure 4 .Figure 5 .

 45 Figure 4. Macroscopic and microscopic jams

The author is partially supported by ANR Project Isotace (ANR-12-MONU-0013).


Alexander Gepperth

Cem Karaoguz

A bio-inspired incremental learning architecture for applied perceptual problems

We present a biologically inspired architecture for incremental learning that remains resource-efficient even in the face of very high data dimensionalities (>1000) that are typically associated with perceptual problems. In particular, we investigate how a new perceptual (object) class can be added to a trained architecture without retraining, while avoiding the wellknown catastrophic forgetting effects typically associated with such scenarios. At the heart of the presented architecture lies a generative description of the perceptual space by a self-organized approach which at the same time approximates the neighbourhood relations in this space on a two-dimensional plane. This approximation, which closely imitates the topographic organization of the visual cortex, allows an efficient local update rule for incremental learning even in the face of very high dimensionalities, which we demonstrate by tests on the well-known MNIST benchmark. We complement the model by adding a biologically plausible short-term memory system, allowing it to retain excellent classification accuracy even under incremental learning in progress. The short-term memory is additionally used to reinforce new data statistics by replaying previously stored samples during dedicated "sleep" phases.

INTRODUCTION

This contribution addresses the issue of incremental learning for applied robotic scenarios. In particular, we target perceptual learning scenarios where data dimensionalities are typically higher (>1000) than current incremental learning algorithms are able to handle. Incremental learning itself is a notoriously ill-defined term, referring to methods that relax some of the classical assumption of machine learning, namely the avail-Talence, France alexander.gepperth@ensta.fr, cem.karaoguz@ensta.fr ability of training and test data obtained from sampling the same unknown distribution. Instead, incremental learning algorithms may be trained on a set of training data, and then retrained with data sampled from a (more or less) different distribution, potentially multiple times. Incremental learning may also refer to learning algorithms that receive their training samples one by one without knowing their number in advance, instead of processing all examples at the same time, which effectively amounts to changing data distributions as well.

In this contribution, we investigate an even stronger form of incremental learning that often occurs in robotic perception problems: the addition of new perceptual classes (e.g., object classes) to a trained classifier. This can happen, for example, by showing a new object to a robot, or a whole set of new objects. This procedure is depicted in Figure 1.

Problem setting and goals of the study

In perceptual problems, the input dimension I is often quite high, and it is not uncommon for the product IO to exceed 10000 where O denotes output dimension- ality. Furthermore, not only execution but also training time is strongly constrained in a robotics application. Imagine a user teaching a new visual percept to a robot; in this setting, training data must be processed sufficiently quickly to allow interaction, apart from the fact that execution of a trained model must be efficient. Furthermore, the amount of available memory is often severely limited as well if an algorithm is to run on a real robot, potentially on an embedded computer. In addition, in a human-robot interaction scenario, training samples are usually provided sequentially and must be reacted upon immediately, and there is no knowing how many other samples a human might yet provide. Lastly, perceptual problems in robotics are often very challenging ones, so any approach needs to be benchmarked on problems of sufficient complexity.

Therefore, the goal of this study is to present a classification architecture capable of being trained with new classes incrementally as depicted in Figure 1, while being compatible with the following constraints:

• training samples are received one by one

• total number of training samples is unknown

• data dimensionality is high

• execution and training time is limited

• memory is limited

• newly added information should be taken into account immediately

• should work well on challenging and representative real-world tasks

Related work and taxonomy of incremental learning approaches

Incremental learning comes in various forms in the literature, and the use of the term is not always consis-tent: first of all, there are approaches while are capable of what we here denote online learning, meaning that they take their training samples one by one without knowing their number in advance. Most notably, this is achieved by multilayer perceptrons (MLPs) but there are also extensions of the support vector machine (SVM) model [START_REF] Bordes | The huller : a simple and efficient online svm[END_REF][START_REF] Syed | Incremental learning with support vector machines[END_REF] that have this capacity. Online learning is necessary for what we consider to be incremental learning (see list of requirements in the previous section) but not sufficient, as we require as well the capacity to adapt to limited changes (often termed "concept drift" and/or "concept shift" [START_REF] Kulkarni | Incremental learning from unbalanced data with concept class, concept drift and missing features: a review[END_REF]) in data statistics while modifying existing models as little as possible.

In the neural network literature, it is well-known that such scenarios lead to what is termed "catastrophic forgetting" [START_REF] Goodfellow | An empirical investigation of catastrophic forgetting in gradient-based neural networks[END_REF], the sudden overall deterioration of model performance. To perform true incremental learning in the aforementioned sense, most approaches perform a local partitioning of the input space and train a separate classification/regression model for each partition [START_REF] Vijayakumar | Locally weighted projection regression: An o(n) algorithm for incremental real time learning in high-dimensional spaces[END_REF][START_REF] Nguyen-Tuong | Local gaussian processes regression for real-time model-based robot control[END_REF][START_REF] Sigaud | On-line regression algorithms for learning mechanical models of robots: A survey[END_REF][START_REF] Butz | Computational complexity of the xcs classifier system[END_REF][START_REF] Cederborg | Incremental local online gaussian mixture regression for imitation learning of multiple tasks[END_REF]. The manner of performing this partitioning is very diverse, ranging from kd-trees [START_REF] Cederborg | Incremental local online gaussian mixture regression for imitation learning of multiple tasks[END_REF] to genetic algorithms [START_REF] Butz | Computational complexity of the xcs classifier system[END_REF] and adaptive Gaussian receptive fields [START_REF] Vijayakumar | Locally weighted projection regression: An o(n) algorithm for incremental real time learning in high-dimensional spaces[END_REF]. Equally, the choice of local models varies between linear models [START_REF] Vijayakumar | Locally weighted projection regression: An o(n) algorithm for incremental real time learning in high-dimensional spaces[END_REF], Gaussian mixture regression [START_REF] Cederborg | Incremental local online gaussian mixture regression for imitation learning of multiple tasks[END_REF] or Gaussian Processes [START_REF] Nguyen-Tuong | Local gaussian processes regression for real-time model-based robot control[END_REF]. Since this article is concerned with high-dimensional perceptual problems, it can be stated for all cited approaches that it is really the partitioning of the input space that is costly in terms of memory. Most notably, covariance matrices used in [START_REF] Vijayakumar | Locally weighted projection regression: An o(n) algorithm for incremental real time learning in high-dimensional spaces[END_REF] are quadratic in the number of input dimensions which makes their use prohibitive in our problem setting.

Bio-inspired approach to incremental learning

As biological incremental learning has reached a high degree of perfection, we explicitly investigated the biological literature for hints as how to this might be achieved. Basing ourselves on observations from the basic sensory cortices, we noted that sensory representations seem to be prototype-based, where prototypesensitive neurons are topologically arranged by similarity [START_REF] Tanaka | Inferotemporal cortex and object vision[END_REF][START_REF] David A Leopold | Norm-based face encoding by single neurons in the monkey inferotemporal cortex[END_REF][START_REF] David A Ross | Not just the norm: Exemplar-based models also predict face aftereffects[END_REF][START_REF] Erickson | Clustering of perirhinal neurons with similar properties following visual experience in adult monkeys[END_REF]. Learning seems to act on these representations in a task-specific way, where more prototypes are allocated to sensory regions where finer discrimination is necessary [START_REF] Daniel B Polley | Perceptual learning directs auditory cortical map reorganization through top-down influences[END_REF], i.e., where more errors occur during learning. Learning is conceivably enhanced through acetylcholine release in case of task failures [START_REF] Norman | The nucleus basalis and memory codes: Auditory cortical plasticity and the induction of specific, associative behavioral memory[END_REF][START_REF] Michael | The role of acetylcholine in learning and memory[END_REF], leading to higher "prototype density" in difficult regions of the sensory space. In particular, learning seems to respect and even generate topological layout of prototypes by changing only a small subset of neural selectivities [START_REF] Rolls | The effect of learning on the face selective responses of neurons in the cortex in the superior temporal sulcus of the monkey[END_REF] at each learning event, namely around those neurons that best matched the presented stimulus [START_REF] Erickson | Clustering of perirhinal neurons with similar properties following visual experience in adult monkeys[END_REF].

We model these findings by using a self-organizing map (SOM) learning to shape the feature preferences of hidden layer neurons in our architecture. SOM is a prototype-based algorithm in the sense that the hidden layer weight vectors "live" in the space of inputs and aim to approximate the probability density in that space. We model the global, task-related error signal by the current classification error that activates SOM learning in case of mismatch or ambiguity. As SOM learning attributes more prototypes to regions where many learning events occur, this will ensure that prototype density increases in difficult regions of the input space. Furthermore, SOM adaptation is stably selfterminating since no more learning will occur when no more errors are made. Inversely, when error rates increase due to the presentation of new input statistics, the hidden layer representation will become plastic until error rates subside again, when a sufficient re-adaptation has been achieved. Lastly, SOM produces a topologically organized representation of the input space, which was the reason to formulate the model in the first place, and modifies weights only locally in case of learning (as observed in biology). It is above all this localized adaption property that makes us choose the SOM algorithm, as this is an essential prerequisite for incremental learning as will be detailed later in this section. On may speculate that the observed topological organization of selectivities in biological neurons serve just for that: to facilitate incremental learning.

Concretely, we created a three-layer neural network, see Figure 2, that learns a set of plastic, topologically organized prototypes in its hidden layer, implementing the approximation scheme for neighbourhood relations. Consequently, a learning scheme that modifies the best-matching prototype and its (2D) neighbours will be strictly local in input space as well, see Fig. 3. A read-out mechanism between hidden and output layer maps local input space regions (i.e., sets of prototypes) to class memberships using simple linear regression learning. Update of hidden layer prototypes is strictly controlled and occurs only when output layer activity is incorrect or ambiguous, on the grounds that prototypes that allow satisfactory classification need not be adapted. Inversely, the mapping from hidden to output layer is adapted only when there are sufficiently similar prototypes in the hidden layer, since excessive dissimilarity means that unknown data is being fed to the network, and appropriate prototypes must be formed before meaningful readout can occur.

On the computational side, all major approaches to incremental learning, (see Section 1.2), perform a partitioning of the input space and learn independent models in each partition. In this way, learning is always localized in the input space, in the sense that a change of statistics in one part of that space will not affect learning in other, distant parts. It is this property that avoids "crosstalk" and therefore catastrophic forgetting that is observed in many connectionist models [START_REF] Goodfellow | An empirical investigation of catastrophic forgetting in gradient-based neural networks[END_REF]. The presented architecture follows this approach using an approximation scheme that strongly simplifies the definition of local regions which is very costly in spaces of high dimensionality I, for example when using a covariance matrix that will contain I 2 entries. To avoid this, local regions are coarsely approximated by hyperspheres whose centers are defined by support points (or prototypes). The quality of this approximation can be controlled by controlling the overall number of prototypes. As such a prototype-based representation approximates the distribution of data points in input space as a whole, it is a generative model [START_REF] Bishop | Pattern recognition and machine learning[END_REF] as it could be used for sampling purposes. As a local region will usually be defined by more than one prototype, and learning should act on all or at least many prototypes of a region, a way needs to be found to adapt "nearby" prototypes together, amounting to the need of an efficient neighbourhood relation between prototypes. We approximate this by placing prototypes onto a two-dimensional lattice where lattice distance expresses closeness in input space, this property being assured by the learning rule for prototype adaptation. This approximation scheme (which is conveniently and implicitly generated by the SOM learning algorithm), which leads to drastically reduced memory requirements for highdimensional problems, is explained in Figure 3.

We also incorporate the interplay between shortterm and long-term memory in our model. There is a large body of literature investigating the roles of the hippocampal and neocortical areas of the brain in learning. Generally speaking, the hippocampus employs a rapid learning rate with separated representations whereas the neocortex learns slowly, building overlapping representations of the learned task [START_REF] Randall | The division of labor between the neocortex and hippocampus[END_REF]. A well-established model of the interplay between the hippocampus and the neocortex suggests that recent memories are first stored in the hippocampal system and they are played back to the neocortex over time [START_REF] Mcclelland | Why there are complementary learning systems in the hippocampus and neocortex: Insights from the successes and failures of connectionist models of learning and memory[END_REF]. This accommodates the execution of new tasks that have not been recently performed as well as the transfer of new task representations from the hippocampus (short-term memory) to the neocortical areas (long-term memory) through slow synaptic changes. Inspired from these biological findings, we investigated the effects of such a setup on our model.

Summarizing, we have tried to incorporate as many facts about incremental learning in biology as possible while keeping the model as simple and efficient as possible. Our modeling takes place at the architectural level, leaving aside the finer details of neural modeling (rate/spike code, dynamic neuron models etc.).

METHODS

As stated in Section 1.3, we implement our approach to incremental learning as a three-layer neural network architecture depicted in Figure 2. To train the hidden layer H of topologically organized prototypes, we use a learning scheme adapted from the well-known self-organizing map (SOM) algorithm [START_REF] Kohonen | Self-organized formation of topologically correct feature maps[END_REF] which is described in Section 2.1, whereas the readout from hidden to output layer O is performed by linear regression (LR) explained in Section 2.2. A particular point are the modulation influences within the architecture, which control and restrict learning in hidden and output layers, see Section 2.3.

We denote neural activity vector in a 2D representation X by z X ( y,t), and weight matrices for SOM and LR, represented by their line vectors attached to target position y = (a, b), by w SOM y . For reasons of readability, we often skip the dependencies on space and time and include them only where ambiguity would otherwise occur. Thus we write z X instead of z X ( y,t) and w SOM instead of w SOM y (t).

Activity generation and prototype training in the hidden layer

The hidden layer H is not intended to reduce the dimensionality of the input but rather to re-encode it in a way that enables incremental learning, see Figure 3, while preserving information. Therefore, instead of reducing the output of the hidden layer to the best-matching unit (as it is usually done for the SOM model), we calculate the (graded) activations of all hidden layer units for performing the following linear regression. Hidden layer activations z H are meant to measure a sparse similarity between input and the prototype associated to a particular unit and are normalized in the [0, 1] interval. They are obtained by first passing all input-prototype distances zH through a Gaussian function with standard deviation κ, and then applying a transfer function TF(•) that sparsifies these similarities. Here, there is a technical point to be observed: since there is no way of knowing a priori the typical input-prototype distances, κ must be adapted to the data during the learning process, making it a dynamic, timedependent quantity:κ ≡ κ(t).

In order not to trigger the adaptation of too many linear regression weights, the transfer function thresholds these similarities, thus ensuring that only those units whose prototypes are very close to the input are active and take part in linear regression training.

Prototype adaptation is performed using the conventional SOM update step except that it takes into account a control signal λ (t) coming from the output level of the hierarchy, which will be described in Section 2.3. In precise terms, the equations for activity generation in the hidden layer in response to input activity z I are as follows:

zH ( y) = ||w SOM y -z I || (1) zH = g κ zH ( 2
)
z H = TF θ ,p zH (3)
The adaptation steps for the weights and the distance metric κ(t) are as follows:

w SOM y (t + 1) = w SOM y + λ (t)ε SOM g σ (|| y -y * ||)(z I -w SOM y ) (4) κ(t + 1) = (1 -τ) κ(t) + τ max y zH ( y) (5) range meaning ∈ [-1, -0.5] certain, incorrect ∈ [-0.5, -0.1] uncertain, incorrect ∈ [0.1, 0.5] uncertain, correct ∈ [0.5, 1] certain, correct
TF θ ,p (z H ) = TF p (z H ) if TF(z H ) > θ 0 otherwise (6)
where

TF p (z H ) = zH p max y (z H ( y,t)) p-1 (7)

Decision making and learning in the output layer

Generation of output layer activities is performed by a simple linear transformation of thresholded hidden layer activities z H , using the linear regression weights w LR . Learning is subsequently modifying these weights to optimize the mapping of hidden layer activities z H to the target representation z T containing the "true" class of a sample.

z P ( y) = w LR y • z H (8)
w LR y (t + 1) = w LR y (t) + 2ε LR z H z P ( y) -z T ( y) (9)
In contrast to the hidden layer learning rate, the learning rate of linear regression, ε LR remains constant at all times.

Feedback control of learning and detection of novelty

Hidden layer and output layer neurons do not adapt their weights all the time, only when it is deemed appropriate. Here, there are two distinct cases to be distinguished: first, when the hidden layer has low overall activity, maybe because the input belongs to a newly added class, then linear regression should not adapt its weights because hidden layer activity is probably not meaningful and it will impair performance for already represented classes. This is achieved automatically by thresholding hidden layer activities in the transfer function T F θ ,p by θ as specified in Equation( 3), leading to zero activity if activity is too low. In this case, linear regression weights will not be updated as this requires non-zero activities in the hidden layer H. Secondly, hidden layer weights w SOM will only be updated when the current estimate of class membership, i.e. the output layer activities z P , is either uncertain or wrong. To measure uncertainty, we first define an uncertainty measure based on the output layer activities, whose basic idea is that a certain estimate of class membership has a clear activity maximum, so a good measure is just to use the bounded difference between first and second maximum:

u(z P ) = max y z P -max2 y z P (10)
This measure, u(z P ), can be combined with the fact whether the activity maximum of z P is in the right place, i.e., in accordance with ground-truth information z T :

m(z P , z T ) = u(z P ), if arg max y z P = arg max y z T -u(z P ), otherwise (11)
Finally, we obtain the modulation measure λ (t) that decides whether hidden layer weights should be trained, by thresholding the confidence measure m(z P , z T ):

λ (t) = 0, m(z P , z T ) > θ m 1, otherwise (12)
Table 1 gives an overview over the intuitive meaning of various values of the confidence measure m. By thresholding m(z P ) we thus allow hidden layer weights to be trained only when the current class estimate is of less-than-perfect quality, either outright incorrect (for θ m ≤ 0) or correct but of significant uncertainty (for θ m > 0).

Incremental learning with short-term memory

We employ a short-term memory mechanism in the form of a simple queue where novel experiences are stored. The novelty of an experience is determined via the comparison of the supervision signal to the prediction z P of the model which we rename here for convenience to z P-SOM : if the model cannot correctly classify the input (corresponding to m(z P-SOM , z T ) < 0), the feature vector z I and its label vector z T are stored in the short-term memory. The short-term memory has limited capacity, hence a new experience overwrites the oldest entry if the short-term memory is full.

In accordance with biological findings, the shortterm memory structure is used in both prediction and learning phases. In the prediction phase, response of the short-term memory representation S, zS to an input z I is computed based on a distance measure, similar to the hidden layer of the model (see Equation 1):

zS (i) = ||w STM i -z I || (13)
where w STM i is the prototype i stored in the short-term memory. The decision whether the model or the shortterm memory should be used for final prediction z P is done by comparing the activities zH and zS :

z P = z P-STM , if min zS < min zH z P-SOM , otherwise (14)
where z P-STM is the label stored in the short-term memory along with prototype i * = arg min i zS . In this way, the short-term memory shall assist the classification task when the system is exposed to samples from a new class of object (i.e. the incremental learning stage).

Mimicking simulated experiences that occur during sleep [START_REF] Shen | Modeling the spontaneous reactivation of experience-specific hippocampal cell assembles during sleep[END_REF], contents of the short-term memory are additionally replayed to the model in certain intervals M which are referred to as sleep phases, where the whole content of the STM is replayed to the system, thus training hidden and output layer weights. In this way, novel experiences do not directly cause synaptic changes in the system, they are rather exposed to the model over the short-term memory. This may help the model to quickly build representations for new classes and enable the use of different dynamics for learning new classes.

Training and testing schedules and dedicated incremental learning measures

The hidden layer training algorithm requires an initialization phase (until T conv , see Section 2.1) and, after that, needs to be trained on P -1 classes for a while to provide a starting point for incremental learning. The initial global prototype ordering is terminated at T 1 , and hidden layer initialization is performed until T conv without modulation (λ (t) ≡ 1) as modulation information can only be obtained from a converged hidden layer via the subsequent linear regression stage. Non-incremental training is conducted until T inc1 with modulation switched on.

When learning without short-term memory, a first incremental training stage from T inc1 to T inc2 , and a second stage from then on until the end of the experiment at T eval . For the first incremental training stage, the threshold parameter θ m , see Section 2.1, is set to a high value: θ m = θ inc m , whereas it is set to zero for all other training or testing stages: θ m = 0. This is done to "protect" the classes that are not in the input data, as a new class would initially lead to quasirandom, weakly localized hidden layer activity which could quickly unlearn previously learned linear regression models. Similarly, the neighbourhood radius σ (t), see Section 2.1, is multiplied by a factor of 10 during the first incremental stage, and slowly decreased (back) to its long-term stable value σ ∞ during the second incremental stage. This allows it to perform incremental learning, which requires a large σ for connected regions sensitive to the new class to form in the hidden network layer, while keeping high precision which requires that the final neighbourhood radius σ ∞ be very small.

When learning with short-term memory, a single incremental training stage is conducted from T inc1 to T eval , maintaining θ m = 0 and σ (t) = σ ∞ . Temporal sequence of initialization, training and evaluation steps for the setup without short-term memory are given in Figure 4. For the setup with short-term memory, temporal organization of training and test phases is shown in Figure 5 where the exposition of novel examples is embedded in sleep phases (shown as green periods) and done implicitly. While the scheme illustrated in Figure 4 is still viable for the system with short-term memory, blending the incremental learning phase into the whole learning process as in Figure 5 may be more plausible for life-long autonomous learning. All other factors of temporal schedules of training and testing are identical to the setup without short-term memory.

EXPERIMENTS

The experiments in this article were conducted on the well-known MNIST handwritten digit recognition benchmark [START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF], representing a task of moderate difficulty and real-world relevance problem with however a large number [START_REF] Tanaka | Inferotemporal cortex and object vision[END_REF] of categories. We evaluate performance as indicated by Figure 1: training out model in a leave-one-out fashion on all but one class and the adding the remaining class. We also conduct a final training phase with all classes included before evaluation on the independent test set to draw conclusions about the efficacy of each approach and possible combinations. For a classification task with P classes, this gives essentially P experiments. Please see Table 2 for relevant properties of the used data set, and Fig. 6 for representative data samples.

Global parameter settings

We use the following fixed parameters for our system, where the total number of neurons in the hidden layer is n × n: n = 30, ε LR = 0.001, ε 0 = 0.1, σ 0 = 0.3n, T 1 = 50000, T conv = 150000, T inc1 = 800.000, T inc2 = 820.000, T eval = 900.000 ε ∞ = 0.001, σ ∞ = 0.05, p = 20, θ m = 0.0, θ inc m = 0.75, θ = 0.75 and τ = 0.001. SOM, STM and LR weight matrices are initialized to random uniform values between -0.001 and 0.001.

Training examples are always randomly and uniformly drawn from the current training set. The capacity of the short-term memory is set to C = 300, the system is going into a sleep-phase in every 5000 iterations, replaying all examples stored in the STM.

Baseline evaluation

To establish a baseline to which we can compare incremental learning performance, we train the proposed architecture on the MNIST data set as described in Section 2, using all classes, while evaluating generalization performance directly after the initialization phase, see Figure 4, yielding a test error rate of 4.8% for MNIST. This is a baseline performance that is quite competitive when compared to other results on this popular benchmark, especially since we are using a rather small hidden layer size of 30x30 elements. With higher hidden layer sizes, the error rate can be improved to < 2%.

Evaluation of incremental learning performance without short-term memory

Incremental learning experiments on the MNIST data set use the scheme explained in Figure 4. Ten experiments are conducted in total each of which uses one of the 10 classes as the excluded class, labeled Inc-0 through Inc-9. Along with the baseline experiment, Figure 7 shows the evolution of test errors.

As expected, overall error increases as soon as samples from the new class are introduced. The rapid decrease of classification error shows that system learns the newly introduced class quickly. The final error values are usually higher than the values obtained right before the introduction of the new class. However, a slight increase in the overall error should be expected since an added class creates a more difficult classification problem. In addition, we find that the second incremental training step in Figure 1, with all classes included, is unavoidable in almost all cases, a notable exception being the Inc-1 experiment on MNIST shown in Fig- ure 7. This is quite understandable, as a new class in most cases resembles an existing one, thus partly occupies the same volume in input data space and therefore necessarily breaks readout models in this region. This usually concerns just a few classes, as we can see in experiment Inc-0 on MNIST from Figure 7: MNIST classes "5" and "6", which are very similar, are seriously impacted by incremental training of class "0" but none of the other classes. The additional training step with all classes essentially corrects the balance between classes that are too similar, which is why it improves performance strongly.

Qualitative effects and biological analogy of incremental learning

Fig. 8 qualitatively compares the impact incremental learning has on hidden layer prototypes for experiment Inc-0 without short-term memory. We can observe several things: first of all, prototypes of a class are not always homogeneously arranged but can formed several "islands", which is normal since prototypes are not organized by class membership but by Euclidean distance. Furthermore, and more surprisingly, prototypes of a class do not always vote most strongly for that class in the linear regression. This is however just an expression of efficiency, as the linear regression uses all the information at its disposal: if a prototype that encodes the digit "1" can be used for classifying the class "7" as well then it is natural that linear regression will try to do this. And lastly, we can clearly see how class "0" intrudes into what was previously a region covered by prototypes of class "6". When observing the development of classification errors for class "6" in Fig. 7, we find a corresponding increase for class "6" during the first phase of incremental learning, reflecting the destruction of many of is prototypes. However, classification performance recovers, probably also because several of the prototypes for "0" actually help classifying "6", as we can plainly see from the color coding.

These results are meant to illustrate the global manner in which incremental learning in the proposed architecture works in practice. They show that new classes are smoothly and strictly locally embedded into existing structures, and always at positions where there is the greatest similarity to existing prototypes, which ensures, as a side effect, that the "invaded" classes may still profit from the new prototypes as they remain very similar.

Effects of short-term memory

Experiments with the short-term memory assisted model were done similar to the previous ones. The capacity of the short-term memory is set to C = 300, the system is going into a sleep-phase in every 5000 iterations. All other parameters of the model were the same as the previous experiments (see Sec. 3.1). Figure 9 shows the results from experiment Inc-6. The plot on the left shows the frequency of short-term memory usage for prediction throughout the whole experiment. As expected, the system utilizes the short-term memory more often when new experiences are presented: in the beginning and at the start of the incremental learning phase. The impact of the short-term memory in prediction performance becomes apparent if the error plots from experimental setups with and without short-term memory are compared: high peaks of error that occur at the start of the incremental learning phase with the system without short-term memory (Figure 7, middle and right) are suppressed by the utilization of the short-term memory (Figure 9, right).

From the perspective of learning, in some experiments with the system where short-term memory is not utilized, the newly introduced class may cause disruptions in another class that is already learned when the classes share similar features. For example, in experiment Inc-0, the newly added class 0 disrupts already learned class 6 (Figure 7, middle). On the other hand, from the experiments conducted with the model assisted by a short-term memory it can be observed that after the introduction of new samples, test errors per class remain consistent while the model gradually takes over the prediction task from the short-term memory. This suggests that the model successfully learns representations for the new class without breaking the old ones.

DISCUSSION

Significance of contribution

We consider the results of this study very significant in the following sense:

• we present an incremental learning architecture that stays resource-efficient at very high input/output dimensionalities, especially w.r.t. memory consumption, which is not the case for any other incremental learning algorithm we are aware of. For a hidden layer with N 2 elements and input/output dimensionalities of I and P, the memory requirements scale roughly with N 2 (I + P) and are thus linear in the sum of input and output dimensions. This means that memory consumption is in no case limited by high problem dimensions but rather by the amount of hidden layer units one is willing to spend for better performance.

• the architecture is conceptually simple, efficient to execute and highly parallelizable as we could show in a previous publication [START_REF] Gepperth | Biologically inspired incremental learning for high-dimensional spaces[END_REF] • we present an intuitive way to benchmark incremental learning and apply it to our architecture using two challenging perceptual problems, with excellent results regarding incremental learning per-formance, which never degrades overall precision by more than a few percents.

• we present a biologically plausible, efficient and easy-to-implement way of simplifying incremental learning even further using a short-term memory system

Discussion of methodology

The architecture without short-term memory performs incremental learning successfully but in a rather complicated and ad hoc way. In particular, the system has to know (i.e., be told) when incremental learning is going to happen to that certain parameters (notable θ m and σ i n f ) can be set properly for successful incremental learning. While this is a valid approach (as in most applications of incremental learning, e.g., in robotics, it is known when new things must be learned), it is neither very elegant nor very easy-to-use: it would be much more appealing if the system could detect conceptually new samples autonomously, and perform incremental learning on them without the user having to prepare samples beforehand in order to appear exclusively (first incremental training phase). The addition of the shortterm memory does just that! It is a very nice side effect that it is in close analogy to biological memory architectures, which may shed some light on the underlying rea- sons such memory architectures have evolved. From a purely technical point of view, however, the short-term memory solves a concrete need while leaving overall classification performance completely unchanged.

Comparison to related work

It was attempted to compare our approach to existing prior art in the field of incremental learning, for which we chose the LWPR algorithm [START_REF] Vijayakumar | Locally weighted projection regression: An o(n) algorithm for incremental real time learning in high-dimensional spaces[END_REF] following the best practices for choosing parameters as laid down in [START_REF] Vijayakumar | A library for locally weighted projection regression[END_REF]. However, when working on MNIST data having 756 dimensions and 10 classes, we found that • performance was unsatisfactory (approximately 15% overall error) when using a single output variable c expressing the class as c ∈ {1, 2, 3, . . . }

• performance improved when using a binary-coded vector of 10 elements, c, to express class membership. However the number of receptive fields was limited to approximately 15 before freezing the computer (3GHz off-the-shelf PC running Linux, 4GB memory) due to lack of memory. In this case, overall error rates of approximately 7% could be achieved which is well above the error rates achieved by our architecture.

For these reasons, we found that a fair comparison was impossible due to the memory requirements of LWPR when faced with high input and output dimensionalities.

As LWPR was not capable to create as many receptive fields as it considered necessary, it could not achieve the best possible performance. We suppose that, using computers with more memory, a fairer comparison could be conducted which will be the subject of subsequent experimental studies. For the time being, we can conclude that LWPR performs much worse than out model in terms of memory consumption when applied to the MNIST dataset, and that its performance, as far as it could be measured, is inferior. We believe that the problem of memory complexity for large input/output dimensionalities will manifest itself for most or all other incremental learning algorithms listed in Sec. 1.2 as all of them perform a partitioning of the input space, often in a fashion that is comparable to the one LWPR uses.

Conclusion and future work

We have presented a neural architecture for incremental learning that is novel in the sense that it can operate in high-dimensional perceptual spaces. The functionality and feasibility of this architecture has been demonstrated on two real-world perceptual classification benchmarks of sufficient difficulty, showing that very high input/output dimensionalities are no obstacle at all, and that the addition of a new perceptual class to a trained model results in only a very small performance impairment. Furthermore, we included a biologically very plausible version of a short-term memory system, which fulfills very concrete functions: first of all, it allows the system to use a much simpler incremental learning scheme with just a single learning phase. Secondly it permits much quicker "reactions" to new classes, giving correct classifications almost instantly where otherwise thousands if iterations would be needed for model convergence. And lastly, it enables the system to concentrate learning on samples from new classes that are embedded into samples of known classes: as the new samples arrive relatively rarely, model convergence would be very slow without STM, which was the reason in the first place to split incremental learning into two phases, with new examples submitted exclusively during the first incremental phase.

Future work will include a careful study of the dynamic interplay between long-term and short-term memory, aiming at simplifying the architecture and making as many parameters as possible self-adaptive. This is a prerequisite for doing extended benchmarks on more challenging, real-world classification problems on which the presented architecture should perform well in as generic a fashion as possible.

Compliance with ethical standards

This article does not contain any studies with human participants or animals performed by any of the authors. Cem Karaoguz has received a research grant from MBDA Missile Systems. Alexander Gepperth and Cem Karagouz declare that they have no conflict of interest.

Figure 1 .

 1 Figure 1. Incremental learning scenario targeted in this study. A) Initial training leaving out a subset of classes (in this paper always just one class). B) optional: incremental retraining with the "missing" classes (for system without short-term memory) C) Re-training with all classes. D) Evaluation of overall classification performance (for all classes) on an independent test set.

Figure 2 .

 2 Figure 2. Architecture overview

Figure 3 .

 3 Figure 3. Illustration of how incremental learning is made possible through a topologically ordered prototype representation. Due to topological ordering, neighbouring prototypes are almost always situated in nearby regions of input space. Therefore, local updates of prototypes will almost always be local in input space as well, thus effectively enabling efficient incremental learning. This is shown here for a subset of prototypes trained on the MNIST database, the best-matching unit (BMU) for a "5" input being indicated by a small red circle. It is obvious that the local 2D update region, indicated by a larger red circle, is indeed local in the input space. The yellow circle indicates a region where this property does not hold (structural defect) but the reader can convince himself that this occurs but rarely.

Figure 4 .Figure 5 .

 45 Figure 4. Temporal organization of training and testing phases without short-term memory (x axis not to scale).

Figure 6 .

 6 Figure 6. Representative examples from the MNIST data set used in this study.

Figure 7 .

 7 Figure 7. Results of incremental learning experiments Inc-0 through Inc-9 conducted on MNIST. Left: development of overall classification performance for Inc-0 through Inc-9. Middle: development of individual class errors for Inc-0. right: development of individual class errors for Inc-1.

Figure 8 .

 8 Figure 8. Qualitative effect of incremental learning on hidden layer prototypes for experiment Inc-0 without STM. Shown is a visual representation of prototypes for each of the n 2 hidden layer units. Left: prototypes just before incremental learning of class 0 at t = T inc1 . Right: prototypes at t = T inc2 . The hue of each prototype corresponds to the class it votes for most strongly, which is read out from the linear regression weights.

Figure 9 .

 9 Figure 9. Results of incremental learning experiment Inc-6 with short-term memory conducted on MNIST. Left: frequency of short-term memory usage for prediction. Middle: development of individual class errors for Inc-6. Right: test error of the same experiment zoomed in to the incremental phase for better visualization. The noise in the curves after 800K iterations is due to higher sampling rate of test error.

Table 1 . Significance of various value ranges of the confidence measure

 1 m(z P ).where g σ (x) is a zero-mean Gaussian function with standard deviation s and y * denotes the position of the best-matching unit (i.e., the one with the highest similarity-to-input) in H. In accordance with standard SOM training practices, the SOM learning rate and radius, ε SOM and σ , are maintained at ε 0 , σ 0 for t < T 1 and are exponentially decreased afterwards in order to attain their long-term values ε ∞ , σ ∞ at t = T conv . TF represents a monotonous non-linear transfer function, TF : [0, 1] → [0, 1] which we model as follows with the goal of maintaining the BMU value unchanged while non-linearly suppressing smaller values:

Table 2 . Data set used in this study and its prin- cipal properties.

 2