Selective Memory Recursive Least Squares: Uniformly Allocated Approximation Capabilities of RBF Neural Networks in Real-Time Learning

When performing real-time learning tasks, the radial basis function neural network (RBFNN) is expected to make full use of the training samples such that its learning accuracy and generalization capability are guaranteed. Since the approximation capability of the RBFNN is finite, training methods with forgetting mechanisms such as the forgetting factor recursive least squares (FFRLS) and stochastic gradient descent (SGD) methods are widely used to maintain the learning ability of the RBFNN to new knowledge. However, with the forgetting mechanisms, some useful knowledge will get lost simply because they are learned a long time ago, which we refer to as the passive knowledge forgetting phenomenon. To address this problem, this paper proposes a real-time training method named selective memory recursive least squares (SMRLS) in which the feature space of the RBFNN is evenly discretized into a finite number of partitions and a synthesized objective function is developed to replace the original objective function of the ordinary recursive least squares (RLS) method. SMRLS is featured with a memorization mechanism that synthesizes the samples within each partition in real-time into representative samples uniformly distributed over the feature space, and thus overcomes the passive knowledge forgetting phenomenon and improves the generalization capability of the learned knowledge. Compared with the SGD or FFRLS methods, SMRLS achieves improved learning performance (learning speed, accuracy and generalization capability), which is demonstrated by corresponding simulation results.

capability), appropriate training methods are essential. For neural network control, the RBFNN is usually employed as an adaptive controller or observer to interact with the plant online, which leads to the considerable attention of the realtime training methods of RBFNNs. Among different realtime training methods, the gradient descent (GD) based and recursive least squares (RLS) based methods have received extensive attention and been applied to various systems [7]- [10].
As one of the most classical training methods of neural networks, the GD method and its variants have achieved good performance in various offline learning scenarios, where the training sets are known a priori [11]. Different from the offline applications of neural networks, the real-time learning, as a special case of online learning, refers to the scenarios where the samples are collected one by one in chronological order. For instance, in neural network control where the RBFNN should keep sensitivity to parameter changes, the GD based update of the neural network will be carried out whenever a new sample is collected, which can also be regarded as a stochastic gradient descent (SGD) based training method with the batch size of one. Generally, the SGD based weight update law of the RBFNN in neural network control is derived from the Lyapunov stability theory based analysis and thus can guarantee the stability and closed-loop control performance of the system [12]- [14]. By contrast, the learning performance of the RBFNN in this case is very limited because the learned knowledge is closely related to the chronological distribution of the samples input to the RBFNN. Specifically, since the objective function of SGD is squared approximation error of the current sampling time, the RBFNN will mostly allocate its approximation capabilities to the latest samples while the knowledge learned from the past samples will be gradually forgotten, which is named passive knowledge forgetting phenomenon in this paper.
Due to passive knowledge forgetting phenomenon, the RBFNN needs to persistently review the past samples to realize the real-time learning (parameter convergence to the optimal value). From the perspective of adaptive control, this is the reason why the weight convergence of the RBFNN based on SGD requires the persistent excitation (PE) condition [15], [16]. When employing an RBFNN as the function approximator, what we expect is to obtain accurate approximation of the unknown function over a large-scale feature space (i.e. input space) region of the RBFNN rather than over local regions filled with the latest samples. Therefore, the training method with memory of the past samples is required to improve the real-time learning performance of the RBFNN.
As a result of its linearly parameterized form, the realtime training of the RBFNN can also adopt the RLS method. While the SGD method optimizes an objective function of the current squared approximation error, the objective function of the RLS method is the sum of squared approximation errors at each sampling time. Therefore, the RBFNN trained through RLS will allocate its approximation capabilities uniformly to each sample and keep memory to all of the past samples during the real-time learning process. Moreover, it is claimed that the RLS based training method has faster learning speed compared with the SGD based method [17]- [20]. However, as the number of samples gradually increases, the approximation capabilities of the RBFNN will be insufficient to express the excessive knowledge such that the RBFNN will be insensitive to new samples, which is also known as the data saturation phenomenon. With data saturation, the RBFNN will allocate its approximation capabilities mostly to the feature space regions where the samples are dense, resulting in limited learning performance.
To overcome the data saturation phenomenon of RLS, different forgetting mechanisms such as the forgetting factor are adopted [21]- [26]. However, while data saturation is suppressed by these forgetting mechanisms, the passive knowledge forgetting phenomenon occurs, i.e. the knowledge carried by the past samples will be gradually forgotten as the realtime learning process goes on [27]. The directional forgetting factor discussed in [19] restricts forgetting to the informationrich subspace, and suppresses the passive knowledge forgetting phenomenon to some extent. However, it is computationally complex and thus inappropriate for real-time applications. Moreover, with these forgetting factor RLS (FFRLS) methods, the learning performance of the RBFNN will still be affected by the spatial distribution of the samples. It is still possible for RBFNN to allocate its approximation capabilities too much to the regions where the samples are relatively dense, leading to poor learning of the knowledge carried by the samples within the regions where the samples are sparse.
For RLS based training methods whose objective function is the sum of the squared approximation errors of each sample, the learning performance of the RBFNN is also influenced by the distribution of samples in the feature space. To obtain the function approximation which has good generalization capability over the feature space, the samples learned by the RBFNN should be uniformly distributed over the feature space while the number of the samples should be large enough [28]. For offline learning scenarios, we can preprocess the training set such that the samples are uniformly distributed. However, for real-time applications where the samples are mostly collected in chronological order with the same sampling period, it is difficult to distribute the samples uniformly. In this case, the RBFNN cannot allocate its approximation capabilities uniformly to the feature space such that the approximation accuracy over different regions of the feature space is different.
From the aforementioned analyisis, the existing real-time training methods for RBFNNs based on SGD or RLS have the following problems to be solved: i) there is no method to suppress passive knowledge forgetting and data saturation at the same time with sufficiently little computational complexity, ii) with these methods, the RBFNN cannot allocate its approximation capabilities uniformly to the feature space, leading to the large difference in the approximation accuracy over different regions of the feature space. To address these problems, this paper proposes an RLS based real-time training method with a mechanism of memory, which is named as the selective memory recursive least squares (SMRLS) method. The design of SMRLS based training method of the RBFNN includes the following three steps: 1) The feature space of the RBFNN is normalized and then evenly discretized into a finite number of partitions, which is named as the finite-state implementation of the feature space. 2) In a real-time learning process, the memory update functions are designed according to a priori knowledge to synthesize samples within the same partition into one sample.
3) The synthesized objective function is obtained by summing up the squared approximation errors of the synthesized samples to replace the original objective function in RLS method. Similar to the ordinary RLS method, a recursive algorithm is proposed to calculate the optimal weights of the RBFNN minimizing the synthesized objective function in the real-time learning process. As an example, specific memory update functions are proposed based on two classical assumptions respectively: i) the latest samples are more reliable than the past ones [27], ii) the measurement is accurate while the function to be approxiamted is time-invariant. Based on these memory update functions, the SMRLS method is simplified such that its computational complexity is reduced to the level of the ordinary RLS method.
The discretization of the feature space and the calculation of the synthesized objective function according to memory update functions are collectively called the memory mechanism of SMRLS. Thanks to ths mechanism, SMRLS obtains the following merits over the SGD and RLS methods: 1) Through the finite-state implementation, the maximum number of the synthesized samples is set to the number of the partitions, so the data saturation phenomenon is suppressed. The knowledge of the samples within the same partition will be synthesized rather than gradually forgotten, so passive knowledge forgetting is also suppressed. This improves the learning speed and accuracy of the RBFNN. 2) Since the samples within the same partition only contribute one synthesized sample to the objective function, the RBFNN allocates its approximation capabilities uniformly to the feature space, i.e. the approximation effects over different regions of the feature space are similar. Therefore, the generalization capability of the learned knowledge is improved. The rest of this paper is organized as follows: Section II introduces preliminaries and formulates the problem. In Section III, the passive knowledge forgetting phenomenon of SGD and FFRLS methods is analyzed. Section IV provides the design process of SMRLS method in detail and analyzes its merits. Section V provides the simulation results to demonstrate the effectiveness of SMRLS. Finally, Section VI concludes the main results.

II. PROBLEM FORMULATION AND PRELIMINARIES
In Section II-A, an open-loop real-time function approximation problem is introduced. Section II-B, II-C, II-D provide a brief introduction to the RBFNN based function approximation, SGD method and RLS method, respectively.

A. Real-Time Function Approximation: An Open-Loop Case
Consider an open-loop real-time approximation task where the output of the approximator will not influence the state orbit of the system. The function to be approximated is described as follows: defined on the compact set Ω x , y ∈ R is the measurement of f (x), and v ∈ R is the measurement error including the measurement noise. The state orbit starting from the initial state x(0) is denoted by µ(x(0)) during the learning task.

Assumption 1.
To facilitate understanding of the theoretical analysis in subsequent sections, it is assumed that the measurement is accurate such that v = 0.
Assume that the output of the approximator isf (x, θ) with the parameter θ, and the following objective function is considered: Let θ * denote the optimal parameter minimizing J(θ), and we have: Let i = 1, 2, . . . , m denote the serial number of the sampling time, m denote the total number of samples, and the real-time approximation task is to constantly improve the approximation of θ * as samples (x(i), y(i)) keep being collected. The following sum of squared error (SSE) objective function is adopted to substitute J(θ) to obtain an estimation of θ * . As analyzed in [28], if the samples (x(i), y(i)) are uniformly distributed over the compact set Ω x under Assumption 1, we obatin: lim which indicates the importance of the uniform distribution of samples for the approximation. Therefore, θ * can be estimated using: if there are sufficient samples uniformly distributed over Ω x . The compact set Ω x consisting of all possible values of x is named as the feature (input) space of the approximator. In practice, the scope of Ω x is designed to be the scope over which the approximationf (x, θ) is expected to be accurate. According to the aforementioned analysis, to obtain an accurate estimation of θ * , the samples (x(i), y(i)) are required to fill the feature space Ω x , which is named as a sufficient distribution. The sufficient distribution can be satisfied using the following methods: i) for active approximation scenarios where the state orbit µ(x(0)) can be designed, µ(x(0)) is designed properly such that the feature space is filled with samples; ii) for passive approximation scenarios, although µ(x(0)) cannot be designed, we can still collect sufficient samples by long-term sampling.
Even if the sufficient distribution of samples is satisfied, for real-time approximation tasks, the approximation method (6) is still difficult to realize because: i) the direct optimization of J f (θ) is computationally complex; ii) the uniform distribution of the samples is hard to guarantee. While there is still no effective method to address ii), some real-time approximation methods have been proposed to address i) including the commonly used SGD and RLS methods [15].
Remark 1. The real-time approximation problem discussed in this section is an open-loop case where the output of the approximator will not influence the state orbit µ(x(0)). In this case, the only issue to be considered for the design of the approximation method is how to accurately approximate the function. However, when the function approximator and the original system constitute a closed-loop system, the stability of the system should also be considered. Although the closed-loop function approximation methods are more complex, they are still based on the open-loop methods. Therefore, the real-time function approximation methods discussed in this paper can also be applied to closed-loop systems with some modifications used to guarantee the stability.

B. RBFNN based Function Approximation
The RBFNN is a class of single hidden layer feedforward neural networks whose activation function is the radial basis function. The output of the RBFNN can be formulated as: where χ ∈ R l is the input vector, f N N (χ) ∈ R is the corresponding output, W = [w 1 , w 2 , . . . , w N ] T ∈ R N is the weight vector, N is the number of neurons, and Φ(χ) = [φ 1 (χ), φ 2 (χ), . . . , φ N (χ)] T ∈ R N is the regressor vector composed of radial basis functions φ i (χ), i = 1, 2, . . . , N . The radial basis function employed in this paper is the Gaussian function as follows: where c i ∈ R l is the center and σ i ∈ R is the receptive field width of the radial basis function φ i (χ).
Consider the function approximation problem in Section II-A. RBFNNs have the ability to approximate any continuous function f (x) : R n → R over the compact set Ω x ⊆ R n to arbitrary accuracy on the premise of sufficient neurons and suitably selected centers and receptive field widths [1]. Then f (x) can be approximated as: where W * ∈ R N is the optimal weight vector of the approximation, and (x) is the bounded approximation error satisfying | (x)| < * , ∀x ∈ Ω x . In this paper, the optimal value of W * ∈ R N is defined as follows: where the compact set Ω x is known as the feature space (input space) of the RBFNN. According to (4), (6), if there are sufficient samples uniformly distributed over Ω x , W * can be estimated by minimizing J f (W ) as follows: Assumption 2. It is assumed that, the RBFNN has a lattice distribution of the neuron centers, which covers the (normalized) feature space of the RBFNN.
The persistent excitation (PE) condition is an important concept in the function approximation problems based on the linearly parameterized approximators including the RBFNN and is defined as follows. [29]): The piecewise-continuous and uniformly-bounded regressor vector Φ(x) : R n → R N , is said to satisfy the PE condition, if there exist positive constants δ, α 1 , α 2 such that

Definition 1. (PE Condition
holds for ∀t 0 0, where I ∈ R N ×N is an identity matrix and α 1 , α 2 are known collectively as the PE levels.
In this paper, the amount of knowledge learned by the RBFNN is measured through the memory of samples, which is defined as follows.
where m > 0 is a small approximation error upper bound.

C. Stochastic Gradient Descent based Training for RBFNNs
To approximate W * in real-time, the following squared error objective function is considered: where k is the current sampling time. The gradient of J GD (W, k) corresponding to W is formulated as: Then W is updated in the direction of −∇J GD (W, k) as follows: where η ∈ R + is the designed learning rate.
The following theorem indicates the rationality of the weight update law (16).
Theorem 1. Assume that the sampling period is small enough. Consider the weight update law (16), if the regressor vector Φ(x) satisfies the PE condition and the learning rate η is designed properly, W will exponentially converge to a small neighborhood around W * .
Proof. Substitute (9), (15) into (16) and we obtain: whereW = W * − W is the estimation error of the optimal weight vector W * . Since the sampling period is small enough, (17) can be reformulated into the continuous form as follows: where η 0 ∈ R N ×N is the designed learning rate which is proportional to η. SinceW = W * − W , we have: Since Φ(x) is bounded and (x) is very small, system (19) can be regarded as a nominal system with a small perturbation. The nominal part of (19) is described as: According to Theorem 2.5.1 in [15], with the PE condition of Φ(x), (20) is exponentially stable. To analyze the convergence of system (19), Lemma 9.2 in [30] is considered. This nonvanishing perturbation lemma indicates that with properly designed η 0 , the small perturbation term η 0 Φ(x) (x) will not lead to large steady-state deviations ofW from zero. Therefore, under the PE condition of Φ(x), W will exponentially converge to a small neighborhood around its optimal value W * with properly designed learning rate η.

D. Recursive Least Squares based Training for RBFNNs
Different from many other neural networks, which are highly nonlinear in the parameters and require nonlinear optimization methods for learning, the RBFNN can learn with linear optimization methods such as the RLS method thanks to its linearly parameterized form.
Consider the following SSE objective function which has been mentioned in (11): Let W * LS (k) ∈ R N denote the optimal weight vector minimizing J LS (W, k), and it can be formulated as: Let R k denote the output sample vector. It is assumed that there are sufficient samples such that the matrix H T k H k ∈ R N ×N is positive definite. Therefore, the optimal weight vector W * LS (k) can be calculated as follows: According to (11), if there are sufficient samples (x(i), y(i)), i = 1, 2, . . . , m uniformly distributed over Ω x , W * can be estimated as follows: When the samples are obtained from the online sampling of the chronological data, the RLS method can be applied to estimate W * LS (k) in real-time. Let P (k) = (H T k H k ) −1 denote the covariance matrix at sampling time k, and two equivalent RLS methods are given as follows: where λ ∈ (0, 1] is the forgetting factor used to suppress the data saturation phenomenon [31]. Let P (0) = p 0 I, W (0) = W 0 denote the initial values of P (k) and W (k) respectively, where I ∈ R N ×N is an identity matrix. According to Theorem 1 in [32], the real objective function of the RLS methods (25) and (26) is formulated as: Therefore, when using (25) or (26) to update W (k), we obtain: Let λ = 1, i.e. the RLS method is applied without the forgetting factor, to make J RLS (W, k) coincide with J LS (W, k) in (21), p 0 should be chosen as large as possible.
Assumption 3. Assume that the matrix H T k H k is positive definite such that W * LS (k) exists. In this paper, the RLS method is initialized with a sufficiently large p 0 such that the weight vector W (k) based on (25) or (26) with λ = 1 satisfies: where * k > 0 is a small estimation error upper bound corresponding to the samples and RBFNN.
The following lemma shows the derivation of the ordinary RLS method which does not adopt a forgetting factor. Lemma 1. Consider the SSE objective function (21), if the following constraints: hold for ∀k ∈ N + such that H T k H k is positive definite, then W (k + 1) can be calculated using W (k) as follows: Proof. According to (23), W (k) and W (k + 1) can be formulated as: Since ∈ R k+1 , we have:

III. PASSIVE KNOWLEDGE FORGETTING AND FINITE MEMORY OF RBFNNS
The passive knowledge forgetting phenomenon of SGD and FFRLS based RBFNNs is introduced in Section III-A. By contrast, the finite memory of the ordinary RLS based RBFNN is analyzed in Section III-B.

A. Passive Knowledge Forgetting Phenomenon
The passive knowledge forgetting phenomenon can be analyzed from the perspective of objective function. For SGD training method, the objective function J GD (W, k) is shown in (14). It should be noted that the squared error term y(i) − W T Φ(x(i)) 2 only appears in J GD (W, k) at sampling time i indicating that the update of the RBFNN will not aim at reducing the error for the update at sampling time i. Therefore, even if the error y(i) − W T Φ(x(i)) 2 is reduced by the update at sampling time i, it may increase in the subsequent training process. Similarly, for FFRLS training method whose objective function J RLS (W, k) is shown in (27), the weight of the squared error term y(i) − W T Φ(x(i)) 2 , i.e. λ k−i , will gradually decrease as the real-time learning process goes on. The term "passive knowledge forgetting" is used to summarize this phenomenon in which the past samples will gradually lose the attention of the approximator as the learning process goes on. Consider the real-time approximation problem of f (x) mentioned in Section II-A. To facilitate the understanding, this passive knowledge forgetting phenomenon is presented with a special case in Fig. 1 T ∈ R 2 is the input of the RBFNN, and the state orbit of x is divided into two segments: i) an explored segment µ a and ii) a new segment µ b to be explored. Let W * a denote the optimal weight vector of the RBFNN along µ a and we have: Let Ω a denote the compact set consisting of all possible values of x along µ a . In this case, f (x) can be locally approximated as follows: It is assumed that after real-time learning along µ a using SGD or FFRLS methods, the weights of the RBFNN have converged to W * a and thus realize locally accurate approximation of f (x) along µ a . However, when the segment µ b is sampled and learned by the RBFNN, the converged weights will gradually deviate from W * a and thus lose the accurate approximation of According to Definition 2, passive knowledge forgetting will make the RBFNN lose its memory of past samples but obtain memory of latest samples. Therefore, when SGD or FFRLS methods are used to train the RBFNN, the knowledge learned by the RBFNN cannot be accumulated effectively as learning goes on, which leads to the unqualified real-time learning performance.

B. Finite Memory of the Ordinary RLS based RBFNN
For ordinary RLS method which does not adopt a forgetting factor, the squared approximation error of each sample always has the same weight in the objective function. Therefore, in real-time learning, the RBFNN will evenly allocate its approximation capabilities to each sample and tries to keep memory of as many samples as possible. The following lemma shows this memory characteristic of the ordinary RLS based RBFNN.  (25) or (26) with λ = 1 is used to train the RBFNN to approximate f (x) under Assumption 3. If the RBFNN has reasonable hyperparameter settings, then it can keep good memory of the past samples as described in Definition 2. Specifically, assume that the serial number of current sampling time is k, and then there exists a small constant * m > 0 such that the following inequation: holds for ∀i k m, where m ∈ N + is a finite integer.
Proof. Under Assumption 1, the approximation error of W * LS defined in (22) can be formulated as follows: where W * is the optimal weight vector defined in (10). Substitute (9) into (39) and we obtain: where * > 0 is the approximation error upper bound that can be arbitrarily reduced by adjusting the hypermeters of the RBFNN. Consider the difference between ordinary RLS and least squares (LS) caused by the initial value p 0 of the covariance matrix P (k) mentioned in Assumption 3, and the approximation error of the RLS based RBFNN satisfies: (41) According to (29), the following inequation is considered: where N is the number of the neurons, * k > 0 is the small estimation error upper bound caused by the initial value p 0 of P (k).
Substitute (40), (42) into (41) and we obtain: Therefore, for ∀i k m, the approximation error of the past sample (x(i), y(i)) satisfies: (44) Since * k , * can be arbitrarily reduced by designing p 0 and the hypermeters of the RBFNN while k m is a finite value, the term kN 2 * 2 k + k * 2 in (44) can be arbitrarily reduced. Let * k = kN 2 * 2 k + k * 2 and this ends the proof.
Lemma 2 indicates that the ordinary RLS based training method makes the RBFNN keep memory of every past sample and thus overcomes the passive knowledge forgetting. However, the ordinary RLS still has the following problems to be addressed, which limits the real-time learning performance of the RBFNN: 1) The approximation error upper bound in (44) is proportional to the number of samples. In practical applications, the memory of the RBFNN is finite, so when the number of the samples increases continuously, the approximation capabilities of the RBFNN will gradually become insufficient to express the excessive knowledge. Therefore, when the number of the samples is large enough, the RBFNN will almost no longer learn new samples to retain the memory of the learned samples, which is named as the data saturation phenomenon.
2) As analyzed in Section II-A, to guarantee the generalization capability of the learned knowledge, the samples included in the objective function should be uniformly distributed over the feature space, i.e. the RBFNN should allocate its finite memory uniformly to the feature space. However, with ordinary RLS, the distribution of samples is determined by the state orbit µ(x(0)) of the input vector x, resulting in the different approximation accuracy over different regions of the feature space.

IV. SELECTIVE MEMORY RECURSIVE LEAST SQUARES
To overcome the aforementioned passive knowledge forgetting and data saturation while allocating the finite memory of the RBFNN uniformly to the feature space, a real-time training method named selective memory recursive least squares (SMRLS) is proposed in this section.

A. Finite-state Implementation of the Feature Space
Although the ordinary RLS method essentially minimizes the objective function (27) with λ = 1 using a recursive method, Lemma 1 indicates that the derivation of RLS actually considers optimization of the SSE objective function (21). Therefore, to obtain a training method uniformly allocating the finite memory of the RBFNN to its feature space, the following SSE objective function is considered: where k is the serial number of the current sampling time, M (k) is the number of samples at sampling time k, and the sampled input x(j), j = 1, 2, . . . , M (k) is uniformly distributed over the feature space Ω x for ∀k ∈ N + . The most prominent characteristic of J LS (W, M (k)) is that its samples are always uniformly distributed over the feature space during the real-time learning process.
To obtain an objective function in the form of (45), the following steps, which are collectively named as the finitestate implementation of the feature space, are adopted: 1) The feature space Ω x of the RBFNN is normalized into a certain scope Ωx. For example, Fig. 2 shows a twodimensional case of the normalization into the scope of 1]. In this case, the original input of the ∈ Ω x is normalized as follows: T ∈ Ωx is the normalized input. The reason for the normalization is to simplify the neuron center distribution of the RBFNN such that f (x) along different state orbits can be approximated using the same neuron distribution satisfying Assumption 2 as follows: whereW * is the optimal weight vector defined as: 2) The normalized feature space Ωx is evenly divided into N P partitions which are mutually disjoint as shown in Fig. 2. Let Ω j x , j = 1, 2, . . . N P denote the N P partitions, and the samples within each partition Ω j x are synthesized into one sample (γ j (k), ϕ j (k)) in real-time. Assume that there are M P (k) partitions which have been sampled before sampling time k and the number of partitions N P is large enough, if the synthesized samples satisfy γ j (k) ∈ Ω j x , j = 1, 2, . . . , M P (k), the uniform distribution of the synthesized samples is guaranteed. To facilitate the analysis, in a certain real-time learning process, the serial number j of a partition is assumed to be the order in which this partition is first sampled. For example, the partition in which the first sample (x(1), y(1)) is collected is denoted by Ω 1 x . With the finite-state implementation, the number of the synthesized samples at sampling time k is equal to the number of the sampled partitions M P (k). Therefore, at sampling time k + 1, the number of the synthesized samples M P (k + 1) satisfies: where Ω j x ⊂ Ωx is the union set of the partitions which have been sampled at least once before sampling time k. Therefore, the SSE objective function of the synthesized samples can be formulated as follows: which is named as the synthesized objective function and coincides with the form of (45).
Since (γ j (k), ϕ j (k)) represents the synthesis of the samples within the partition Ω j x at sampling time k, its value will be updated only when the samples within Ω j x increase. Therefore, a partition judgment function is defined as follows: The partition judgment function G(j,x(k)) is used to determine in which partition the current sample (x(k), y(k)) locates and thus to determine which synthesized sample should be updated. Based on the partition judgment function, the update laws of γ j (k) and ϕ j (k) are proposed as follows: γ j (k + 1) = γ j (k) + G(j,x(k + 1))S γ (γ j (k),x(k + 1)) ϕ j (k + 1) = ϕ j (k) + G(j,x(k + 1))S ϕ (ϕ j (k), y(k + 1)) (52) where S γ (·) and S ϕ (·) are memory update functions representing the update mechanism of the synthesized samples. Specially, if the new knowledge (x(k+1), y(k+1)) is believed to be more reliable than the old knowledge (γ j (k), ϕ j (k)), the memory update functions S γ (·), S ϕ (·) can be designed as: Remark 2. It should be noted that the implementation of the memory update law (52) requires the initialization of the synthesized sample (γ j (k), ϕ j (k)). In this paper, the initial value of (γ j (k), ϕ j (k)) is set to the first sample within the partition Ω j x during the online learning process. Let k j k denote the sampling time when partition Ω j x is first sampled, and the initialization of (γ j (k), ϕ j (k)) is described as follows: where j = 1, 2, . . . , M P (k) represent the serial numbers of the partitions which have been sampled at least once before sampling time k.
Through the normalization and partitioning operations, the finite memory of the RBFNN is uniformly allocated to the normalized input space Ωx and the maximum quantity of the synthesized samples is limited to the total number of the partitions, i.e. N P . With the finite-state implementation, the aforementioned uniform sample distribution and data saturation problems of the ordinary RLS method are addressed.

B. Recursive Optimization for Synthesized Objective Function
To calculate the weight vector W (k) minimizing J S (W, M P (k)) in real-time with little computational complexity, a method named selective memory recursive least squares (SMRLS) is proposed. The derivation of SMRLS is shown in detail as follows.
Assumption 5. It is assumed that the sample (x(k + 1), y(k + 1)) is collected from the partition Ω ā x at sampling time k + 1. Consider the objective function (50) and the memory update functions (52). Different from the SSE objective function of the ordinary LS method in (21) whose samples continuously increase at each sampling time, the objective function (50) of SMRLS is about the synthesized samples (γ j (k), ϕ j (k)), j = 1, 2, . . . , M P (k) whose number does not necessarily increase as the learning process goes on. At the sampling time k + 1, the behaviors of J S (W, M P (k)) can be classified into two categories: addition or replacement. Under Assumption 5, these two behaviors can be explained as follows: 1) The addition behavior means that if the sample (x(k + 1), y(k + 1)) is the first sample within the partition Ω ā x , a new synthesized sample (γ a (k + 1), ϕ a (k + 1)) is generated and a new term ϕ a (k + 1) − W T Φ(γ a (k + 1)) 2 will be added to the synthesized objective function (50). In this case, a = M P (k + 1) = M P (k) + 1.
2) The replacement behavior means that if the sample (x(k + 1), y(k + 1)) is not the first sample within the partition Ω ā x , the synthesized sample (γ a (k), ϕ a (k)) will be updated to (γ a (k + 1), ϕ a (k + 1)) according to (52) and the squared error term ϕ a (k) − W T Φ(γ a (k)) 2 in (50) will also be replaced by a new term ϕ a (k + 1) − W T Φ(γ a (k + 1)) 2 . In this case, a M P (k + 1) = M P (k). To distinguish these two behaviors, a function of the latest sampled input vectorx(k + 1) is defined as follows: where Ω M P (k) is the union set of the partitions which have been sampled at least once before sampling time k.
To obtain the update law of SMRLS, the update from sampling time k to k + 1 is analyzed. According to (50), the SMRLS based weight vector W (k) and W (k+1) are assumed to satisfy the following constraints: Based on the aforementioned addition and replacement behaviors of J S (W, M P (k)), the update from W (k) to W (k + 1) can be classified into two cases as follows: 1) Ifx(k + 1) / ∈ Ω M P (k) , which means that the sample (x(k + 1), y(k + 1)) is the first sample within the partition Ω ā x , the aforementioned addition behavior occurs. Therefore, we have a = M P (k + 1) = M P (k) + 1, and the synthesized objective function J S (W, M P (k + 1)) can be described as: Consider (56), (57), and the constraints of the update can be summarized as: which coincides with the constraints of the ordinary RLS as analyzed in Lemma 1. Therefore, ifx(k + 1) / ∈ Ω M P (k) , the update from W (k) to W (k + 1) can adopt the ordinary RLS based update law (31), (32) and be formulated as follows: S (k + 1) = P −1 S (k) + Φ(γa(k + 1))Φ T (γa(k + 1)) Φ(γa(k + 1)) = [φ1(γa(k + 1)), . . . , φN (γa(k + 1))] T (59) where P S (k) is the corresponding covariance matrix. According to the analysis in Section II-D, the initial value of P S (k) is set to p 0 I, where I ∈ R N ×N is an identity matrix and p 0 > 0 is as large as possible.
According to Remark 2, the value of (γ a (k + 1), ϕ a (k + 1)) is set to (x(k + 1), y(k + 1)), and thus (59) becomes: (60) 2) Ifx(k + 1) ∈ Ω M P (k) , which means that (x(k + 1), y(k + 1)) is not the first sample within the partition Ω ā x , the aforementioned replacement behavior occurs. Therefore, we have a M P (k + 1) = M P (k), and the synthesized objective function J S (W, M P (k + 1)) can be described as: To obtain the recursive update law of W (k), a virtual synthesized objective function, which involves the synthesized samples of the partitions apart from Ω ā x , is proposed as follows: Since J V S (k + 1) = J V S (k), substitute (62) into (61), and we have: Similar to J S (W, M P (k + 1)), the synthesized objective function J S (W, M P (k)) can be described as follows: denote the virtual weight vector minimizing the objective function J V S (k) as follows: Consider (56), (63), (64), (65), and we have the following update constraints: J S (W, M P (k + 1)) = J V S (k + 1) + ∆J S (k + 1).

Remark 3.
It should be noted that the implementation of (70) requires storage of the past synthesized samples of each partition. At sampling time k + 1, each partition Ω j x , j = 1, 2, . . . , M P (k + 1) should hold the memory of (γ j (k), ϕ j (k)) to realize the recursive calculation. Assume that the number of partitions N P is large enough such that Φ(γ a (k + 1)) − Φ(γ a (k)) is very small and Φ(γ a (k)) can be approximated by Φ(γ a (k + 1)). Substitute Φ(γ a (k)) = Φ(γ a (k + 1)) into (70), and we obtain: (72) which only requires the storage of ϕ j (k). To reduce the computational complexity, (72) can be used to replace (70) in the recursive optimization process. The rationality of this simplification is based on the premise that the number of partitions N P is large enough, so when N P is small, this simplification will reduce the learning accuracy of SMRLS.

Remark 4.
Since the maximum number of the samples in SMRLS is limited to N P , i.e. the number of partitions, the accuracy of SMRLS is related to the partitioning operation. When N P is too small, the accuracy of SMRLS will be affected while too many partitions will increase the computational complexity. Therefore, N P should be designed properly according to specific learning tasks and partitioning methods other than evenly partitioning deserve further exploration.

C. SMRLS Based on Common Assumptions
The implementation of SMRLS requires the specific design of the memory update functions in (52) based on a priori knowledge of the system. Different memory update functions will lead to very different characteristics of the SMRLS based methods. For real-time function approximation scenarios, there is a classical assumption that the latest samples are more reliable than the past ones, which is also the basis of the SGD and FFRLS methods [27]. Interestingly, SMRLS can also be simplified with this assumption and the derivation of the simplified method is shown in detail as follows.
Assumption 6. In the SMRLS based real-time function approximation process, the latest sample within the partition Ω j x , j = 1, 2, . . . , M P (k) can represent the unknown function more accurately compared with the past samples within the same partition.
(76) Similar to the analysis of Remark 3, when N P is large enough, SMRLS-R can also be simplified as follows: (77) which reduces the computational complexity of SMRLS-R.
In many practical systems, the unknown functions to be approximated are time-invariant, and thus the following assumption is also considered.
Assumption 7. The function to be approximated by the RBFNN is time-invariant and the measurement is accurate.
Under Assumption 7, each sample within the same partition is regarded to have the same effect to characterize the unknown function. Therefore, the first sample within the partition Ω j x , j = 1, 2, . . . , N P is taken as the memory of Ω j x and not updated in the subsequent learning process. In this case, the memory of each partition is fixed after the learning of the first sample within the partition, and thus the method is named as fixed SMRLS (SMRSL-F). Consider (71), and we only need to stop the update when ρ(x(k+1)) = 1 to obtain the description of SMRLS-F: ∈ Ω M P (k) : 2) if ρ(x(k + 1)) = 1, i.e.x(k + 1) ∈ Ω M P (k) : (78) In this case, the memory update functions in (52) are regarded as S γ (γ j (k),x(k+1)) = S ϕ (ϕ j (k), y(k+1)) = 0 after the initializatin of (γ j (k), ϕ j (k)). Compared with the ordinary RLS method, SMRLS-F obtains lower computational complexity while its accuracy is influenced by the number of partitions according to Remark 4.

D. Merits of Selective Memory Recursive Least Squares
Compared with the classical RLS or SGD based training methods, SMRLS has the following merits: 1) Selective Memory of the RBFNN: Since the derivation of SMRLS is based on the ordinary RLS method, the SMRLS based RBFNN also has the finite memory mentioned in Lemma 2 such that all of the knowledge (samples) will be synthesized and learned by the RBFNN during the real-time learning process. In this way, with properly designed memory update functions, the passive knowledge forgetting phenomenon is effectively suppressed.
In addition, the finite-state implementation of the feature space also suppresses the data saturation problem of the ordinary RLS method. According to (71), if ρ(x(k + 1)) = 0, the eigenvalues of P −1 S (k + 1) satisfies: where λ i (P −1 S (k)), i = 1, 2, . . . , N represent the eigenvalues of P −1 S (k) and tr(·) is the trace operator. Since P −1 S (k) is positive definite, i.e. λ i (P −1 F (k)) > 0, the eigenvalues of P −1 S (k) will gradually increase as more partitions are explored and thus the eigenvalues of P S (k) will gradually decrease. The term P S (k+1)Φ(x(k+1)) in case 1) of (71), which represents the learning rate, satisfies the following constraint: PS(k + 1)Φ(x(k + 1)) 2 λ 2 max (PS(k + 1)) Φ(x(k + 1)) 2 N λ 2 max (PS(k + 1)) = N/λ 2 min (P −1 S (k + 1)) (80) where N is the number of the neurons, λ max (·) and λ min (·) represent the maximum and minimum of the eigenvalues, respectively. According to (71), the learning rate P S (k + 1)Φ(x(k + 1)) will gradually decrease as more partitions are explored, which means that the finite memory of the RBFNN is being consumed. With the ordinary RLS method, this memory consumption happens whenever a new sample is learned by the RBFNN and thus the weights W (k) will nearly stop updating when the finite memory is depleted, i.e. the data saturation occurs. By contrast, with the SMRLS method, the maximum number of the synthesized samples is set to N P , so the finite memory of the RBFNN will not be depleted during the real-time learning process. Therefore, SMRLS also effectively suppresses the data saturation phenomenon.
Since both the passive knowledge forgetting and the data saturation are suppressed through the selective memory, the real-time learning speed and accuracy of the RBFNN are improved.

2) Uniformly Allocated Approximation Capabilities:
As analyzed in Section II-A, the uniform distribution of the samples over the feature space is essential for the generalization capability of the learned knowledge. With SMRLS, if the number of partitions N P is large enough, the synthesized samples (γ j (k), ϕ j (k)), j = 1, 2, j = 1, 2, . . . , M P (k) can be regarded to be uniformly distributed over the feature space. Therefore, the SMRLS based RBFNN allocates its finite approximation capabilities uniformly to the feature space and thus obtains good generalization capability.

V. SIMULATION STUDIES
In this section, the RBFNN is employed to approximate an unknown function of a single inverted pendulum-cart system, whose state space description is formulated as follows: .
The state variables x 1 and x 2 represent the angular position and velocity of the pendulum, respectively, g = 9.8m/s 2 is the gravitational acceleration, m c = 0.1kg is the weight of the cart, m = 0.02kg is the weight of the pendulum, and l = 0.2m is half the length of the pendulum. The RBFNN is used to approximate f I (x) in real-time.
In this section, the following real-time learning methods are compared: SMRLS-R (76), ordinary RLS (25)  It should be noted that this paper focuses on the realtime learning problem of the RBFNN rather than the control. Therefore, the RBFNN is not involved in the closed-loop system and it is assumed that the output of the unknown function f I (x) to be approximated can be measured directly as follows: The objective of the real-time learning task is to obtain an accurate approximation of f I (x) over the normalized feature space Ωx using the increasing samples (x(k), y(k)).

A. Simulation Studies for Data Saturation
To test the suppression effect of different learning methods to data saturation, the following two real-time learning scenarios are considered.

(A) Abrupt parameter perturbation
Consider the following state orbit of x 1 (t) and corresponding measurement: where the value of l in (82) is changed abruptly as below: The training of RBFNNs is carried out over [0, 100s]. Fig. 3 and Fig. 4 demonstrate the convergence process of the approximation errors and weights, respectively. In the first 50s of the training, the RBFNNs have approximated the unknown function f I (x) with the parameter l = 0.2, and thus when the parameter abruptly changes, the data saturation phenomenon greatly reduces the convergence speed of the ordinary RLS method. Compared with SGD and FFRLS methods which can suppress data saturation, SMRLS-R achieves faster learning speed to the changed function. According to the derivation of SMRLS in Section IV, the SMRLS-R based weight convergence can be completed within only one period of the state orbit (2πs for case (A)) if the sampling period is small enough such that all of the partitions along the state orbit can be sampled at least once in a single period.
To test the accuracy of the approximation, the weights at 100s are recorded as the learned knowledge to approximate f I (x) with the latest parameter l = 0. 3. Let e f = f I (x) − W T Φ(x) denote the approximation error. As shown in Fig. 5, while all the four methods achieve the basic approximation of f I (x), the accuracy of FFRLS and SMRLS-F methods is obviously better than the other two methods. The reason why SMRLS-R is less accurate than FFRLS lies in the limitation of the maximum number of samples, i.e. N P , mentioned in Remark 4. The increase of N P will improve the approximation accuracy but also increases the computational complexity and the requirement for sampling period.

(B) Decreasing measurement error
Consider the following state orbit of x 1 (t) over [0, 100s] and corresponding measurement: where v(k) = e −0.2k sin k is the decreasing measurement error. Case (B) represents the scenarios where the output of the unknown function cannot be measured directly and requires extra estimation with the converging measurement error. For example, in a class of control systems, the output of the unknown function to be approximated can only be obtained using an observer and the objective of the RBFNN is to approximate the function with the measurements. Fig. 6 shows the convergence of the approximation error in the real-time training process, indicating that SMRLS-R still has the fastest convergence speed. The approximation effect using the learned knowledge is shown in Fig. 7, where FFRLS has the highest learning accuracy and the accuracy of SMRLS-R is higher than the ordinary RLS method.
According to the results in case (A) and (B), SMRLS-R suppresses the data saturation phenomenon effectively, and thus has faster convergence speed for time-varying learning scenarios compared with the other classical real-time training methods. However, the accuracy of its approximation is determined by the partitioning operation.

(C) Non-repetitive state orbit
To test the suppression effect of SMRLS to passive knowledge forgetting, the following non-repetitive state orbit is considered: where the duration of the learning task is set to 100s.
Since case (C) coincides with Assumption 7, SMRLS-F is also compared with the other methods in this case. Fig. 8 shows that FFRLS has the highest real-time tracking accuracy of f I (x) along the non-repetitive state orbit while SMRLS-R, ordinary RLS and SMRLS-F have similar tracking accuracy. However, high tracking accuracy does not mean high accuracy of the learned knowledge. As shown in Fig. 9, the weights of FFRLS maintain a relatively large updating speed over the training process indicating that its high-accuracy tracking of f I (x) may be at the cost of forgetting the past samples. The weight vectors at 100s are recorded as the learned knowledge and Fig. 10 demonstrates the learning accuracy of the methods. Consistent with the analysis in Section III, the passive knowledge forgetting of FFRLS and SGD improves their real-time tracking accuracy of f I (x) while reducing the accuracy of the learned knowledge because the past samples are gradually forgotten. By contrast, the SMRLS and ordinary RLS methods sufficiently utilize both the past and the latest samples such that the accuracy of the learned knowledge is guaranteed. Moreover, the results also indicate that SMRLS-F has very similar learning performance to SMRLS-R in case (C) and thus becomes a computationally efficient alternative of SMRLS-R under Assumption 7.
C. Simulation Studies for Approximation Capability Allocation (D) Non-uniformly distributed state orbit One of the most important merits of SMRLS is its uniform allocation of the approximation capabilities such that the RBFNN will pay average attention to its approximation effect in different regions of the feature space and thus obtains good generalization capability. To test this merit, the random Non-Uniform Rational B-Splines (NURBS) state orbit µ(x(0)) in Fig. 11(a) is adopted and the value of f I (x) along this orbit is shown in Fig. 12.
Since the learning of SGD is less effective than the RLS based methods and SMRLS-F has very similar learning effect to SMRLS-R, only SMRLS-R, ordinary RLS and FFRLS methods are compared in case (D). Fig. 11(b)-(d) and Fig. 13 demonstrate the update of the weight vectors and approximation errors, respectively. In the learning process, FFRLS has higher real-time tracking accuracy to f I (x) than SMRLS-R and ordinary RLS methods. In real-time approximation tasks including case (D), the stability of the weights reflects the generalization capability of the learned knowledge. Specifically, if the weights hardly update over a time interval, it means that this set of weights can approximate f I (x) accurately during this interval. Therefore, the dramatical update of the weights shown in Fig. 11(d) indicates the weak generalization capability of FFRLS. Fig. 13 demonstrates the accuracy of the learned knowledge along the same random NURBS state orbit. The four subplots (i) − (iv) of Fig. 13(b) are corresponding to the four segments in red of Fig. 11(a). While ordinary RLS mostly allocates the approximation capabilities of RBFNN to the regions where the samples are dense (i.e. the segments in blue of Fig. 11(a)) and FFRLS mostly allocates the approximation capabilities to the latest samples, SMRLS-R allocates the approximation capabilities uniformly to the feature space. Therefore, the FFRLS based approximation is more accurate over the regions filled with the latest samples than the other regions, the ordinary RLS approximates f I (x) more accurately over the regions where the samples are dense, and SMRLS-R has relatively uniform approximation accuracy along the state orbit. Particularly, SMRLS-R effectively improves the approximation accuracy over the regions (i) − (iv) in Fig.  11(a) where the samples are sparse compared with FFRLS and ordinary RLS methods.
The aforementioned simulation results indicate that the SMRLS method can suppress data saturation and passive knowledge forgetting at the same time and thus improves the real-time learning performance of the RBFNN. Moreover, through SMRLS, the finite approximation capabilities of the RBFNN are uniformly allocated to the feature sapce such that the learned knowledge has good generalization capability.

VI. CONCLUSION
In this paper, the real-time learning problem of RBFNNs is studied and a phenomenon named passive knowledge forgetting is analyzed to explain the limited learning performance of RBFNNs caused by the forgetting mechanism in SGD and RLS based training methods. To suppress both passive knowledge forgetting and data saturation, the SMRLS method is proposed by discretizing the feature space and optimizing a synthesized objective function with the recursive method. Through SMRLS, the RBFNN can allocate its finite approximation capabilities uniformly to its feature space and thus obtains improved real-time learning performance.
In addition to the RBFNN, SMRLS is also applicable for real-time training of other single hidden layer feedforward neural networks. Moreover, as a result of its good learning performance and acceptable computational complexity, SMRLS has considerable potentials in the field of adaptive control, system identification and machine learning.