Machine Learning

• Observed data \(\mathbf{x}\), called input variables or predictors.
• Data to be predicted \(y\), called output variables or responses.

Usually, the input variables is composed of \(d\) features represented as a \(d\)-dimensional vector of vector space \(\mathcal{X}\): \(\mathbf{x}=[x_1, \ldots,x_d]^\top\in\mathcal{X}\subset\mathbb{R}^d\). For simplicity, in the following, we will refer to \(\mathbb{R}^d\), but be aware that for some situations, we might work with a subspace of it. For instance, we can have strictly positive measurements, i.e., \(x_i \geq 0, \forall i\in{1,\ldots,d}\), and this kind of constraint could be included in the learning problem.

• If \(y\) are quantitative data, the learning problem is usually called regression and \(y\in\mathcal{Y}\subset\mathbb{R}\).
• If \(y\) are categorical data, the learning problem is usually called classification and \(y\in\{C_1, \ldots, C_C\}\) with \(C_i\) refers to the \(i^{th}\) category.

These are very close problem: the algorithm might differ only by the performance measure used to learn the predictive model.

The link between the input and output variables is modeled by the prediction rule, denoted by the function \(f\). This function takes as input the predictors and outputs response … Mathematically, it writes:

The (learnable) parameters of the function are denoted \(\boldsymbol{\theta}\) and are optimized (or fit) with the learning set, i.e., they are automatically tuned to increase the performance measure. Once the optimal set of parameters are obtained, the function is usually denoted \(\hat{f}\), indicating that it was fitted on the data. We might also define hyperparameters of \(f\) as the set parametes of the function that is fixed before the learning set and are not optimized.

We call training set \(\mathcal{S}\) the set of sample available to learn the prediction rule \(f\). For a training set of size \(n\), three different cases exist:

1. Supervised learning: \(\mathcal{S} = \{(\mathbf{x}_1, y_1), \ldots, (\mathbf{x}_n, y_n)\}\), both input and output variables are available.
2. Unsupervised learning: \(\mathcal{S} = \{(\mathbf{x}_1), \ldots, (\mathbf{x}_n)\}\), only input variables are available.
3. Semi-supervised learning: \(\mathcal{S} = \mathcal{S}_s \cup \mathcal{S}_u = \{(\mathbf{x}^s_1, y^s_1), \ldots, (\mathbf{x}^s_{n_s}, y^s_{n_s})\} \cup \{\mathbf{x}^u_1,\ldots, \mathbf{x}^u_{n_u}\}\), this is a mixed scenario often encountered in practice. Some samples are provided with the output variables, but most of the sample are available only with input variable.

The table tab:braking:distance contains the braking distance for various cars and driver. The file can be download here. Load the data and try a linear regression (with different polynomial degree) to related the speed with the braking distance.

Tips: you can load the file using the following snippet:

import pandas as pd
df = pd.read_csv("/home/mfauvel/Documents/Enseignement/INRA/MasterSigma/data/braking_distance.csv", index_col=0) # Path need to be adapted !
print(df.head())

       Speed: 0  Speed: 40  Speed: 50  Speed: 70  Speed: 90  Speed: 100  Speed: 110  Speed: 130
Car #                                                                                          
0      4.488990   12.60200    7.93177    32.4445    60.5546     56.0518     76.0007     101.365
1     -0.561011   11.11080   11.79650    35.3948    57.9588     63.6805     78.7270     107.571
2     -3.026890    1.03175   25.18210    28.4790    48.5566     61.8898     76.6653     104.722
3     -6.701320   11.58230    8.69358    25.5270    49.7054     59.2811     70.4830     107.165
4     -5.491230   11.01720    8.43203    26.6295    52.0120     66.0341     81.1067     115.733

As discussed in 2.3, if we fit a linear model with the mean square error, we have an explicit solution for 1-dimensional data. The result is still valid for any dimension. To see this, denote \(\mathbf{X} = \big[\mathbf{x}_1, \ldots, \mathbf{x}_n\big]^\top\) a \(n\times d+1\) matrix that collects all the input samples augmented with one variable set to 1 (for the intercept), and let \(\mathbf{y}=[y_1, \ldots, y_n]^\top\) the vector containing all the outputs samples. Then the MSE can be written as

where \(\mathbf{w}= [w_1, \ldots, w_d, b]^\top\). We can derive w.r.t \(\mathbf{w}\) and by setting the derivative to zero, we obtain the general expression for fitting a linear model with the MSE:

This is exactly what the fit method does in scikit-learn. The “hat” on \(\mathbf{w}\) indicates that it is an estimation, w.r.t. collected data.

Interpretability of a linear model comes from the possibility to relate the values of \(\hat{\mathbf{w}}\) to the output: a very low value of \(w_i\) indicates a very low relationship¹ of the i\(^{th}\) input variable with the output variable, providing input variables have been standardized (see section 5).

In the following, we will apply the linear regression to the California Housing dataset from the scikit-learn datasets (https://scikit-learn.org/stable/datasets/real_world.html#california-housing-dataset). From the data set description:

the output variable is the median house value for California districts, expressed in hundreds of thousands of dollars ($100,000). The input variables are

• MedInc median income in block group.
• HouseAge median house age in block group.
• AveRooms average number of rooms per household.
• AveBedrms average number of bedrooms per household.
• Population block group population.
• AveOccup average number of household members.
• Latitude block group latitude.
• Longitude block group longitude.

where a block group is the smallest geographical unit for which the U.S. Census Bureau publishes sample data (a block group typically has a population of 600 to 3,000 people).

import pandas as pd
import geopandas as gpd
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# Open data and create the geo-data-frame
cal_housing = fetch_california_housing(as_frame=True)
df = pd.concat([cal_housing.data, cal_housing.target], axis=1)
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.Longitude, df.Latitude))


# Inspect the first 5 rows of the dataframe
print(gdf.head())

   MedInc  HouseAge  AveRooms  AveBedrms  Population  AveOccup  Latitude  Longitude  MedHouseVal                     geometry
0  8.3252      41.0  6.984127   1.023810       322.0  2.555556     37.88    -122.23        4.526  POINT (-122.23000 37.88000)
1  8.3014      21.0  6.238137   0.971880      2401.0  2.109842     37.86    -122.22        3.585  POINT (-122.22000 37.86000)
2  7.2574      52.0  8.288136   1.073446       496.0  2.802260     37.85    -122.24        3.521  POINT (-122.24000 37.85000)
3  5.6431      52.0  5.817352   1.073059       558.0  2.547945     37.85    -122.25        3.413  POINT (-122.25000 37.85000)
4  3.8462      52.0  6.281853   1.081081       565.0  2.181467     37.85    -122.25        3.422  POINT (-122.25000 37.85000)

We can use pandas info() function to get some information about the data:

# Print some information about the data set_model_parameters
gdf.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 20640 entries, 0 to 20639
Data columns (total 10 columns):
 #   Column       Non-Null Count  Dtype   
---  ------       --------------  -----   
 0   MedInc       20640 non-null  float64 
 1   HouseAge     20640 non-null  float64 
 2   AveRooms     20640 non-null  float64 
 3   AveBedrms    20640 non-null  float64 
 4   Population   20640 non-null  float64 
 5   AveOccup     20640 non-null  float64 
 6   Latitude     20640 non-null  float64 
 7   Longitude    20640 non-null  float64 
 8   MedHouseVal  20640 non-null  float64 
 9   geometry     20640 non-null  geometry
dtypes: float64(9), geometry(1)
memory usage: 1.6 MB

The very first steps of any machine learning analysis is to split the data into two sets: One used to train/learn the model, and the second one used to compute the accuracy of the model. Scikit-learn propose several tools to do such operations (see https://scikit-learn.org/stable/modules/classes.html#module-sklearn.model_selection). In this course, we will make use a lot of the train_test_split function:

states = gpd.read_file("data/usa-states-census-2014.shp")
# Split the data into train and test data set
gdf_train, gdf_test = train_test_split(gdf, train_size=0.5, random_state=10)

# Plot the target value
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(15, 5), edgecolor='w')
for i, (gdf, data) in enumerate(zip([gdf_train, gdf_test], ["train", "test"])):
    states[states.NAME == "California"].plot(color="gray", alpha=0.5, ax=axs[i])
    gdf.plot(
        kind="scatter",
        x="Longitude",
        y="Latitude",
        s=gdf["Population"] / 100,
        c="MedHouseVal",
        colorbar=True,
        ax=axs[i],
    )
    axs[i].set_title(f"{data}")

For classification problems, we have a categorical variable to estimate instead of a having quantitative variable, it means the algorithm should predict categories (“Dogs”, “Cats”, “Monkey” …). For such output variable, the MSE is not a good measure of performance: how do you define the square error between category “Dogs” and “Cats” ?

One solution usually adopted for linear models is to consider the estimation of the posterior probability \(p(y_i=c | \mathbf{x}_i)\) to each class \(c\in\big\{1, \ldots, C\big\}\) and decide the class membership using the maximum a posteriori (MAP) rule: \[\hat{c}_i = \arg\max_{c} p(y_i=c | \mathbf{x}_i).\]

The logistic regression (despite its name it is a classifier algorithm) models the logarithm of the posterior probability as a linear function of the inputs variables (Hastie, Tibshirani, and Friedman 2009): \[\log\big(p(y=c | \mathbf{x})\big) = \mathbf{w}_c^\top\mathbf{x} + b_c + K\] where \(K\) is a constant independent of the class. You should recognise the expression of a linear model of regression: this is exactly what it is. The logistic regression performs a regression on the logarithm of the posterior probability of each class: we have to learn \(C\) regression problems. However, these problems are not independent since the sum of the posterior probabilities should be exactly 1: this is the role of the constant \(K\):

thus

\[p(y_i=c | \mathbf{x}_i) = \frac{\exp\Big\{\mathbf{w}_{c}^\top\mathbf{x}_i + b_{c}\Big\}}{\sum_{c'=1}^C \exp\Big\{\mathbf{w}_{c'}^\top\mathbf{x}_i + b_{c'}\Big\}}.\]

This ratio of exponential is known as the softmax function and it is used (widely in neural network) to transform a set of unconstrained variables to a set of variables between \([0, 1]\) that sums to one (as probability distribution).

To train our logistic regression model, we need to estimate \((\mathbf{w}_c, b_c)\ \forall c\in\{1, \ldots, C\}\). From our training set, we have the following information:

The output variable \(y_i\) is vector of size the number of class, and contains the information above, usually on-hot-encoded (https://machinelearningmastery.com/why-one-hot-encode-data-in-machine-learning/): \[\mathbf{y}_i=\begin{bmatrix}0\\ \vdots\\1\\ \vdots\\0\end{bmatrix} \text{ with }\mathbf{y}_{ij}=1 \text{ if } \mathbf{x}_i\in c_j \text{ and } 0 \text{ otherwise}.\]

For such observations (membership probabilities), the loss conventionally used is the negative of the log likelihood, which can be written as (the regularization is already included):

There is no explicit solution and the optimization problem is solved by gradient descent. For linear model, there are many existing algorithms, with different regularization and optimization strategies. The scikit-learn library has implemented a lot of them, you can read a description here: https://scikit-learn.org/stable/modules/linear_model.html#logistic-regression.

We will illustrate how the logistic regression works on a synthetic example.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LogisticRegression as LR
from sklearn.metrics import accuracy_score as  OA

# Create toy 2D data 
X, y = make_blobs(n_samples=100,
                  centers=3,
                  n_features=2,
                  random_state=0)

# Split
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.5, random_state=0)

# Fit model
model = LR(random_state=0, C = 1)
model.fit(X_train, y_train)
y_pred_test = model.predict(X_test)

That’s all, as usual with scikit learn. So far, we have generated some toy data set with 2 input variables and 3 classes. The generated data have been split into a training and testing set and we learn the model with default parameters.

Now since we are in two dimensional space, it is easy to plot the prediction of the model. We can output two information: the posterior probability (the class membership) and the classification. Sometimes they are referred to as soft and hard assignment values. Scikit learn offers two function: predict_proba and predict for the former and the latter respectively. The following code generates a grid of points in \(\mathbb{R}^2\) and compute their class membership, then it displays it using contours plot.

# Generate points for the bounding box
x_min, x_max = X[:, 0].min() - 0.5, X[:, 0].max() + 0.5
y_min, y_max = X[:, 1].min() - 0.5, X[:, 1].max() + 0.5
h = 0.02  # step size in the mesh
xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))

# Predict posterior probability and class
Z = model.predict_proba(np.c_[xx.ravel(), yy.ravel()])
y_pred = model.predict(np.c_[xx.ravel(), yy.ravel()])

# Plot the results
fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15, 5))
for c in range(3):
    cs = axs[c].contour(xx, yy, Z[:,c].reshape(xx.shape), [0.25, 0.5, 0.75])
    axs[c].scatter(X_test[:, 0],
                        X_test[:, 1],
                        c=y_test)
    axs[c].clabel(cs, inline=1, fontsize=10)

axs[-1].contourf(xx, yy, y_pred.reshape(xx.shape))
axs[-1].scatter(X_test[:, 0],
                X_test[:, 1],
                c=y_test)

The figure 9 shows the classification result: the decision boundaries are linear w.r.t. the input variables (on the right plot, boundaries are straight lines). It is off course possible to perform polynomial transformation before the linear model to obtain non-linear decision function.

A decision tree is a (simple & naive) model that partition the input variables space into zones where the prediction (regression or classification) is constant. The learning step consists in finding these zones. Basically, the algorithm recursively search for each input variable a threshold that satisfy a purity index, i.e., that splits the data into into pure subset (in terms of output variables), and construct a set of if ... then ... else binary rules that defines these zones. Since each zone contains “pure” output variables, i.e. very similar values (same class number for instance) a constant prediction make sense.

Let’s us continue with the toy regression example to better understand how a decision tree behave. The figure 11 shows a decision tree fits with two size of tree, here controlled by the max_depth parameter.

import numpy as np
from sklearn.tree import DecisionTreeRegressor
import matplotlib.pyplot as plt
from sklearn import tree

# Load data
data = np.load("./data/toy_data.npz")
x, y, xt, yt = data["x"], data["y"], data["xt"], data["yt"]

# Fit model & plot model
plt.figure(figsize=(5,5))
plt.scatter(x, y, c='b')
plt.scatter(xt[::10,0], yt[::10,0], c='k')
for max_depth, col in zip([10, 2],
                          ['r', 'g']):
    model = DecisionTreeRegressor(max_depth=max_depth)
    model.fit(x, y)

    # Plot decision function
    t = np.linspace(-5, 5, 1000).reshape(-1, 1)
    plt.plot(t, model.predict(t), col)
plt.legend(["Train points","Test points","Depth=10", "Depth=3"])

From the figure, it can be seen that a too deep tree might overfit the data. We will see in section 4 how to select it. Yet, in the next section, we will show that pragmatically, it is not of high importance.

Sklearn provides function to plot decision tree: this allow to interpret how the decision is done and which variables are the most important. This this the plot_tree function (https://scikit-learn.org/stable/modules/generated/sklearn.tree.plot_tree.html#sklearn.tree.plot_tree). The figure 12 shows a decition tree with a max_depth of size 2.

A decision tree can be either used for classification or regression.

Decision tree can overfit and its learning capacity is limited (e.g., because of constant prediction). However, by combining many trees, build on a subset of the data and a subset of the input variables it is possible to achieve much better results. This is the idea of Random Forest (RF), a powerful ensemble method (Hastie, Tibshirani, and Friedman 2009).

The algorithm of a random forest can be summarized as:

1. Select the number of trees and the parameters of each tree (max_depth etc ..)
2. For each tree
   • Select a subset of the sample
   • Fit the decision tree recursively by selecting a subset of input variables for each node
3. For the prediction of one sample:
   • Compute the decision for each tree
   • Average the results for regression or do a majority vote for classification.

RF is usually fast and little sensitive to hyperparameter setting: if the number of trees is high enough, good results will be achieved without the need to investigate several values.

The combination of several decision trees does not maintain the interpretability of decision tree. Yet, RF comes with the ability to rank input features according to how they improved the purity criterion at each node of the decision trees. In sklearn, this is accessible using the feature_importance_ property of the model.

As usual, the behavior of RF is illustrated on the toy example. Here for classification.

from sklearn.ensemble import RandomForestClassifier as RF

# Learn
model = RF()
model.fit(X_train, y_train)
y_pred_test = model.predict(X_test)
Z = model.predict_proba(np.c_[xx.ravel(), yy.ravel()])
y_pred = model.predict(np.c_[xx.ravel(), yy.ravel()])

# Plot the results
fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15, 5))
for c in range(3):
    cs = axs[c].contour(xx, yy, Z[:,c].reshape(xx.shape), [0.25, 0.5, 0.75])
    axs[c].scatter(X_test[:, 0],
                        X_test[:, 1],
                        c=y_test)
    axs[c].clabel(cs, inline=1, fontsize=10)

axs[-1].contourf(xx, yy, y_pred.reshape(xx.shape))
axs[-1].scatter(X_test[:, 0],
                X_test[:, 1],
                c=y_test)

One of the most simple algorithm, yet effective but time consuming, is the sequential forward feature selection. It can be applied to any supervised regression/classification algorithm. It works as follows:

1. Starts with an empty pool \(F\) of selected features,
2. Select the feature \(f_1\) that provides the best value for the selected criterion and add it to \(F\).
3. Select the feature \(f_2\) such that the couple of features \((f_1,f_2)\) provides the best value for the selected criterion and add it to \(F\).
4. Select the feature \(f_3\) such that the triplet of features \((f_1,f_2,f_3)\) …
5. …
6. The algorithm stops either if the increase of the criterion is too low or if the maximum number of features is reached.

A variant of the algorithm is obtained by sequentially remove one feature from the set of input feature, no surprise it is called sequential backward feature selection. In practice this version if less effective than the forward one.

Let’s see how it works for the California housing data set.

import numpy as np
from sklearn.datasets import fetch_california_housing
from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt

# Open data
california_housing_data = fetch_california_housing()
X, y = california_housing_data.data, california_housing_data.target
feature_names = np.asarray(california_housing_data.feature_names)

# Split data
X_train, X_test, y_train, y_test = train_test_split(X, y)

We now use the class SequentialFeatureSelector from Scikit-learn (https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.SequentialFeatureSelector.html#sklearn.feature_selection.SequentialFeatureSelector). We need to also set a regression model. For this example a Knn is used, but any models from scikit is working.

from sklearn.feature_selection import SequentialFeatureSelector as  SFS
from sklearn.neighbors import KNeighborsRegressor as KNR

knn = KNR(n_neighbors=3)
sfs = SFS(knn, n_features_to_select=3, n_jobs=-1)

sfs.fit(X_train, y_train)
print(f"Selected feature: {feature_names[sfs.get_support()]}")

Selected feature: ['MedInc' 'Latitude' 'Longitude']

We can compute the score on the test set and compare it when no feature selection is done.

knn.fit(sfs.transform(X_train), y_train)
print(f"Regression score on the test set with feature selection: {knn.score(sfs.transform(X_test), y_test)}")
knn.fit(X_train, y_train)
print(f"Regression score on the test set without feature selection: {knn.score(X_test, y_test)}")

Regression score on the test set with feature selection: 0.7539353810565258
Regression score on the test set without feature selection: 0.11211148340676591

In this part, we will cover the famous Principal Component Analysis (PCA), an unsupervised feature extraction algorithm.

PCA is a Linear transformation used to reduce the dimensionality of the data (Jolliffe 2002). The new feature \(z_i\) is found by projecting the original data (i.e., a linear combination of), such as the projected feature account for most of the variability of the data:

• \(z_1,~z_2,~z_3,\ldots\) are mutually uncorrelated (i.e., \(\langle z_i, z_j\rangle=0, \forall j\neq i\)),
• \(\text{var}(z_i)\) is as large as possible,
• \(\text{var}(z_1)>\text{var}(z_2)>\text{var}(z_3)>\ldots\)

The PCA algorithm is as follow:

• It starts by searching \(\mathbf{v}_1\) such as \(\max_{\mathbf{v}_1}\text{var}(z_1)\) with

  \begin{eqnarray}\label{eq:lagrangian} \text{var}(z_1) & = & \text{var}(\langle\mathbf{v}_1,\mathbf{x}\rangle)\\ &=& \mathbf{v}_1^\top\boldsymbol{\Sigma}\mathbf{v}_1 \end{eqnarray}

  This problem is indetermined: if \(\hat{\mathbf{v}}_1\) maximizes the variance, \(\alpha\hat{\mathbf{v}}_1\) too! Therefore, we add a constraint on \(\mathbf{v}_1\): \(\langle\mathbf{v}_1,\mathbf{v}_1\rangle=1\). Thus, we have a constraint optimization problem, whose Lagrangian can be written as  (Boyd and Vandenberghe 2004):

  \begin{eqnarray} \mathcal{L}(\mathbf{v}_1,\lambda_1) = \mathbf{v}_1^\top\boldsymbol{\Sigma}\mathbf{v}_1 + \lambda_1(1- \mathbf{v}_1^\top\mathbf{v}_1) \end{eqnarray}

  By computing the derivative w.r.t \(\mathbf{v}_1\) and set it to zeros, we have:

  \begin{eqnarray} \frac{\partial\mathcal{L}}{\partial\mathbf{v}_1} = 2\boldsymbol{\Sigma}\mathbf{v}_1-2\lambda_1\mathbf{v}_1 \end{eqnarray}

  Hence \(\mathbf{v}_1\) is an eigenvector of the covariance matrix of \(\mathbf{x}\):

  \begin{eqnarray} \boldsymbol{\Sigma}\mathbf{v}_1 =\lambda_1\mathbf{v}_1 \end{eqnarray}

  To get a maximal variance, \(\mathbf{v}_1\) should be the eigenvector corresponding to the largest eigenvalues (to see this, we insert the above equation into \ref{eq:lagrangian})!

  \begin{eqnarray} \text{var}(z_1) = \mathbf{v}_1^\top\boldsymbol{\Sigma}\mathbf{v}_1 = \lambda_1 \mathbf{v}_1^\top\mathbf{v}_1 = \lambda_1 \end{eqnarray}

  Any linear algebra library has a function to find the eigenvectors/eigenvalues of a square positive definite matrix (see for instance https://numpy.org/doc/stable/reference/generated/numpy.linalg.eigh.html)

• Once \(\mathbf{v}_1\) if found, the algorithm no searches for \(\mathbf{v}_2\) such as \(\max\text{var}(z_2)\) and \(\langle\mathbf{v}_2,\mathbf{v}_2\rangle=1\) and \(\langle\mathbf{v}_1,\mathbf{v}_2\rangle=0\). Again we can write the Lagragian:

  \begin{eqnarray} \mathcal{L}(\mathbf{v}_2,\lambda_2,\beta_1) = \mathbf{v}_2^\top\boldsymbol{\Sigma}\mathbf{v}_2 + \lambda_2(1- \mathbf{v}_2^\top\mathbf{v}_2) + \beta_1(0 - \mathbf{v}_2^\top\mathbf{v}_1) \end{eqnarray}

  and compute the derivative w.r.t \(\mathbf{v}_2\), and try to vanish it:

  \begin{eqnarray} \frac{\partial\mathcal{L}}{\partial\mathbf{v}_2} &=& 2\boldsymbol{\Sigma}\mathbf{v}_2-2\lambda_2\mathbf{v}_2-\beta_1\mathbf{v}_1\\ \boldsymbol{\Sigma}\mathbf{v}_2 &=& \lambda_2\mathbf{v}_2+2\beta_1\mathbf{v}_1 \end{eqnarray}

  At optimality, \(\langle\mathbf{v}_1,\mathbf{v}_2\rangle=0\). Left-multiplying by \(\mathbf{v}_1^\top\) the above equation:

  \begin{eqnarray} \mathbf{v}_1^\top\boldsymbol{\Sigma}\mathbf{v}_2 &=& 2\beta_1 \\ \lambda_1\mathbf{v}_1^\top\mathbf{v}_2 &=& 2\beta_1 \\ 0 &=& 2\beta_1 \end{eqnarray}

  Hence, we have

  \begin{eqnarray} \boldsymbol{\Sigma}\mathbf{v}_2 =\lambda_2\mathbf{v}_2. \end{eqnarray}

  \(\mathbf{v}_2\) is the eigenvector corresponding the second largest eigenvalues. For axis \(k\ge 3\), we can repeat the same procedure and get that \(\mathbf{v}_k\) is the eigenvector corresponding the \(k^{\text{th}}\) largest eigenvalues.

In practice, the eigensolver returns all the eigenvector/eigenvalues in one call. Basically, the algorithm reduces to

1. Empirical estimation the mean value:

   \begin{eqnarray} \boldsymbol{\mu} = \frac{1}{n}\sum_{i=1}^n\mathbf{x}_i \end{eqnarray}

2. Empirical estimation the covariance matrix:

   \begin{eqnarray} \boldsymbol{\Sigma} = \frac{1}{n-1}\sum_{i=1}^n(\mathbf{x}_i-\boldsymbol{\mu})(\mathbf{x}_i-\boldsymbol{\mu})^\top \end{eqnarray}

3. Compute the eigenvalues/eigenvectors and keep the \(p\) first… How to choose \(p\) ? Take a look to the explained variance, i.e., the variance of the data explained by the \(p\) first eigenvectors, it can be computed as

   \[\frac{\sum_{i=1}^p\lambda_i}{\sum_{i=1}^d\lambda_i}\]

   In practice, we choose \(p\) such as a given amount of the variance is kept (e.g., 90% or 95%).

4. Tips for high dimensional data set: if the number of samples is lower than the number of feature, see (Manolakis, Lockwood, and Cooley 2016) page 420.

We now compare how PCA works on the California housing data (see https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html for a list of options):

from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

# Two cases are investigated: same number of feature than with sequential feature selector
# or by thresholing the cummulative variance
for param in [3, 0.90, 0.95, 0.99, 0.999]:
    model = Pipeline([
        ("Scaler", StandardScaler()), # You can try to remove this part and see the difference
        ("Feature Extraction", PCA(n_components = param)),
        ("KNN", KNR(n_neighbors=3))
    ])
    model.fit(X_train, y_train)
    print(f"Regression score with the number of component {model.named_steps['Feature Extraction'].n_components_}: {model.score(X_test, y_test)}")

Regression score with the number of component 3: 0.28924002334403465
Regression score with the number of component 5: 0.47881899995519794
Regression score with the number of component 6: 0.523310120515489
Regression score with the number of component 7: 0.6129552696601799
Regression score with the number of component 8: 0.6368428619732054