Model selection is the process of adjusting the parameters (and hyper-parameters) of a prediction model so that, after training, it will perform well on future samples. This adjustment is generally guided by the performance of the model on a test set, either once-off, or in a cross-validation fashion.
In this blog, we previously compared approaches based on cross-validation (working with random selection), and on deterministic splitting using the Kennard-Stone algorithm.
Deterministic methods have some advantage, especially if the size of the dataset is limited, in that they can be designed to maximise coverage of the feature space, hence the robustness of the resulting model.
The SPXY algorithm is a variant of the Kennard-Stone algorithm, designed to take into account both spectra (predictors) and response (i.e. the quantity to be predicted from the spectra). If we denote the spectra array by X and the response array by Y, SPXY stands for “Sample set Partitioning based on joint X-Y distances”. The algorithm was first proposed by Galvao et al., Talanta 67, 736 (2005).
Mathematical introduction to the SPXY algorithm
Before working through the equation, it’s useful to review the Kennard-Stone (KS) splitting algorithm, which forms the basis of the SPXY algorithm. You can also check our previous post on the topic, for more info and code implementation.
KS starts by calculating the Euclidean distance between each pair of spectra. Assume that we have Msamples, each denoted by an index m = 1, ..., M and N wavelengths, each denoted by an index n=1, ..., N. Then, the Euclidean distance between samples i and j is
d_{ij} = \sqrt{\sum_{n = 1}^{N} (x_{in} - x_{jn})^2},where the sum is done on the wavelength index.
The second step is to find the most distant pair of spectra, that is the pair that maximises the Euclidean distance calculated above.
The third (iterative) step is to find all other samples for the training set, by iteratively selecting the remaining samples having the greatest distance from the spectra already selected. The iteration stops when a predetermined number of samples have been selected.
The SPXY algorithm follows exactly the same logical steps, but it takes into account the Euclidean distances in both X and Y space. Let’s rename the distance in the formula above as d^{X}_{ij}, to make it clear that is the distance between the spectra i and j. Assume that the response values associated to these spectra are y_{i} and y_j respectively. Then the corresponding Euclidean distance of the response is simply
d^{Y}_{ij} = \vert y_{i} - y_{j} \vert,since each response is single valued. In the most general form of multi-valued response, where each measured response has K components, the formula for the Euclidean distance is
d^{Y}_{ij} = \sqrt{\sum_{k = 1}^{K} (\tilde{y}_{ik} - \tilde{y}_{jk})^2},where I put a “tilde” symbol (~) on the y to indicate that the different components must by scaled (for instance with a standard scaler) before calculating the Euclidean distance, since the value and range of these responses will be different, in general.
The last step of the SPXY is to sum the normalised version of d^{X}_{ij} and d^{Y}_{ij}:
d_{ij} = d^{X}_{ij} / d^{X}_{\mathrm{max}} + d^{Y}_{ij} / d^{Y}_{\mathrm{max}} .The normalisation is done by dividing each pairwise distance by the maximum distances d^{X}_{\mathrm{max}} and d^{Y}_{\mathrm{max}} respectively, ensuring that spectra and response contributions are assigned equal weights. Of course, you can imagine a generic weighted sum
d_{ij} = \alpha_{} d^{X}_{ij} / d^{X}_{\mathrm{max}} + (1 - \alpha) d^{Y}_{ij} / d^{Y}_{\mathrm{max}} ,whereby the relative weight between spectra and response is tuned by a factor \alpha bounded between 0 and 1. In the case in which \alpha = 1 , we recover the conventional KS algorithm, and when \alpha = 0.5 , we have (up to a constant) the SPXY algorithm.
Calculating pairwise distances
The major step towards a practical implementation of the SPXY algorithm, is the efficient calculation of pairwise distances. For this task, we use the combination of two functions available from SciPy: pdist and squareform .
Let’s start with importing the required libraries and loading a dataset
|
1 2 3 4 5 6 7 8 9 |
import pandas as pd import numpy as np from scipy.spatial.distance import pdist, squareform url = 'https://raw.githubusercontent.com/nevernervous78/nirpyresearch/master/data/Incombustible_Material_NIR.csv' data = pd.read_csv(url) X = data.values[:,2:] y = data["TIC Value"].values |
For the record, the data comprises 62 NIR spectra of mixes of coal and calcium carbonate, with varying “total incombustible content” (TIC). The data is raw and unmodified, and was collected using a ChemWiz-ADK NIR spectrometer (data courtesy of S. D’Hyon).
Now let’s calculate the pairwise distances with a straightforward application of the pdist function
|
1 2 |
pairwise_dist = pdist(X, metric='euclidean') print(pairwise_dist.shape) |
The function returns a 1D array with 1891 elements. This is equal to the number of pairwise combinations of M=62 samples, which is
M (M-1)/2.
We now apply the squareform function to pairwise_dist :
|
1 2 |
dist_x = squareform(pairwise_dist) print(m.shape) |
This transform the initial 1D array into a MxM matrix (2D array), which in our case, is a 62×62 elements array. The distance matrix is easier to interpret, since it shows the distances in a grid. For example, let’s print the first 5 elements with the command print(dist_x[:5, :5]) . The result is
|
1 2 3 4 5 |
[[0. 4.60603004 5.84176011 3.03871023 2.92160501] [4.60603004 0. 5.32967143 3.45723176 5.4538194 ] [5.84176011 5.32967143 0. 5.92814338 6.51508812] [3.03871023 3.45723176 5.92814338 0. 3.49743732] [2.92160501 5.4538194 6.51508812 3.49743732 0. ]] |
The top-left elements is 0, which is the distance of the first spectrum from itself. Moving one step to the right, we get the distance of the first spectrum from the second, and so on. Since the order in which the two spectra are taken doesn’t matter in the Euclidean distance, dist_x is a symmetric matrix.
Wrapping this all up, the distance metric for the SPXY algorithm is calculated as follows.
|
1 2 3 4 5 6 7 8 9 10 11 12 |
# Calculate pairwise distances for the spectra dist_x = squareform(pdist(X, metric='euclidean')) # Calculate pairwise distances for the responses dist_y = squareform(pdist(y.reshape(-1, 1), metric='euclidean')) # Normalise the results dist_x = dist_x / np.max(dist_x) dist_y = dist_y / np.max(dist_y) # Calculate XY distances dist_xy = dist_x + dist_y |
Note that the y array is originally 1D, therefore it must be reshaped to 2D form to be a correct input to pdist .
A home-made implementation of the SPXY algorithm
Now that we have all the pairwise distances, we can start the SPXY selection process. We begin by setting up the required lists and finding the indices corresponding to the maximum pairwise distance
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
# Test size must be whitin 0 and 1 test_size = 0.25 # Get the number of training and test samples n_samples = X.shape[0] n_test = int(n_samples * test_size) n_train = n_samples - n_test # Initialize indices train_idx = [] # Training set. Initially empty remaining_idx = list(range(n_samples)) # The remaining indices. Initially, all of them. # Find the indices of the pair with maximum distance max_idx = np.unravel_index(np.argmax(dist_xy), dist_xy.shape) # Add the indices to the training set list train_idx.extend([max_idx[0], max_idx[1]]) # Remove these indices from the list of the remaining spectra remaining_idx.remove(max_idx[0]) remaining_idx.remove(max_idx[1]) |
Now we have the first two elements of the training set. Next with need to calculate the distance of each remaining spectra from this pair of elements. By ‘distance’, in this case, we mean the distance to the closest element of the pair. In other word, for each remaining spectra we need to 1) Calculate the two distances to the selected pair, 2) take the smallest one and store it, and 3) repeat for all the remaining spectra. Finally, we’ll have to select the maximum of all the minimum distances calculated previously
(Similarly, when we will have selected k spectra, the ‘distance’ of a remaining spectrum from the set of k spectra, is the distance of such a spectrum from the closest element. )
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
# Iteratively select samples while len(train_idx) < n_train: # Calculate minimum distance of each remaining spectrum to each spectrum in the training set min_distances = [] for idx in remaining_idx: distances_to_train = [dist_xy[idx, train_pt] for train_pt in train_idx] min_distances.append(min(distances_to_train)) # Select the element with maximum minimum-distance selected = remaining_idx[np.argmax(min_distances)] train_idx.append(selected) remaining_idx.remove(selected) # Remaining samples are assigned to the test set test_idx = remaining_idx |
Done! Now that we have the indices, the training-test split can be done by selecting the corresponding elements.
|
1 2 |
X_train, X_test = X[train_idx], X[test_idx] y_train, y_test = y[train_idx], y[test_idx] |
An existing (and fast) implementation of the SPXY algorithm
For an existing implementation of the SPXY algorithm (as well as the KS and a host of other algorithms), I can heartily recommend the library astartes, developed by Jackson Burns and coworkers. This library was developed as a drop-in replacement for the test_train_split routine of scikit-learn.
import astartes as ast X_train, X_test, y_train, y_test, train_idx, test_idx = ast.train_test_split(X,y,sampler = ‘spxy’, return_indices=True)
By changing the sampler keyword, you can apply different algorithms. For a complete list of the available samplers, check the Astartes guide.
References
[1] R. Galvão et al., A method for calibration and validation subset partitioning, Talanta, 67, 736-740 (2005).
[2] See also the open-access article: Y. Xu and R. Goodacre, On Splitting Training and Validation Set: A Comparative Study of Cross-Validation, Bootstrap and Systematic Sampling for Estimating the Generalization Performance of Supervised Learning, Journal of Analysis and Testing, 2, 249-262 (2018).
**
The feature image is a posterized rendering of an original photo by the author.
