Régression linéaire par morceaux#

Le paragraphe Régression linéaire étudie le lien entre le coefficient R^2 et la corrélation pour finalement illustrer une façon de réaliser une régression linéaire par morceaux. L’algorithme s’appuie sur un arbre de régression pour découper en morceaux ce qui n’est pas le plus satisfaisant car l’arbre cherche à découper en segment en approximant la variable à régresser Y par une constante sur chaque morceaux et non une droite. On peut se poser la question de comment faire pour construire un algorithme qui découpe en approximant Y par une droite et non une constante. Le plus dur n’est pas de le faire mais de le faire efficacement. Et pour comprendre là où je veux vous emmener, il faudra un peu de mathématiques.

Une implémentation de ce type de méthode est proposée dans la pull request Model trees (M5P and co) qui répond à au problème posée dans Model trees (M5P) et originellement implémentée dans Building Model Trees. Cette dernière implémentation réestime les modèles comme l’implémentation décrite au paragraphe Implémentation naïve d’une régression linéaire par morceaux mais étendue à tout type de modèle.

Exploration#

Problème et regréssion linéaire dans un espace à une dimension#

Tout d’abord, une petite illustration du problème avec la classe PiecewiseRegression implémentée selon l’API de scikit-learn.

../_images/piecenaive.png

Cette régression par morceaux est obtenue grâce à un arbre de décision. Celui-ci trie le nuage de points (X_i, Y_i) par ordre croissant selon les X, soit X_i \leqslant X_{i+1}. L’arbre coupe en deux lorsque la différence des erreurs quadratiques est maximale, erreur quadratique obtenue en approximant Y par sa moyenne sur l’intervalle considéré. On note l’erreur quadratique :

\begin{array}{rcl} C_(i,j) &=& \frac{1}{j - i + 1} \sum_{i \leqslant k \leqslant j} Y_i \\ D_(i,j) &=& \frac{1}{j - i + 1} \sum_{i \leqslant k \leqslant j} Y^2_i \\ E_(i,j) &=& \frac{1}{j - i + 1} \sum_{i \leqslant k \leqslant j} ( Y_i - C(i,j))^2 = \frac{1}{j - i + 1} \sum_{i \leqslant k \leqslant j} Y_i^2 - C(i,j)^2 = D(i,j) - C(i,j)^2 \end{array}

La dernière ligne applique la formule \var{X} = \esp{X^2} - \esp{X}^2 qui est facile à redémontrer. L’algorithme de l’arbre de décision coupe un intervalle en deux et détermine l’indice k qui minimise la différence :

\Delta_k = E(1, n) - (E(1, k) + E(k+1, n))

L’arbre de décision optimise la construction d’une fonction en escalier qui représente au mieux le nuage de points, les traits verts sur le graphe suivant, alors qu’il faudrait choisir une erreur quadratique qui corresponde aux traits oranges.

../_images/piecenaive2.png

Il suffirait donc de remplacer l’erreur E par celle obtenue par une régression linéaire. Mais si c’était aussi simple, l’implémentation de sklearn.tree.DecisionTreeRegressor la proposerait. Alors pourquoi ? La raison principale est que cela coûte trop cher en temps de calcul. Pour trouver l’indice k, il faut calculer toutes les erreurs E(1,k) E(k+1,n), ce qui coûte très cher lorsque cette erreur est celle d’une régression linéaire parce qu’il est difficile de simplifier la différence :

\begin{array}{rcl} \Delta_{k} - \Delta_{k-1} &=& - (E(1, k) + E(k+1, n)) + (E(1, k-1) + E(k, n)) \\ &=& E(1, k-1) - E(1, k) + E(k, n) - E(k+1, n) \end{array}

Arbre de régression constante

On s’intéresse au terme E(1, k-1) - E(1, k) dans le cas le nuage de points est représenté par une constante sur chaque segment. C’est l’hypothèse faite par l’algorithme classique de construction d’un arbre de régression (segments verts sur le premier dessin) :

\begin{array}{rcl} C_(1,k-1) - C_(1,k) &=& \frac{1}{k-1} \sum_{1 \leqslant i \leqslant k-1} Y_i - \frac{1}{k} \sum_{1 \leqslant i \leqslant k} Y_i \\ &=& (\frac{1}{k-1} - \frac{1}{k}) \sum_{1 \leqslant i \leqslant k-1} Y_i - \frac{Y_k}{k} \\ &=& \frac{1}{k(k-1)} \sum_{1 \leqslant i \leqslant k-1} Y_i- \frac{Y_k}{k} \\ &=& \frac{1}{k} C(1,k-1) - \frac{Y_k}{k} \end{array}

On en déduit que :

\begin{array}{rcl} E(1, k-1) - E(1, k) &=& \frac{1}{k} D(1,k-1) - \frac{Y_k^2}{k} + (C_(1,k-1) - C_(1,k))(C_(1,k-1) + C_(1,k)) \\ &=& \frac{1}{k} D(1,k-1) - \frac{Y_k^2}{k} + \pa{\frac{1}{k} C(1,k-1) - \frac{Y_k}{k}} \pa{\frac{Y_k}{k} - \frac{1}{k} C(1,k-1) + 2 C(1,k-1)} \end{array}

On voit que cette formule ne fait intervenir que C(1,k-1), D(1,k-1), Y_k, elle est donc très rapide à calculer et c’est pour cela qu’apprendre un arbre de décision peut s’apprendre en un temps raisonnable. Cela repose sur la possibilité de calculer le critère optimisé par récurrence. On voit également que ces formules ne font pas intervenir X, elles sont donc généralisables au cas multidimensionnel. Il suffira de trier les couples (X_i, Y_i) selon chaque dimension et déterminer le meilleur seuil de coupure d’abord sur chacune des dimensions puis de prendre le meilleur de ces seuils sur toutes les dimensions. Le problème est résolu.

Le notebook Custom Criterion for DecisionTreeRegressor implémente une version pas efficace du critère MSE et compare la vitesse d’exécution avec l’implémentation de scikit-learn. Il implémente ensuite le calcul rapide de scikit-learn pour montrer qu’on obtient un temps comparable. Le résultat est sans équivoque. La version rapide n’implémente pas \Delta_{k} - \Delta_{k-1} mais plutôt les sommes \sum_1^k w_i Y_i, \sum_1^k w_i Y_i^2 dans un sens et dans l’autre. En gros, le code stocke les séries des numérateurs et des dénominateurs pour les diviser au dernier moment.

Arbre de régression linéaire

Le cas d’une régression est plus complexe. Prenons d’abord le cas où il n’y a qu’un seule dimension, il faut d’abord optimiser le problème :

E(1, n) = \min_{a,b} = \sum_{k=1}^n (a X_k + b - Y_k)^2

On dérive pour aboutir au système d’équations suivant :

\begin{array}{rcl} \frac{\partial E(1,n)}{\partial a} &=& 0 = \sum_{k=1}^n X_k(a X_k + b - Y_k) \\ \frac{\partial E(1,n)}{\partial b} &=& 0 = \sum_{k=1}^n a X_k + b - Y_k \end{array}

Ce qui aboutit à :

\begin{array}{rcl} a(1, n) &=& \frac{\sum_{k=1}^n X_kY_k - \pa{\sum_{k=1}^n X_k}\pa{\sum_{k=1}^n Y_k} } {\sum_{k=1}^n X_k^2 -\pa{\sum_{k=1}^n X_k}^2 } \\ b(1, n) &=& \sum_{k=1}^n Y_k - a \pa{\sum_{k=1}^n X_k} \end{array}

Pour construire un algorithme rapide pour apprendre un arbre de décision avec cette fonction de coût, il faut pouvoir calculer a(1, k) en fonction de a(1, k-1), b(1, k-1), X_k, Y_k ou d’autres quantités intermédiaires qui ne font pas intervenir les valeurs X_{i<k} < Y_{i<k}. D’après ce qui précède, cela paraît tout-à-fait possible. Mais dans le cas multidimensionnel, il faut déterminer le vecteur A qui minimise \sum_{k=1}^n \norme{Y - XA}^2 ce qui donne A = (X'X)^{-1} X' Y. Si on note M_{1..k} la matrice M tronquée pour ne garder que ses k premières lignes, il faudrait pouvoir calculer rapidement :

A_{k-1} - A_k = (X_{1..k-1}'X_{1..k-1})^{-1} X'_{1..k-1} Y_{1..k-1} - (X_{1..k}'X_{1..k})^{-1} X'_{1..k} Y_{1..k}

La documentation de sklearn.tree.DecisionTreeRegressor ne mentionne que deux critères pour apprendre un arbre de décision de régression, MSE pour sklearn.metrics.mean_squared_error et MAE pour sklearn.metrics.mean_absolute_error. Les autres critères n’ont probablement pas été envisagés. L’article [Acharya2016] étudie la possibilité de ne pas calculer la matrice A_k pour tous les k. Le paragraphe Streaming Linear Regression utilise le fait que la matrice A est la solution d’un problème d’optimisation quadratique et propose un algorithme de mise à jour de la matrice A (cas unidimensionnel). Cet exposé va un peu plus loin pour proposer une version qui ne calcule pas de matrices inverses.

Implémentation naïve d’une régression linéaire par morceaux#

On part du cas général qui écrit la solution d’une régression linéaire comme étant la matrice A = (X'X)^{-1} X' Y et on adapte l’implémentation de scikit-learn pour optimiser l’erreur quadratique obtenue. Ce n’est pas simple mais pas impossible. Il faut entrer dans du code cython et, pour éviter de réécrire une fonction qui multiplie et inverse une matrice, on peut utiliser la librairie LAPACK. Je ne vais pas plus loin ici car cela serait un peu hors sujet mais ce n’était pas une partie de plaisir. Cela donne : piecewise_tree_regression_criterion_linear.pyx C’est illustré toujours par le notebook DecisionTreeRegressor optimized for Linear Regression.

Aparté sur la continuité de la régression linéaire par morceaux#

Approcher la fonction y=f(x) + \epsilon quand x et y sont réels est un problème facile, trop facile… A voir le dessin, précédent, il est naturel de vouloir recoller les morceaux lorsqu’on passe d’un segment à l’autre. Il s’agit d’une optimisation sous contrainte. Il est possible également d’ajouter une contrainte de régularisation qui tient compte de cela. On exprime cela comme suit avec une régression linéaire à deux morceaux.

E = \sum_{X_i \leqslant t} (a_1 X_i + b_1 - y)^2 + \sum_{X_i \geqslant t} (a_2 X_i + b_2 - y)^2 + \lambda (a_1 t + b_1 - a_2 t - b)^2

Le cas multidimensionnel est loin d’être aussi simple. Avec une dimension, chaque zone a deux voisines. En deux dimensions, chaque zone peut en avoir plus de deux. La figure suivante montre une division de l’espace dans laquelle la zone centrale a cinq voisins.

../_images/voisin.png

Peut-on facilement approcher une fonction z = f(x,y) + \epsilon par un plan en trois dimensions ? A moins que tous les sommets soient déjà dans le même plan, c’est impossible. La zone en question n’est peut-être même pas convexe. Une régression linéaire par morceaux et continue en plusieurs dimensions n’est pas un problème facile. Cela n’empêche pas pour autant d’influencer la détermination de chaque morceaux avec une contrainte du type de celle évoquée plus haut mais pour écrire la contrainte lorsque les zones sont construites à partir des feuilles d’un arbre de décision, il faut déterminer quelles sont les feuilles voisines. Et ça c’est un problème intéressant !

Régression linéaire et corrélation#

On reprend le calcul multidimensionnel mais on s’intéresse au cas où la matrice X'X est diagonale qui correspond au cas où les variables X_1, ..., X_C ne sont pas corrélées. Si X'X = diag(\lambda_1, ..., \lambda_C) = diag(\sum_{k=1}^n X^2_{k1}, ..., \sum_{k=1}^n X^2_{kC}), la matrice A s’exprime plus simplement A = D^{-1} X' Y. On en déduit que :

a_c = \frac{\sum_{k=1}^n X_{kc} Y_k}{\sum_{k=1}^n X^2_{kc}} = \frac{\sum_{k=1}^n X_{kc} Y_k}{\lambda_c}

Cette expression donne un indice sur la résolution d’une régression linéaire pour laquelle les variables sont corrélées. Il suffit d’appliquer d’abord une ACP (Analyse en Composantes Principales) et de calculer les coefficients a_c associés à des valeurs propres non nulles. On écrit alors X'X = P'DP où la matrice P vérifie P'P = I.

Idée de l’algorithme#

On s’intéresser d’abord à la recherche d’un meilleur point de coupure. Pour ce faire, les éléments (X_i, y_i) sont triés le plus souvent selon l’ordre défini par une dimension. On note E l’erreur de prédiction sur cette échantillon E = \min_\beta \sum_k (X_k \beta - y_k)^2. On définit ensuite E(i, j) = \min_\beta \sum_{k=i}^j (X_k \beta - y_k)^2. D’après cette notation, E = E(1,n). La construction de l’arbre de décision passe par la détermination de k^* qui vérifie :

\begin{array}{rcl} E(1,k^*) + E(k^*+1, n) &=& \min_k E(1,k) + E(k+1, n) \\ &=& \min_k \pa{ \min_{\beta_1} \sum_{l=1}^k (X_l \beta_1 - y_l)^2 + \min_{\beta_2} \sum_{l=k+1}^n (X_l \beta_2 - y_l)^2} \end{array}

Autrement dit, on cherche le point de coupure qui maximise la différence entre la prédiction obtenue avec deux régressions linéaires plutôt qu’une. On sait qu’il existe une matrice P qui vérifie :

PP' = 1 \text{ et } (XP)'(XP) = P'X'XP = D = Z'Z

Où D=diag(d_1, ..., d_C) est une matrice diagonale. On a posé Z = XP, donc d_a = <Z_a, Z_a>. On peut réécrire le problème de régression comme ceci :

\beta^* = \arg\min_\beta \sum_i \norm{ y_i - X_i\beta} = \arg\min_\beta \norm{Y - X\beta}

Comme X = ZP' :

\norm{Y - X\beta} = \norm{Y - X\beta} = \norm{Y - ZP'\beta} = \norm{Y - Z\gamma}

Avec \gamma = P'\beta. C’est la même régression après un changement de repère et on la résoud de la même manière :

\gamma^* = (Z'Z)^{-1}Z'Y = D^{-1}Z'Y

La notation M_i désigne la ligne i et M_{[k]} désigne la colonne. On en déduit que le coefficient de la régression \gamma_k est égal à :

\gamma_k = \frac{<Z_{[k]},Y>}{<Z_{[k]},Z_{[k]}>} = \frac{<(XP')_{[k]},Y>}{<(XP')_{[k]},(XP')_{[k]}>}

On en déduit que :

\norm{Y - X\beta} = \norm{Y - \sum_{k=1}^{C}Z_{[k]}\frac{<Z_{[k]},Y>}{<Z_{[k]},Z_{[k]}>}} = \norm{Y - \sum_{k=1}^{C}(XP')_{[k]}\frac{<(XP')_{[k]},Y>}{<(XP')_{[k]},(XP')_{[k]}>}}

Algorithme A1 : Arbre de décision optimisé pour les régressions linéaires

On dipose qu’un nuage de points (X_i, y_i) avec X_i \in \R^d et y_i \in \R. Les points sont triés selon une dimension. On note X la matrice composée des lignes X_1, ..., X_n et le vecteur colonne