Radial Basis Function

A classical RBF network is almost always described as having:

• Input Layer: Receives the input vector $x\in {\mathbb{R}}^n$.
• Feature Layer: A non-linear feature map that calculates the activation based on the distance from the input to defined “center” points.
• Output Layer: Computes a simple, weighted linear sum of the feature layer’s outputs.

To compute the function of an RBF network, we define $m$ “center points” in the input space: $c_i\in {\mathbb{R}}^n$ (where $1\le i\le m$).

The Feature Map (Hidden Layer):

For each center $c_i$, the neuron computes a component ${\phi }_i:{\mathbb{R}}^n\to \mathbb{R}$ using a Gaussian function:

${\phi }_i(x)=\exp \left(-\frac{\Vert x-c_i{\Vert }_2^2}{{\sigma }_i^2}\right)$

• $c_i$: The center of the receptive field.

• $\Vert x-c_i{\Vert }_2^2$: The squared Euclidean distance between the input data and the center.

• ${\sigma }_i^2$: Controls the width (variance) of the receptive field.

• If the input $x$ is exactly equal to the center $c_i$, the distance is $0$. The math becomes $\exp (0)=1$. The neuron outputs its maximum value.

• As $x$ moves further away from $c_i$, the negative exponent grows larger. $\exp (-\text{large number})$ approaches $0$.

The Final Output: The network produces its final prediction by calculating the weighted linear combination of all $m$ feature maps:

$\text{Output}={\sum }_{i=1}^m{w}_i\cdot {\phi }_i(x)$

---

While the Gaussian $\exp$ function is the standard, other mathematical distance functions can be used for the feature mapping, such as:

• Multiquadratics: $f(t)=\sqrt{t^2+c^2}$, where $c>0$ is constant.
• Inverse Multiquadratics: $f(t)=\frac{1}{\sqrt{t^2+c^2}}$ (Where $t$ represents the distance $\Vert x-c\Vert$)

---

When the number of RBF centers ($m$) is kept strictly smaller than the number of data points ($T$) ( $m\ll T$ ) exact interpolation is impossible. The network must instead find the optimal weights $w$ that minimize the overall prediction error.

• The Objective Function $E(w)$: We define the standard Sum of Squared Errors objective function $E(w)=\frac{1}{2}{\sum }_{t=1}^T{\left[\underset{\text{RBF Network Output}}{\underbrace{\sum _{i=1}^m{w}_i\exp \left(-\frac{\Vert x(t)-c_i{\Vert }_2^2}{{\sigma }^2}\right)}}-\underset{\text{Target}}{\underbrace{y(t)}}\right]}^2$

• The hidden layer of an RBF network is non-linear, but the output layer is purely a linear sum, finding these weights is actually much, much easier than Backpropagation

---

Before we even touch the weights or the objective function, we just look at the raw input data. We use a grouping algorithm (like K-Means Clustering) to find natural “clusters” or hotspots in the data, and we drop a center ($c$) in the middle of each cluster.

Assuming we have already chosen our $m$ center points ($c_1\dots c_m$), the task of finding the perfect weights for interpolation can be written as a single matrix multiplication problem.

Let:

• $Y$: A column vector $(T\times 1)$ containing all the true target measurements $y(1)\dots y(T)$.
• $w$: A column vector $(m\times 1)$ containing the unknown weights $w_1\dots w_m$ that we need to solve for.
• $M$: The Design Matrix (or Activation Matrix) of size ${\mathbb{R}}^{T\times m}$.

The interpolation problem requires us to find a weight vector $w$ such that:

$\left(\begin{array}{c}{y}_{(1)}\\ y_{(2)}\\ ⋮\\ y_{(T)}\end{array}\right)=M\cdot \left(\begin{array}{c}{w}_1\\ ⋮\\ w_m\end{array}\right)$

The matrix $M$ stores the activation value of every hidden neuron for every data point in the training set.

$M=\left(\begin{array}{cccc}{\phi }_1(x(1))& {\phi }_2(x(1))& \dots & {\phi }_m(x(1))\\ {\phi }_1(x(2))& {\phi }_2(x(2))& \dots & {\phi }_m(x(2))\\ ⋮& ⋮& \ddots & ⋮\\ {\phi }_1(x(T))& {\phi }_2(x(T))& \dots & {\phi }_m(x(T))\end{array}\right)$

• Rows (Data Points): Each row represents a single test point $x(t)$ from the dataset.
• Columns (Hidden Neurons): Each column represents one specific hidden RBF neuron (defined by its center $c_i$).

Calculating the optimal weight vector $w$ depends entirely on the ratio of training data points ($T$) to chosen RBF centers ($m$).

Case 1: T<m (Under-determined)

• The system has fewer equations than unknowns.
• There are infinitely many exact solutions.

Case 2: $T=m$ (Square System / Strict Interpolation)

• The number of data points exactly matches the number of centers. As long as all chosen centers $c_i$ are distinct, matrix $M$ is guaranteed to be non-singular (invertible).
• $w=M^{-1}\cdot Y$ Case 3: $T>m$ (Over-determined / Approximation)
• The system has more data points than centers.
• Matrix $M$ is rectangular ($T\times m$) and therefore has no standard inverse. Exact interpolation is impossible.
• We must find the weights that minimize the Sum of Squared Errors. We solve this using the Moore-Penrose Pseudo-Inverse ($M^{+}$): $w=(M^{T}M{)}^{-1}{M}^T\cdot Y$

The scalar summation form of the error function is:

$E(w)=\frac{1}{2}{\sum }_{t=1}^T((Mw{)}_t-y(t){)}^2$

Using the dot product property of vectors ($v^{T}v=\sum v_i^2$) with error vector be $(y-Mw)$:

$E(w)=\frac{1}{2}(y-Mw{)}^T(y-Mw)$

$E(w)=\frac{1}{2}\left(y^{T}y-2y^{T}Mw+w^T{M}^{T}Mw\right)$

• $y\in {\mathbb{R}}^{T\times 1}$: Column vector of true target values.
• $M\in {\mathbb{R}}^{T\times m}$: The RBF activation matrix.
• $w\in {\mathbb{R}}^{m\times 1}$: Column vector of the unknown weights.

The Necessary Optimality Condition:

$\frac{\partial E}{\partial w}(w^{\ast })=-M^{T}y+M^{T}M{w}^{\ast }=0$

$M^{T}M{w}^{\ast }=M^{T}y$

Solving for $w^{\ast }$: If the resulting square $m\times m$ matrix $(M^{T}M)$ is invertible

$w^{\ast }=\underset{\text{Pseudo-Inverse of }M}{\underbrace{(M^{T}M{)}^{-1}{M}^T}}y$

$(M^{T}M{)}^{-1}{M}^T$, is known as the Moore-Penrose Pseudo-Inverse (often denoted as $M^{+}$). It allows us to mathematically find the “line of best fit” for a rectangular matrix.

• Singular Value Decomposition: In general, if the matrix $M^{T}M$ is singular one can follow similar approaches to define a pseudo-inverse using the Singular Value Decomposition (SVD) of the matrix.
• Regularization: Another common way to regularize the problem (fix the singularity/instability) is to add a small mathematical penalty $\lambda$ to the diagonal of the matrix before inverting it: $w^{\ast }=(M^{T}M+\lambda I_d{)}^{-1}{M}^{T}y\text{for }\lambda >0$ (Where $I_d$ is the Identity Matrix. In machine learning, this trick is commonly known as Ridge Regression!)

---

The RBF Learning Algorithm

Choose the Centers ($c_i$)

• The center points of the hidden neurons must be fixed before weight calculation begins.
• Method: This is usually done using various heuristics. A common, simple approach is to randomly select a subset of the actual test points $x(t)$ to serve as the centers. More advanced methods involve clustering algorithms like K-Means.
  • Orthogonal Least Squares (OLS) Reduction: (Chen, Billings, Cowan 1989/1991). A forward-selection algorithm that builds the RBF hidden layer iteratively. It selects data points to serve as centers one by one, choosing the point that provides the maximum reduction in the residual error at each step. Gram-Schmidt Orthogonalization: The OLS algorithm orthogonalizes the chosen feature vectors at each step. This guarantees that each newly added center captures unique variance in the data that previous centers failed to explain. Like Principal Component Analysis OLS acts as a dimensionality reduction technique, identifying the most critical “components” (centers) that represent the underlying structure of the data.

Choose the Spread/Width (${\sigma }_i$)

• Once the centers are placed, you must define the variance for each Gaussian receptive field.
• Method: This determines how strongly nearby points activate the neuron. Often, a single global $\sigma$ is chosen based on the average distance between the chosen centers.
• The $d_{max}$ Heuristic:

${\sigma }_i=\sigma =\frac{d_{max}}{\sqrt{m}}$

- $d_{max}$: The maximal Euclidean distance between any two chosen centers in the network.
- $m$: The total number of chosen centers.

Solve for the Weights ($w^{\ast }$)

• With the hidden layer locked, the network must find the optimal linear weights connecting the hidden layer to the output.
• The optimal weights can be found using the Pseudo-Inverse equation (): $w^{\ast }=(M^{T}M{)}^{-1}{M}^{T}y$.
• Computing a matrix inverse is computationally prohibitive (scaling at $O(N^3)$ complexity).
• To avoid computing the inverse, we start with a randomized initial weight vector $w^{(0)}$ and use an iterative numerical optimization algorithm (such as Steepest Descent / Least Mean Squares) to approximate $w^{\ast }$.

Orthogonal Least Squares (OLS) algorithm The goal of OLS is to compute an “energy function” (a measure of variance or information) for the network. By evaluating this function, we can determine exactly how much a specific test point contributes to reducing the overall error if it is selected to be a center.

The Pruning Strategy:

1. Start by defining the activation matrix $M$ as if every single test point is a center (yielding a $T\times T$ square matrix).
2. Evaluate the energy contribution of each column.
3. Systematically remove the columns (centers) that contribute the least energy, reducing the dimension from $T$ down to $m$.

For any given choice of weight vector $w$, the network outputs a prediction vector $\hat{y}=Mw$. The difference between the true target values ($y$) and the network’s prediction is the error vector ($e$). $y=M\cdot w+e$

• The Subspace: The columns of the design matrix $M$ span a specific subspace in ${\mathbb{R}}^T$. This subspace represents all possible outputs the current RBF network can mathematically generate.
• The Projection ($\hat{y}$): The optimal network prediction $\hat{y}$ is exactly the orthogonal projection of the true target vector $y$ onto the column space of $M$.
• The Error ($e$): Because it is an orthogonal projection, the optimal error vector $e$ is perfectly perpendicular (orthogonal) to the column space of $M$.

Because $\hat{y}$ is a projection, we can use the Pythagorean theorem (or Gram-Schmidt orthogonalization) to cleanly separate the “energy” (variance) of $\hat{y}$ into individual chunks provided by each individual column of $M$. We then simply select the centers whose columns provide the largest chunks of energy.

To analytically solve the RBF network without explicitly computing the inverse of the computationally heavy matrix $(M^{T}M)$, we use the Gram-Schmidt orthogonalization algorithm to perform a QR Decomposition.

We decompose the activation matrix $M\in {\mathbb{R}}^{T\times m}$ into two separate matrices:

$M=Q\cdot R$

• $Q\in {\mathbb{R}}^{T\times m}$: An orthogonal matrix where every column is pairwise orthogonal to the others.

• $R\in {\mathbb{R}}^{m\times m}$: An upper triangular matrix with entries of $1$ on the main diagonal.

$Q^{T}Q=\left(\begin{array}{ccc}{h}_1& 0& 0\\ 0& h_2& 0\\ 0& 0& h_m\end{array}\right)$

(Where $h_i=q_i^T{q}_i\ne 0$, representing the squared norm/energy of each orthogonal vector).

$y=(Q\cdot R)\cdot w+e$ a new intermediate weight vector $g=R\cdot w$

$y=Q\cdot g+e$ To find the optimal vector $g$ that minimizes the squared error, we apply the standard pseudo-inverse derivation

$g^{opt}=(Q^{T}Q{)}^{-1}{Q}^{T}y$

$Q^{T}Q$ is purely diagonal, its inverse $(Q^{T}Q{)}^{-1}$ is take the reciprocal of the diagonal entries $\text{diag}(h_1\dots h_m{)}^{-1}$.

$(Q^{T}Q{)}^{-1}=\left(\begin{array}{ccc}\frac{1}{h_1}& 0& 0\\ 0& \frac{1}{h_2}& 0\\ 0& 0& \frac{1}{h_m}\end{array}\right)$

Because the columns of $Q$ ($q_1\dots q_m$) are pairwise orthogonal, the optimal intermediate weight vector $g^{opt}$ can be calculated element by element.

$Q^{T}y=\left(\begin{array}{c}{q}_1^T\cdot y\\ q_2^T\cdot y\\ ⋮\\ q_m^T\cdot y\end{array}\right)$ $Q^{T}Q=\left(\begin{array}{cccc}{q}_1^T{q}_1& 0& \dots & 0\\ 0& q_2^T{q}_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & q_m^T{q}_m\end{array}\right)$

$(Q^{T}Q{)}^{-1}=\left(\begin{array}{cccc}\frac{1}{q_1^T{q}_1}& 0& \dots & 0\\ 0& \frac{1}{q_2^T{q}_2}& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & \frac{1}{q_m^T{q}_m}\end{array}\right)$

$g^{\text{opt}}=\left(\begin{array}{ccc}\frac{1}{q_1^T{q}_1}& 0& 0\\ 0& \ddots & 0\\ 0& 0& \frac{1}{q_m^T{q}_m}\end{array}\right)\cdot \left(\begin{array}{c}{q}_1^{T}y\\ ⋮\\ q_m^{T}y\end{array}\right)$

The $i$-th component is computed as: $g_i^{\text{opt}}=\left(\frac{1}{q_i^T\cdot q_i}\right)\cdot (q_i^T\cdot y)$

$g_i^{opt}=\frac{q_i^T\cdot y}{q_i^T\cdot q_i}$

For the squared error to be truly minimized, the optimal prediction vector ($Q\cdot g^{opt}$) must be perfectly orthogonal to the residual error vector ($e$).

• Mathematically: $(Q\cdot g^{opt})\perp e$

• Consequence: Their dot product is zero. $(Q\cdot g^{opt}{)}^T\cdot e=0$.

We evaluate the total variance (energy) of the target vector by taking its inner product with itself: $y^{T}y$.

$y^{T}y=(Q\cdot g^{opt}+e{)}^T(Q\cdot g^{opt}+e)$

$y^{T}y=(Q\cdot g^{opt}{)}^T(Q\cdot g^{opt})+\underset{=0}{\underbrace{(Q\cdot g^{opt}{)}^{T}e}}+\underset{=0}{\underbrace{e^T(Q\cdot g^{opt})}}+e^{T}e$

$y^{T}y=(g^{opt}{)}^T{Q}^{T}Q{g}^{opt}+e^{T}e$

Since $Q^{T}Q$ is a diagonal matrix containing elements $h_i$, the final scalar Energy Equation is:

$y^{T}y={\sum }_{i=1}^m{h}_i(g_i^{opt}{)}^2+e^{T}e$

The total explained variance of the network is the sum of the individual contributions: ${\sum }_{i=1}^m{h}_i(g_i^{opt}{)}^2$.

To find the specific variance contribution of a single candidate center $i$, we substitute the known matrix definitions into the term:

$h_i\cdot (g_i^{opt}{)}^2=(q_i^T\cdot q_i)\cdot {\left(\frac{q_i^T\cdot y}{q_i^T\cdot q_i}\right)}^2$

Simplifying this expression yields the exact energy contribution formula:

${\text{Contribution}}_i=\frac{(q_i^T\cdot y{)}^2}{q_i^T\cdot q_i}$

The Orthogonal Least Squares (OLS) algorithm uses this derived metric as its core heuristic for forward selection:

• Selection: At each step, evaluate the contribution formula for all remaining candidate test points. Select the test point that yields the absolute highest value.

• Stopping Criteria: The algorithm iteratively adds centers until one of two conditions is met:

  1. A pre-determined maximum number of centers ($m$) is reached.

  2. The cumulative explained variance provides a “reasonable approximation” of the total system variance ($||y|{|}^2$).

Because the formula relies on the vectors ($q_i$) being orthogonal, adding a new center fundamentally changes the system. When a test point is selected as a center, the Gram-Schmidt algorithm must be applied to orthogonalize all remaining candidate vectors against the newly chosen center. This ensures that the contribution formula only ever measures new, unexplained variance in subsequent steps.