Extract trees from a random forest in python

You may need to extract trees from a classifier for various reasons. In my case, I thought that the feature of xgboost ntree_limit was quite convenient when cross validating a gradient boosting method over the number of trees.

What it does is that it only uses the first ntree_limit trees to perform the prediction (instead of using all the fitted tree).

predict(data, output_margin=False, ntree_limit=0, pred_leaf=False, 
        pred_contribs=False, approx_contribs=False, 
        pred_interactions=False, validate_features=True, training=False)

And it is also available as an extra argument of .predict() if you use the scikit-learn interface :

ypred = bst.predict(dtest, ntree_limit=bst.best_ntree_limit)

Indeed, by doing so, if you want to find the optimal number of trees for your model, you do not have to fit the model for 50 trees, and then predict, then fit it for 100 trees and then predict. You may fit the model once and for all for 200 trees and then, playing with ntree_limit you can observe the performance of the model for various number of trees.

The RandomForest, as implemented in scikit-learn does not show this parameter in its .predict() method. However, this is something we can quickly fix. Indeed, the RandomForest exposes estimators_. You can modify it (beware, this is a bit hacky and may not work for other versions of scikit-learn).

rf_model = RandomForestRegressor()
rf_model.fit(x, y)

estimators = rf_model.estimators_

def predict(w, i):
    rf_model.estimators_ = estimators[0:i+1]
    return rf_model.predict(x)

And that’s it, the predict method now only looks at the first i trees ;)

An introduction to model stacking

What is stacking and why using it

Introduction

In the case of supervised learning, stacking is a process that enables to improve the performance of a predictor. It can be used for classification and regression problems. If you took part in statistics competitions, you may already be familiar with it, but the resources about this technique are quite scarce on the internet.

From blending to stacking

As with blending, where a simple average between models with similar performance often proposes a model whose performance is higher than the one of each model in the blend, stacking combine models in a way that is dependant on the training set.

Stacking can be thought as “the sequel” of blending. Imagine you have two models with a similar performance on a dataset. The simplest blend consists in averaging the two models. However, you may become curious and propose different weight for each model. Per example, a blending with a weight of 0.7 for the first model predictions and 0.3 for the second may be better than the default 50/50 weights.

Now let’s suppose you have models and are looking for the best weights for these n models. The problem you are facing now becomes to performing a linear regression, doesn’t it ?

The only difference is that you cannot use the prediction on the training set directly (where some models like Random Forests usually have a perfect accuracy). This is where stacking comes in!

When do you need stacking ?

Staging is usually the winning solution of many data science competitions.

Per example, the Homesite quote conversion

Quick overview for now about the NMA approach:

10 variations of the dataset in total (factor combinations, factors mapped to response rates, replacing correlated pairs by differences etc) lots of models (xgboost, keras, ranger, logreg, even occasional svm - although that took forever) trained on various datasets and different params; stored as lvl1 metafeatures mix lvl1 metafeatures with: xgboost, nnet, hillclimbing, glmnet and ranger, stack - 5 lvl2 metafeatures mix the lvl2 metafeatures with hillclimbing bag at each stage as much as time permitted

I will come back to the notions of lv1 / lvl2… stages

Or, in the BNP Paribas Cardif Claims management:

We also produced many different base level models without much Feature engineering, just different input format types (like load all categorical variables as counts, or as onehot encoding etc).

Our ensemble was consisted of 223 models. Faron did a lot of work in removing noise and discarding many of these in order to get to our bets score with a lvl2 ensemble of geomean weights between an ET , 2NN and 2 Xgmodels.

And there are plenty of other examples. So basically, stacking comes in when the accuracy of your classifier or regressor is the essence of your problem. It makes the interpretability of the model really low and is harder to implement and deploy than a simple machine learning pipe.

Principles of stacking

The idea, to give a correct weight to each model is to perform a cross validation on the training set and return a dataset for which each element correspond to the unseen fold prediction. On this new dataset, you can fit the new model.

Per example, in the case of a regression problem, if you have rows in your data set, features and models, this step turns your training data from a matrix to a matrix.

In the case of a muli class problem, if you have rows in your data set, features, classes and models, this step turns your training data from a matrix to a matrix.

An example

The dataset

I will use the MNIST dataset, under a CSV format, which can be found here: on Kaggle and a logloss penalty. Of course, you will be able to play with other metrics / datasets using the code below!

The stacking / CV class

The class below can be used for a multi-class learning problem. Some minor adaptations may be required, per example for regression problems.

from sklearn.model_selection import KFold
import datetime
import pandas as pd
import numpy as np
from time import time

class bcolors:
    HEADER = '\033[95m'
    OKBLUE = '\033[94m'
    OKGREEN = '\033[92m'
    WARNING = '\033[93m'
    FAIL = '\033[91m'
    ENDC = '\033[0m'
    BOLD = '\033[1m'
    UNDERLINE = '\033[4m'


class ModelStager:

    def __init__(self, penalty, n_folds,
                 verbose=1, shuffle=True, random_state=1):
        self._penalty = penalty
        self._n_folds = n_folds
        self._verbose = verbose
        self._random_state = random_state
        self._shuffle = shuffle

    def _print(self, input_str):
        time = datetime.datetime.now().strftime("%Y-%m-%d %H:%M")
        print(bcolors.HEADER + "[ModelStager | " + time + "] " + bcolors.ENDC + str(input_str)) 

    def fit(self, X, y, model):
        kfold = KFold(n_splits=self._n_folds, shuffle=self._shuffle,
                  random_state=self._random_state)

        cv_scores = []
        oof_predictions = pd.DataFrame(index=X.index, columns=range(y.nunique()))

        fold_idx = 0

        for tr_idx, val_idx in kfold.split(X):

            X_tr = X.iloc[tr_idx]
            X_val = X.iloc[val_idx]

            y_tr = y.iloc[tr_idx]
            y_val = y.iloc[val_idx]

            if self._verbose:
                self._print("Data_tr shape : " + str(X_tr.shape))

            fold_idx = fold_idx + 1
            t = time()

            model.fit(X_tr, y_tr)

            validation_prediction = model.predict_proba(X_val)

            oof_predictions.iloc[val_idx] = validation_prediction

            cv_score_model = self._penalty(y_val, validation_prediction)
            cv_scores.append(cv_score_model)

            if self._verbose:
                self._print("Fold %.0f : TEST %.5f | TIME %.2fm (1-fold)" %
                            (fold_idx, cv_score_model, (time() - t) / 60))

        self._print("TEST AVERAGE : %.5f" % (np.mean(cv_scores)))

        return oof_predictions

As you can see, the ModelStager also performs cross validation. All the magic happens in oof_predictions, which is in charge of keeping track of the out-of-fold prediction and returning it. As mentioned earlier, it shares the index with X, and the columns correspond to the number of classes.

All the bcolors and custom printing function are just things I am used to work with, no need to bother about it.

Example

Random Forest and Extra trees

If you append this at the bottom of the previous class, you may re run the operations.

Two models are proposed, and their ensemble below (using a logistic regression).

if __name__ == "__main__":
    
    from sklearn.metrics import log_loss
    from sklearn.linear_model import LogisticRegression
    from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
    from xgboost import XGBClassifier

    train_data = pd.read_csv("./mnist_train.csv", nrows=5000)
    X = train_data.drop(["label"], axis=1)
    y = train_data["label"]
    
    stager = ModelStager(log_loss, 5)

    print("RF model")
    model_rf = RandomForestClassifier()
    stage1_rf = stager.fit(X, y, model_rf)

    print("ET model")
    model_et = ExtraTreesClassifier()
    stage1_et = stager.fit(X, y, model_et)

    print("Stage 1 : (RF, ET) -> logistic model")
    stage1_rf_et = pd.concat([stage1_rf, stage1_et], axis=1)
    stager.fit(stage1_rf_et, y, LogisticRegression())

Results

RF model
[ModelStager | 2021-01-07 13:54] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:54] Fold 1 : TEST 0.48133 | TIME 0.05m (1-fold)
[ModelStager | 2021-01-07 13:54] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 2 : TEST 0.44262 | TIME 0.05m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 3 : TEST 0.46714 | TIME 0.05m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 4 : TEST 0.45846 | TIME 0.05m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 5 : TEST 0.45377 | TIME 0.05m (1-fold)
[ModelStager | 2021-01-07 13:55] TEST AVERAGE : 0.46066
ET model
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 1 : TEST 0.44834 | TIME 0.04m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 2 : TEST 0.44679 | TIME 0.04m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 3 : TEST 0.43367 | TIME 0.04m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 4 : TEST 0.43551 | TIME 0.04m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:55] Fold 5 : TEST 0.42378 | TIME 0.04m (1-fold)
[ModelStager | 2021-01-07 13:55] TEST AVERAGE : 0.43762
Stage 1 : (RF, ET) -> logistic model
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 20)
[ModelStager | 2021-01-07 13:55] Fold 1 : TEST 0.20850 | TIME 0.01m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 20)
[ModelStager | 2021-01-07 13:55] Fold 2 : TEST 0.16870 | TIME 0.00m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 20)
[ModelStager | 2021-01-07 13:55] Fold 3 : TEST 0.21278 | TIME 0.00m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 20)
[ModelStager | 2021-01-07 13:55] Fold 4 : TEST 0.20536 | TIME 0.00m (1-fold)
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 20)
[ModelStager | 2021-01-07 13:55] Fold 5 : TEST 0.18016 | TIME 0.00m (1-fold)
[ModelStager | 2021-01-07 13:55] TEST AVERAGE : 0.19510

This is quite a huge boost ;) to be honest, this part is a little bit of an artifact, as ensemble of decision trees are usually quite bad at predicting probabilities, the logloss is artificially high. And the logistic regression corrects this phenomenon.

Random Forest, Extra Trees and Gradient Boosting

As stated above, the main performance gain comes from using an algorithm that is better at optimizing logloss.

Gradient boosting methods (most notably, xgboost) are good at predicting probability. This is illustrated when we perform the cross validation of a gradient boosting model over the original dataset.

print("XGB model")
model_xgb = XGBClassifier(use_label_encoder=False)
stage_1_xgb = stager.fit(X, y, model_xgb)

Yields the following results

XGB model
[ModelStager | 2021-01-07 13:55] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:56] Fold 1 : TEST 0.21077 | TIME 1.05m (1-fold)
[ModelStager | 2021-01-07 13:56] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:57] Fold 2 : TEST 0.16564 | TIME 1.15m (1-fold)
[ModelStager | 2021-01-07 13:57] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:58] Fold 3 : TEST 0.25023 | TIME 1.09m (1-fold)
[ModelStager | 2021-01-07 13:58] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 13:59] Fold 4 : TEST 0.24772 | TIME 1.03m (1-fold)
[ModelStager | 2021-01-07 13:59] Data_tr shape : (4000, 784)
[ModelStager | 2021-01-07 14:00] Fold 5 : TEST 0.18703 | TIME 1.12m (1-fold)
[ModelStager | 2021-01-07 14:00] TEST AVERAGE : 0.21228

Though much better than the single RandomForestClassifier or ExtraTreesClassifier alone, it does not beat the staged model.

Now let’s add the xgb features to the stage 1:

print("Stage 1 : (FR, ET, XGB) -> logistic model")
stage1_rf_et_xgb = pd.concat([stage1_rf, stage1_et, stage_1_xgb], axis=1)
stager.fit(stage1_rf_et_xgb, y, LogisticRegression())

And once again, the performance increases.

Stage 1 : (FR, ET, XGB) -> logistic model
[ModelStager | 2021-01-07 14:00] Data_tr shape : (4000, 30)
[ModelStager | 2021-01-07 14:00] Fold 1 : TEST 0.19343 | TIME 0.01m (1-fold)
[ModelStager | 2021-01-07 14:00] Data_tr shape : (4000, 30)
[ModelStager | 2021-01-07 14:00] Fold 2 : TEST 0.15602 | TIME 0.01m (1-fold)
[ModelStager | 2021-01-07 14:00] Data_tr shape : (4000, 30)
[ModelStager | 2021-01-07 14:00] Fold 3 : TEST 0.20996 | TIME 0.01m (1-fold)
[ModelStager | 2021-01-07 14:00] Data_tr shape : (4000, 30)
[ModelStager | 2021-01-07 14:00] Fold 4 : TEST 0.20830 | TIME 0.01m (1-fold)
[ModelStager | 2021-01-07 14:00] Data_tr shape : (4000, 30)
[ModelStager | 2021-01-07 14:00] Fold 5 : TEST 0.17053 | TIME 0.01m (1-fold)
[ModelStager | 2021-01-07 14:00] TEST AVERAGE : 0.18765

And if we include the gradient boosting predictions in the stage 1 features, the logloss drops from 0.19510 to 0.18765

More stage 1 features

I only presented 3 models in the stage 1. I could have added plenty of others, such as nearest neighbors, linear models… However, I strongly recommend to play with the code below and try to add these models, I am pretty sure than much better scores can be obtained ;)

Or I could also have performed some feature engineering for some models and not for others. As you can see, the number of combination becomes really huge. The “secret” to have a good performance after stacking is to have models that are as different (surprisingly, the performance of each model is not that important) as possible.

Beyond the linear model

I mostly referred to the stage 1 as a weighting operation, but it does not have to be a linear model. You can also use other model on top of your stage one features! Per example, another gradient boosting model, or a neural network. You can even repeat the above to produce stage 2 features, and train another model on this stage 2.

The code

from sklearn.model_selection import KFold
import datetime
import pandas as pd
import numpy as np
from time import time

class bcolors:
    HEADER = '\033[95m'
    OKBLUE = '\033[94m'
    OKGREEN = '\033[92m'
    WARNING = '\033[93m'
    FAIL = '\033[91m'
    ENDC = '\033[0m'
    BOLD = '\033[1m'
    UNDERLINE = '\033[4m'


class ModelStager:

    def __init__(self, penalty, n_folds,
                 verbose=1, shuffle=True, random_state=1):
        self._penalty = penalty
        self._n_folds = n_folds
        self._verbose = verbose
        self._random_state = random_state
        self._shuffle = shuffle

    def _print(self, input_str):
        time = datetime.datetime.now().strftime("%Y-%m-%d %H:%M")
        print(bcolors.HEADER + "[ModelStager | " + time + "] " + bcolors.ENDC + str(input_str)) 

    def fit(self, X, y, model):
        kfold = KFold(n_splits=self._n_folds, shuffle=self._shuffle,
                  random_state=self._random_state)

        cv_scores = []
        oof_predictions = pd.DataFrame(index=X.index, columns=range(y.nunique()))

        fold_idx = 0

        for tr_idx, val_idx in kfold.split(X):

            X_tr = X.iloc[tr_idx]
            X_val = X.iloc[val_idx]

            y_tr = y.iloc[tr_idx]
            y_val = y.iloc[val_idx]

            if self._verbose:
                self._print("Data_tr shape : " + str(X_tr.shape))

            fold_idx = fold_idx + 1
            t = time()

            model.fit(X_tr, y_tr)

            validation_prediction = model.predict_proba(X_val)

            oof_predictions.iloc[val_idx] = validation_prediction

            cv_score_model = self._penalty(y_val, validation_prediction)
            cv_scores.append(cv_score_model)

            if self._verbose:
                self._print("Fold %.0f : TEST %.5f | TIME %.2fm (1-fold)" %
                            (fold_idx, cv_score_model, (time() - t) / 60))

        self._print("TEST AVERAGE : %.5f" % (np.mean(cv_scores)))

        return oof_predictions


if __name__ == "__main__":
    
    from sklearn.metrics import log_loss
    from sklearn.linear_model import LogisticRegression
    from sklearn.ensemble import RandomForestClassifier, ExtraTreesClassifier
    from xgboost import XGBClassifier

    train_data = pd.read_csv("./mnist_train.csv", nrows=5000)
    X = train_data.drop(["label"], axis=1)
    y = train_data["label"]
    
    stager = ModelStager(log_loss, 5)

    print("RF model")
    model_rf = RandomForestClassifier()
    stage1_rf = stager.fit(X, y, model_rf)

    print("ET model")
    model_et = ExtraTreesClassifier()
    stage1_et = stager.fit(X, y, model_et)

    print("Stage 1 : (RF, ET) -> logistic model")
    stage1_rf_et = pd.concat([stage1_rf, stage1_et], axis=1)
    stager.fit(stage1_rf_et, y, LogisticRegression())

    print("XGB model")
    model_xgb = XGBClassifier(use_label_encoder=False)
    stage_1_xgb = stager.fit(X, y, model_xgb)

    print("Stage 1 : (FR, ET, XGB) -> logistic model")
    stage1_rf_et_xgb = pd.concat([stage1_rf, stage1_et, stage_1_xgb], axis=1)
    stager.fit(stage1_rf_et_xgb, y, LogisticRegression())
P-adic numbers visualization

What are p-adic numbers ?

P-adic and rationnal numbers

P-adic numbers are an original way to look at the (limit of sequence of) elements in .

More precisely, just like represents the limits of Cauchy sequences in endowed with the distance : , represents the limits of Cauchy sequences in with another distance : , where is detailed below.

P-adic valuations

For a prime number, define as the exponent of in the decomposition of in a product of prime factors. Also define

Then is a distance on integers.

In with the distance , note that the sequence converges towards .

Why they matter

Various results can be proved using p-adic numbers. I discovered them in “Introduction to number theory”, where they are used to determine whether an ellipse has rationnal points. They also enable to give a meaning to

Visualization

The idea

A p-adic number can be written where the sum might be infinite. Though it seems weird because the terms are growing, note that the sequence actually tends to really quickly in

A traditionnal way to picture p-adic numbers is with co-centric circles, like below:

Representation of p-adic integers

All the credit goes to: Heiko Knopse for this illustration, more are available on his site

My idea is to take this idea to the limit. Formally, for , the complex number is associated to .

is a parameter between and used to ensure convergence.

Results

Representing some integers

Some integers in zp

Convergence

Convergence of a sequence in zp

Addition

An interesting property is that . It is illustrated below. As you can see, addition in the p-addic representation shifts numbers to the right.

Convergence of a sequence in zp

Code

from cmath import *


class PAddicRepresenter:

    def __init__(self, p, l, output_length=30):
        self._p = p
        self._l = l
        self._output_length = output_length

    def to_plane(self, n):
        l = self._l
        p = self._p
        decomposed_int = self._completed_int_to_base(n)
        complex_coordinates = sum(
            [l ** n * exp(1j * c * 2 * pi / p) for n, c in enumerate(decomposed_int)])
        return complex_coordinates.real, complex_coordinates.imag

    def transform_sample(self, ns):
        xs, ys = [], []

        for n in ns:
            x, y = self.to_plane(n)
            xs.append(x)
            ys.append(y)

        return xs, ys

    def _int_to_base(self, n):
        p = self._p
        i = 0
        decomposition = []
        while n > 0:
            residual = n % p
            n = (n - residual) / p
            decomposition.append(residual)
        return decomposition

    def _completed_int_to_base(self, n):
        decomposed_int = self._int_to_base(n)
        return decomposed_int + [0] * (self._output_length - len(decomposed_int))

The first visualization being obtaining using the following:

import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (8,8)
from PAddicRepresenter import PAddicRepresenter


n_points = 3**10
p = 3
small_sample_size = 55
l = 0.45

par = PAddicRepresenter(p, l)

xs, ys = par.transform_sample(range(n_points))

fig, ax = plt.subplots()

ax.hist2d(xs, ys, bins = 500, cmap = 'Greys')

ax.scatter(xs[0:small_sample_size], ys[0:small_sample_size], c='black')
for i in range(small_sample_size):
    ax.annotate(str(i), (xs[i] - 0.03 , ys[i] + 0.05))
 
plt.axis('off')
plt.show()

Learning more

For those interested in number theory, I strongly recommend :

Number Theory 1: Fermat’s Dream by Kazuya Kato, Nobushige Kurokawa and Takeshi Saito

Number Theory 2: Introduction to Class Field Theory by the same authors which requires more knowledge in algebra and group theory.

OCaml Bigarray vs array

A simple benchmark

I was recently wondering if I could speed up the access to the element of a matrix in OCaml, using Bigarrays. As the benchmark below shows, it turns out that no, Bigarray was actually slower than matrix operations.

However, a couple of interesting observations can also be made: inlining manually provided a huge increase of the performance for nested arrays, and Owl is actually the fastest solution available (but not by a lot).

open Owl

let n = 10000


let array_at m i j =
  Array.unsafe_get (Array.unsafe_get m i ) j


let sum_array_no_inline n m =
    let a = ref 0. in
    for i = 0 to (n-1) do
      for j = 0 to (n-1) do
        a := !a +. array_at m i j 
      done;
    done;
    !a

let sum_array n m = 
  let a = ref 0. in
  for i = 0 to (n-1) do
    for j = 0 to (n-1) do
      a := !a +. (Array.unsafe_get (Array.unsafe_get m i ) j )
    done;
  done;
  !a


let sum_big_array n m =
  let a = ref 0. in
  for i = 0 to (n-1) do
    for j = 0 to (n-1) do
      a := !a +. (Bigarray.Array2.unsafe_get m i j);
    done;
  done;
  !a


let sum_owl_array n m =
  let a = ref 0. in
  for i = 0 to (n-1) do
    for j = 0 to (n-1) do
      a := !a +. Mat.get m i j;
    done;
  done;
  !a


let sum_owl_array_lib n m =
  Mat.sum' m


let time f x =
  let start = Unix.gettimeofday () in 
  let res = f x in 
  let stop = Unix.gettimeofday () in 
  Printf.printf "Execution time: %fs\n%!" (stop -. start);
  res


let () =
  let arr = Array.init n (fun i -> Array.init n (fun j -> (Random.float 1.0))) in
  
  print_string "[ArrayNoInline] ";
  let a1 = time (sum_array_no_inline n) arr in

  print_string "[Array] ";
  let a2 = time (sum_array n) arr in

  let big_arr = Bigarray.Array2.of_array Bigarray.float32 Bigarray.c_layout arr in
  print_string "[BigArray] ";
  let b = time (sum_big_array n) big_arr in

  let owl_arr = Mat.of_arrays arr in
  print_string "[OwlArray] ";
  let c = time (sum_owl_array n) owl_arr in
  print_string "[OwlArrayLib] ";
  let d = time (sum_owl_array_lib n) owl_arr in

  print_string "\n";
  print_float a1;
 
  print_string "\n";
  print_float a2;
 
  print_string "\n";
  print_float b;

  print_string "\n";
  print_float c;

  print_string "\n";
  print_float d;
  ()

And the results

[ArrayNoInline] Execution time: 0.432230s
[Array] Execution time: 0.105445s
[BigArray] Execution time: 3.037937s
[OwlArray] Execution time: 2.177349s
[OwlArrayLib] Execution time: 0.080217s
Python plot 3d scatter and density

It is often easy to compare, in dimension one, an histogram and the underlying density. This is quite useful when one want to visually evaluate the goodness of fit between the data and the model. Unfortunately, as soon as the dimesion goes higher, this visualization is harder to obtain. Here, I will present a short snippet rendering the following plot:

3d scatter plot and density

The heatmap is flat, on top of it, a wireframe is plotted and the sampled points are constrained to have the same height as the wireframe, so that their density is more visual.

Feel free to use the snippet below :)

from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.stats import multivariate_normal


# Sample parameters
mu = np.array([0, 0])
sigma = np.array([[0.7, 0.2], [0.2, 0.3]])
rv = multivariate_normal(mu, sigma)
sample = rv.rvs(500)

# Bounds parameters
x_abs = 2.5
y_abs = 2.5
x_grid, y_grid = np.mgrid[-x_abs:x_abs:.02, -y_abs:y_abs:.02]

pos = np.empty(x_grid.shape + (2,))
pos[:, :, 0] = x_grid
pos[:, :, 1] = y_grid

levels = np.linspace(0, 1, 40)

fig = plt.figure()
ax = fig.gca(projection='3d')

# Removes the grey panes in 3d plots
ax.xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))

# The heatmap
ax.contourf(x_grid, y_grid, 0.1 * rv.pdf(pos),
            zdir='z', levels=0.1 * levels, alpha=0.9)

# The wireframe
ax.plot_wireframe(x_grid, y_grid, rv.pdf(
    pos), rstride=10, cstride=10, color='k')

# The scatter. Note that the altitude is defined based on the pdf of the
# random variable
ax.scatter(sample[:, 0], sample[:, 1], 1.05 * rv.pdf(sample), c='k')

ax.legend()
ax.set_title("Gaussian sample and pdf")
ax.set_xlim3d(-x_abs, x_abs)
ax.set_ylim3d(-y_abs, y_abs)
ax.set_zlim3d(0, 1)

plt.show()

Learning more

Data Visualization with Python for Beginners and Matplotlib 3.0 Cookbook are complete references for using Matplotlib and Seaborn

Best casual readings in mathematics

If you are reading this blog, you probably have a science degree ;) after graduating the opportunities to use mathematics in the day job are often limited to a small subset of what was learnt. The books below, varying in difficulty, are an occasion to practice or to have “mathematical recreations”.

I try to rate the level of technicality of the book *** are the math books you would expect to read to prepare for a master’s degree exam, ** still require a pencil and a paper to explore the details and * or less contain very few technicalities.

Miscellanous

These are among my favorites: they require very little knowledge about any specific topic in mathematics and are pure moments of cleverness, smart arguments and nice illustrations.

(**) Proofs from the book also available in PDF. When Erdos proved something and the proof looked clumsy to him, he was saying “This is not the proof from the book” or “Let’s look for the proof of the book”. And indeed, when proving statements, some proofs seem more natural than others: usually, the shortest and most convincing ones. This book is an attempt - a successful one - to gather them. The infinity of prime numbers, D’Alembert Gauss theorem, the Law of quadratic reciprocity are just a few examples of all the results presented in this book.

Some proofs may be profitable to the reading of Number Theory 1: Fermat’s Dream.

(**) The Art of Mathematics: Coffee Time in Memphis more than a hundred exercises, with hints and detailed solutions. The difficulty varies greatly from an exercise to another and the solutions come from many different fields. Ideal for long trips ;)

Algebra an number theory

If you like number theory, the following books are a must have. The first volume is easily accessible, however, the following ones will require a working knowledge of Galois theory. I love the way authors present the intuitions behind the proofs and the main steps to go through before actually “jumping” into the proof.

(**) Number Theory 1: Fermat’s Dream

(***) Number Theory 2: Introduction to Class Field Theory If the previous book presented some facts that looked “magic”, this one focuses on explaining why they happen. It is much harder than the previous read. I would strongly advise this read to those who loved the previous one.

(***) Number Theory 3: Iwasawa Theory and Modular Forms The same advice applies ;)

If you do not know about Galois theory, these two references may help. They do not qualify as “casual readings” but they will help understand the “Number theory saga”.

(**) Galois Theory for Beginners: A Historical Perspective

(***) Galois Theory by Emil Artin Though this course is quite old, the book gives a clear presentation of the topic

Applications of mathematics in real life

Music

(**) Music: A Mathematical Offering by Dave Benson covers many topics about the interplay between mathematics and music.

Finance

(*) A Practitioner’s Guide to Asset Allocation by William Kinlaw, Mark P. Kritzman, David Turkington. The mathematical details are very light in this book, the focus is put on the models, the history and controversies around models and some actual data about typical correlations, returns of asset classes (this is scarcer than one would think for a book of finance!). If you want to invest by yourself and know about how to diversify, this is a very good starting point.

(**) Portfolio Optimization and Performance Analysis by Jean-Luc Prigent. If you liked the previous book and want to dig (much deeper) in portfolio optimization, this book is a detailed analysis of the existing models.

History

(*) Euler: The Master of Us All great book about the works of Euler, as the title indicates. The chapter are organized according to the branches of mathematics Euler contributed to (all of them, at his time), and the proofs are the proofs he presented at the time.

Théorème vivant (in French) this one is hard to describe. Do not expect to understand precisely the contents of Cedric Villani’s work by reading this book. Likewise, the equations come with few explanations. It is more like a diary of a researcher.

LightGBM on the GPU

LightGBM is currently one of the best implementations of gradient boosting. I will not go in the details of this library in this post, but it is the fastest and most accurate way to train gradient boosting algorithms. And it has a GPU support. Installing something for the GPU is often tedious… Let’s try it!

Setting up LightGBM with your GPU

I will assume a nVidia GPU. I personnally have a GeForce GTX 745, with the Driver Version: 410.48. If you do not have a GPU already, be careful in the model you chose. When buying a GPU, you have to make sure the “compute capability” is high enough with respect to the software you plan to use. Per example, rapids.ai needs at leats a NVIDIA Pascal™ GPU or better with compute capability 6.0+. I am not going to discuss about rapids.ai in this post, but if you plan to install LightGBM on your GPU, you will soon enough want to play with rapids.ai as well.

Your simplest choice is probably : GTX 1660 Ti which was released in february 2019 and has a compute capability of 7.5

Among older GPUs which have a compute capability of 6+, the prices change quite often but you could make a good deal below.

Nvidia TITAN Xp, GeForce GTX 1080 Ti, GTX 1080, GTX 1070 Ti, GTX 1070, GTX 1050,

But keep in mind that with these cards, the support may be abandoned soon enough…

Test if LightGBM supports GPU

If you can run the following python script:

from lightgbm import LGBMClassifier
from sklearn.datasets import make_moons

model = LGBMClassifier(boosting_type='gbdt', num_leaves=31, max_depth=- 1, learning_rate=0.1, n_estimators=300, device = "gpu")

train, label = make_moons(n_samples=300000, shuffle=True, noise=0.3, random_state=None)

model.fit(train, label)

Without this message:

[LightGBM] [Fatal] GPU Tree Learner was not enabled in this build.
Please recompile with CMake option -DUSE_GPU=1
Traceback (most recent call last):
  File "seq.py", line 11, in <module>
    model.fit(train, label)
  File "/home/kerneltrip/anaconda3/lib/python3.7/site-packages/lightgbm/sklearn.py", line 800, in fit
    callbacks=callbacks)
  File "/home/kerneltrip/anaconda3/lib/python3.7/site-packages/lightgbm/sklearn.py", line 595, in fit
    callbacks=callbacks)
  File "/home/kerneltrip/anaconda3/lib/python3.7/site-packages/lightgbm/engine.py", line 228, in train
    booster = Booster(params=params, train_set=train_set)
  File "/home/kerneltrip/anaconda3/lib/python3.7/site-packages/lightgbm/basic.py", line 1666, in __init__
    ctypes.byref(self.handle)))
  File "/home/kerneltrip/anaconda3/lib/python3.7/site-packages/lightgbm/basic.py", line 47, in _safe_call
    raise LightGBMError(decode_string(_LIB.LGBM_GetLastError()))
lightgbm.basic.LightGBMError: GPU Tree Learner was not enabled in this build.
Please recompile with CMake option -DUSE_GPU=1

Then, you do not need this tutorial ;)

Setup guide

Though there is some information here, following the instructions did not do the job for me (hence this detailed guide).

GPU drivers

First, you need to have you drivers set up.

sudo add-apt-repository ppa:graphics-drivers/ppa 
sudo apt update 

You may find your device and the drivers using:

:~$ ubuntu-drivers devices
== /sys/devices/pci0000:00/0000:00:01.0/0000:01:00.0 ==
vendor   : NVIDIA Corporation
[...]
model    : GM107 [GeForce GTX 745]
driver   : nvidia-340 - distro non-free
driver   : nvidia-384 - distro non-free
driver   : nvidia-410 - third-party non-free recommended
driver   : nvidia-396 - third-party non-free
[...]

And then, you can run the following, where 410 replaces the recommended version of the driver:

sudo apt-get update
sudo apt-get install --no-install-recommends nvidia-410
sudo apt-get install --no-install-recommends nvidia-opencl-icd-410 nvidia-opencl-dev opencl-headers

LGBM dependencies

The officials instructions are the following, first the prerequisites:

sudo apt-get install --no-install-recommends git cmake build-essential libboost-dev libboost-system-dev libboost-filesystem-dev

(For some reason, I was still missing Boost elements as we will see later)

Building LGBM for the GPU

Time to download LightGBM

git clone --recursive https://github.com/microsoft/LightGBM
cd LightGBM
mkdir build ; cd build

Let’s try:

cmake -DUSE_GPU=1 ..

Unfortunately, this does not work.

CMake Error at /usr/local/share/cmake-3.17/Modules/FindPackageHandleStandardArgs.cmake:164 (message):
  Could NOT find OpenCL (missing: OpenCL_LIBRARY OpenCL_INCLUDE_DIR)

The official instructions helped, but :

cmake -DUSE_GPU=1 -DOpenCL_LIBRARY=/usr/local/cuda/lib64/libOpenCL.so -DOpenCL_INCLUDE_DIR=/usr/local/cuda/include/ ..

Still does not work. Now the errors about OpenCL are replaced with others related to Boost :

CMake Error at /usr/local/share/cmake-3.17/Modules/FindPackageHandleStandardArgs.cmake:164 (message):
  Could NOT find Boost (missing: Boost_INCLUDE_DIR filesystem system)
  (Required is at least version "1.56.0")
Call Stack (most recent call first):
  /usr/local/share/cmake-3.17/Modules/FindPackageHandleStandardArgs.cmake:445 (_FPHSA_FAILURE_MESSAGE)
  /usr/local/share/cmake-3.17/Modules/FindBoost.cmake:2145 (find_package_handle_standard_args)
  CMakeLists.txt:121 (find_package)

This might be overkill but:

sudo apt-get install libboost-all-dev

Did the job.

The following additional packages will be installed:
  icu-devtools libboost-atomic-dev libboost-atomic1.58-dev libboost-atomic1.58.0
  libboost-chrono-dev libboost-chrono1.58-dev libboost-chrono1.58.0 libboost-context-dev
  libboost-context1.58-dev libboost-context1.58.0 libboost-coroutine-dev libboost-coroutine1.58-dev
  libboost-coroutine1.58.0 libboost-date-time-dev libboost-date-time1.58-dev libboost-dev
  [...]

Finally:

cmake -DUSE_GPU=1 -DOpenCL_LIBRARY=/usr/local/cuda/lib64/libOpenCL.so -DOpenCL_INCLUDE_DIR=/usr/local/cuda/include/ ..

Outputs:

/cuda/lib64/libOpenCL.so -DOpenCL_INCLUDE_DIR=/usr/local/cuda/include/ ..
-- OpenCL include directory: /usr/local/cuda/include
-- Found Boost: /usr/include (found suitable version "1.58.0", minimum required is "1.56.0") found components: filesystem system 
-- Performing Test MM_PREFETCH
-- Performing Test MM_PREFETCH - Success
-- Using _mm_prefetch
-- Performing Test MM_MALLOC
-- Performing Test MM_MALLOC - Success
-- Using _mm_malloc
-- Configuring done
-- Generating done
-- Build files have been written to: /home/kerneltrip/Codes/LightGBM/build

And I can run:

make -j$(nproc)

Should look like this :

Scanning dependencies of target lightgbm
Scanning dependencies of target _lightgbm
[  3%] Building CXX object CMakeFiles/_lightgbm.dir/src/boosting/gbdt_model_text.cpp.o
[  3%] Building CXX object CMakeFiles/lightgbm.dir/src/main.cpp.o
[  4%] Building CXX object CMakeFiles/lightgbm.dir/src/application/application.cpp.o
[...]
[ 98%] Built target _lightgbm
[100%] Linking CXX executable ../lightgbm
[100%] Built target lightgbm

Everything was successfully built! Time to set up the python. The repo should look like this on your machine :

~/Codes/LightGBM$ tree -d -L 1
.
├── build
├── compute
├── docker
├── docs
├── examples
├── helpers
├── include
├── pmml
├── python-package
├── R-package
├── src
├── swig
├── tests
└── windows

The official instructions recommend the following operations, but I would not recommend them.

sudo apt-get -y install python-pip
sudo -H pip install setuptools numpy scipy scikit-learn -U
cd python-package/
sudo python setup.py install --precompile
cd ..

install python-pip has the habit of conflicting with the pip that you may have.

Instead, install the missing packages step by step.

In my case, I use conda and I was only missing setuptools.

conda install setuptools
python setup.py install --precompile

Did the job! Now…

from lightgbm import LGBMClassifier
from sklearn.datasets import make_moons


model = LGBMClassifier(boosting_type='gbdt', num_leaves=31, max_depth=- 1, learning_rate=0.1, n_estimators=300, device = "gpu")

train, label = make_moons(n_samples=300000, shuffle=True, noise=0.3, random_state=None)

model.fit(train, label)

Run without any issues ! You can observe the GPU usage with glances[gpu] (this will be fast though)

An error message that did not appear on CPU:

Unfortunately, you may find this error message with the GPU (on some datasets), which you did not have on the CPU :(

    raise LightGBMError(decode_string(_LIB.LGBM_GetLastError()))
lightgbm.basic.LightGBMError: Check failed: (best_split_info.left_count) > (0) at /home/kerneltrip/Codes/LightGBM/src/treelearner/serial_tree_learner.cpp, line 613 .

The most recent information I could get is : https://github.com/microsoft/LightGBM/issues/2742 Apparently, the issue happened on the CPU, then it was fixed, but not on the GPU version. This issue has only been raised some hours ago, let’s hope it will be fixed soon enough.

Using postMessage between iframes

There were plenty of resources regarding the use of postMessage here and there, however, none focused on reproducing the bugs locally. Just “copy this line in the iframe, this line in the page hosting the iframe and everything will work.”

I had the same issue and wanted to fix it and make sure it was running without long build, deploys etc

Here we will start from scratch, using python to serve pages locally. First we will serve two pages on the same port, then ensure postMessage between these pages works, break it with serving it on different ports and finally fix the iframe communication.

Note that I do not focus on the origin of the event checks below, but developer.mozilla.org does.

Sending data to a parent frame with postMessage

Step 1 : creating and serving two pages

Let’s start with serving a page, containing an iframe locally, at the same address. For this you will need a web browser, and a working python distribution. In what follows, everything will be served on port 4000.

index.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>I am the main page</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>

	<body>
		<h1>I am the main page</h1>  
		<iframe id="iframe" src="http://localhost:4000/iframe.html" height="600" width="800">
		</iframe>
	</body>
</html>

iframe.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>Iframe</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>

	<body>
		<h1>I am the iframe</h1>  
		<p>I should be in the site</p>
	</body>
</html>

Side note, I remembered that serving a page locally involved a “python http.server”. The following trick should be known by everybody using terminal, if you want to avoid googling again and again ;)

user@user:~$ history | grep http.se
 1283  python3 -m http.server 4000
 [...]

The command I was looking for was:

python3 -m http.server 4000

Opening localhost in your browser should look like…

iframe sample

Step 2 : postMessage to parent frame

Let’s add some javascript in the pages so that the iframe can send data to the parent frame. Basically, what happens below is that the iframe manages to run a function located in the parent window.

You can check it by clicking the button once the pages are served.

index.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>I am the main page</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>

	<body>
		<h1>I am the main page</h1>  
		<iframe id="iframe" src="http://localhost:4000/iframe.html" height="600" width="800">
		</iframe>
	</body>
</html>

iframe.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>Iframe</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>


	<body>
		<h1>I am the iframe</h1>  
		<p>I should be in the site</p>
		<button id="my-btn">Send to parent</button>
	</body>
</html>

<script>
	document.getElementById('my-btn').addEventListener('click', handleButtonClick, false);
	function handleButtonClick(e) {
		window.parent.postMessage('message');
		console.log("Button clicked in the frame");
	}	
</script>

Step 2: emulating cross domain, locally

Triggering the error

Now let’s try to serve the frame and the site on “different domains”. For this, we need to move the pages in different folders, so that it looks like:

.
├── iframe
│   └── index.html
└── index
    └── index.html

We need to run two instances of http.serve, the main page will be served on port 4000 and the iframe on port 4200.

user@user:~/Codes/iframe-communication/index$ python3 -m http.server 4000

and

user@user:~/Codes/iframe-communication/iframe$ python3 -m http.server 4200

While the contents of the pages will be (note that the port to the iframe has to be changed):

index/index.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>I am the main page</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>


	<body>
		<h1>I am the main page</h1>  
		<iframe id="iframe" src="http://localhost:4200" height="600" width="800">
		</iframe>
	</body>
</html>

<script>
	window.addEventListener('message', function() {
		alert("Message received");
		console.log("message received");
	}, false);
</script>

iframe/index.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>Iframe</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>


	<body>
		<h1>I am the iframe</h1>  
		<p>I should be in the site</p>
		<button id="my-btn">Send to parent</button>
	</body>
</html>

<script>
	document.getElementById('my-btn').addEventListener('click', handleButtonClick, false);
	function handleButtonClick(e) {
		window.parent.postMessage('message');
		console.log("Button clicked in the frame");
	}	
</script>

Clicking on the button, you should see the following error in your console:

(index):24 Failed to execute 'postMessage' on 'DOMWindow': 
The target origin provided ('null') does not match the recipient 
window's origin ('http://localhost:4000').

Fixing the error

The message is clear enough, there is a missing argument : the iframe should now the address of the parent.

Writing the following:

index/index.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>I am the main page</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>


	<body>
		<h1>I am the main page</h1>  
		<iframe id="iframe" src="http://localhost:4200" height="600" width="800">
		</iframe>
	</body>
</html>

<script>
	window.addEventListener('message', function() {
		alert("Message received");
		console.log("message received");
	}, false);
</script>

iframe/index.html

<!doctype html>

<html lang="en">
	<head>
		<meta charset="utf-8">

		<title>Iframe</title>
		<meta name="description" content="Iframe communication tutorial">
		<meta name="author" content="TheKernelTrip">

	</head>


	<body>
		<h1>I am the iframe</h1>  
		<p>I should be in the site</p>
		<button id="my-btn">Send to parent</button>
	</body>
</html>

<script>
	document.getElementById('my-btn').addEventListener('click', handleButtonClick, false);
	function handleButtonClick(e) {
		window.parent.postMessage('message', "http://localhost:4000");
		console.log("Button clicked in the frame");
	}	
</script>
Resources in topological data analysis

Topological Data Analysis is a promising field in data analysis. Indeed, the hypothesis made on the data are much weaker than the hypothesis behind “usual” machine learning algorithms. Quarantine may be a great time for learning it ;)

This list is not exhaustive, I may have missed important resources. If so, please let me know in the comments!

Articles

Mathematical topic presentations

The following papers are great presentations of the topic:

TDA applications

Books

Videos

A private company, Ayasdi has published many videos on the topic. Most of them are related to real world applications of TDA.

And their youtube channel has many others.

Software

Python

R

tSNE vs PCA

Introduction

Both PCA and tSNE are well known methods to perform dimension reduction. The question of their difference is often asked and here, I will present various points of view: theoretical, computational and emprical to study their differences. As usual, one method is not “better” in every sense than the other, and we will see that their successes vastly depend on the dataset and that a method may preserve some features of the data, while the other do not.

A little bit of wording and notations

PCA stands for Principal Component Analysis

whereas tSNE stands for Stochastic Neighbor Embedding, the t itself referring to the Student-t kernel.

As “usual” we will call the number of observations and the number of features.

Theoretical aspects

tSNE

Quoting the authors of [1]:

Stochastic Neighbor Embedding (SNE) starts by converting the high-dimensional Euclidean distances between data points into conditional probabilities that represent similarities. The similarity of data point to datapoint is the conditional probability , that would pick as its neighbor if neighbors were picked in proportion to their probability density under a Gaussian centered at . For nearby datapoints, relatively high,i whereas for widely separated datapoints, will be almost infinitesimal.

In a nutshell, this maps your data in a lower dimensional space, where a small distance between two points means that they were close in the original space.

The method goes on defining:

$$p_{j \vert i} = \frac{\exp(-\lVert x_i - x_j \rVert^2 / 2\sigma_i^2)}{\sum_{k \neq i} \exp(-\lVert x_k - x_j \rVert^2 / 2\sigma_i^2)}$$

As the probabilities in the original space, and the same probability, but in the space of represented points.

$$q_{j \vert i} = \frac{\exp(-\lVert y_i - y_j \rVert ^2)}{\sum_{k \neq i} \exp(-\lVert y_k - y_j \rVert^2)}$$

Then, the divergence between these distribution is minimized using a gradient descent over the elements:

$$C = \sum_{i,j}p_{j \vert i} \log \frac{p_{j \vert i}}{q_{j \vert i} } $$

PCA

On the other hand, PCA does not rely on a probabilistic approach to represent the points. Indeed, the approach is more geometric (or related to linear algebra): PCA looks for axis that explains as much variance as possible of the data.

Following the wikipedia article, we can have a better understanding of how it performs dimension reduction:

The transformation maps a data vector from an original space of variables to a new space of variables which are uncorrelated over the dataset. However, not all the principal components need to be kept. Keeping only the first principal components, produced by using only the first eigenvectors, gives the truncated transformation

$$\mathbf{T}_L = \mathbf{X} \mathbf{W}_L$$

where the matrix now has rows but only columns. In other words, PCA learns a linear transformation where the columns of matrix form an orthogonal basis for the features (the components of representation ) that are decorrelated. By construction, of all the transformed data matrices with only columns, this score matrix maximises the variance in the original data that has been preserved, while minimising the total squared reconstruction error or .

What we see here is that the total squared reconstruction error is minimized. Therefore, if two points are far away in the original space, their distance will have to be important as well once they are mapped to a lower dimensional space.

In that case, this means that after a PCA, two points are far away from each other if they were far away from each other in the original data set.

Feasability

Computational complexity

If you are not familiar with this topic, I already wrote about it in computational complexity and machine learning.

PCA can be performed quite quickly : it consists in evaluating the covariance matrix of the data, and performing an eigen value of this matrix. Covariance matrix computation is performed in operations while its eigen-value decomposition is operations. So, the complexity of PCA is . Considering that is fixed, this is .

On the other hand, tSNE can benefit form the Barnes-Hut approximation [2], making it usable in high dimensions with a complexity of (more informations can be found in the article).

Therefore, it is not obvious if a method has to be preferred to another with large datasets.

Parameters

The PCA is parameter free whereas the tSNE has many parameters, some related to the problem specification (perplexity, early_exaggeration), others related to the gradient descent part of the algorithm.

Indeed, in the theoretical part, we saw that PCA has a clear meaning once the number of axis has been set.

However, we saw that appeared in the penalty function to optimize for the tSNE. Besides, as a gradient descent is required, all the usual parameters of the gradient descent will have to be specified as well (step sizes, iterations, stopping criterions…).

Looking at the implementations, the arguments of the constructors in scikit-learn are the following:

class sklearn.decomposition.PCA(n_components=None, 
		copy=True, whiten=False, svd_solver='auto', 
		tol=0.0, iterated_power='auto', 
		random_state=None)

The parameters mostly relate to the solvers, but the method is unique. On the other hand…

class sklearn.manifold.TSNE(n_components=2, perplexity=30.0, 
		early_exaggeration=12.0, learning_rate=200.0, 
		n_iter=1000, n_iter_without_progress=300, 
		min_grad_norm=1e-07, metric='euclidean', init='random', verbose=0, 
		random_state=None, method='barnes_hut', angle=0.5, n_jobs=None)

However, default parameters usually provide nice results. From the original paper:

The performance of SNE is fairly robust to changes in the perplexity, and typical values are between 5 and 50.

Therefore, both in terms of feasability or efforts demanded by the calibration procedure, there is no reason to prefer a method to the other.

Empirical analysis

Figures

The following shows the reduction from a three dimensional space to a two dimensional space using both methods. The colors of the points are preserved, so that one can “keep track” of them.

Gaussian blobs

Gaussian blobs

Fig. 1: Gaussian blobs in three dimensions

Gaussian blobs - PCA

Fig. 2: Gaussian blobs after PCA

Gaussian blobs

Fig. 3: Gaussian blobs after tSNE

Both methods were successful at this task: the blobs, that were initally well separated, remain well separated in the lower dimensional space. An interesting phenomenon, which validates what the theoretical arguments predicted is that in the case of the PCA, the (light) blue and cyan points are far away from each other, whereas they appear to be closer when the tSNE is used.

The swiss roll

Swiss roll

Fig. 4: Swiss roll in three dimensions

Swiss roll - PCA

Fig. 5: Swiss roll after PCA

Swiss roll

Fig. 6: Swiss roll after tSNE

Somehow the roll is broken by the tSNE, which is weird because one would expect the red dots to be close to the orange dots… On the other hand, a linear classifier would be more successful on the data represented with the tSNE than with the PCA.

Digit dataset

Swiss roll - PCA

Fig. 7: Digits after PCA

Swiss roll

Fig. 8: Digits after tSNE

The tSNE works amazingly well on this data set, and exhibits a neat separation of most of the digits!

Hand dataset

Hand

Fig. 9: Hand data

Hand - PCA

Fig. 10: Hand - PCA

Hand - tSNE

Fig. 11: Hand - tSNE

Here we note that the fingers “remain together” with the tSNE.

Other observations

Other observations could be inferred as well, per example, the size of a cluster does not mean much with the tSNE, while it has a meaning in the case of the PCA.

Gaussian blobs

Fig. 12: Gaussian blobs in three dimensions

Gaussian blobs - PCA

Fig. 13: Gaussian blobs after PCA

Gaussian blobs

Fig. 14: Gaussian blobs after tSNE

Code

Shall you want to explore new datasets, feel free to use the following code!

import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as p3
from sklearn.datasets import make_swiss_roll, make_blobs, load_digits
from sklearn import decomposition
from sklearn.manifold import TSNE

# Generate data (swiss roll dataset)
n_samples = 2500

def generate_swiss_roll_data(n_samples):
    noise = 0.05
    X, _ = make_swiss_roll(n_samples, noise)
    # Make it thinner
    X[:, 1] *= .5
    distance_from_y_axis = X[:, 0] ** 2 + X[:, 2] ** 2 
    X_color = plt.cm.jet(distance_from_y_axis / np.max(distance_from_y_axis + 1))
    return X, X_color, "Swiss roll"

def generate_gaussian_blobs(n_samples):
    X, y = make_blobs(n_samples=n_samples, centers=5, n_features=3,random_state=3)
    X_color = plt.cm.jet(y / np.max(y + 1))
    return X, X_color, "Gaussian blobs"

def generate_gaussian_blobs_different(n_samples):
    X, y = make_blobs(n_samples=n_samples, centers=2, n_features=3,random_state=3)
    X[y == 1, :] *= 2
    X_color = plt.cm.jet(y / np.max(y + 1))
    return X, X_color, "Gaussian blobs different sizes"

def generate_digits(n_samples):
    X, y = load_digits(n_class = 10, return_X_y = True)
    X_color = plt.cm.jet(y / np.max(y + 1))
    return X, X_color, "Digits"

def generate_hand(n_samples):
    X = np.loadtxt("./Hand.csv", delimiter=",")
    z_component =  X[:, 2] - X[:, 0]
    X_color = plt.cm.jet(z_component / np.max(z_component + 1))
    return X, X_color, "Hand"


def produce_plots(data_generator, data_generator_name):
    X, X_color, title = data_generator(n_samples)
    fig = plt.figure()
    ax = p3.Axes3D(fig)
    ax.view_init(7, -80)
    ax.scatter(X[:, 0], X[:, 1], X[:, 2],
            color = X_color,
            s=20, edgecolor='k')
    ax.set_xlabel('$x$')
    ax.set_ylabel('$y$')
    ax.set_zlabel('$z$')
    plt.title(title)
    plt.savefig(data_generator_name + '.png')

    # Fit and plot PCA
    X_pca = decomposition.PCA(n_components=2).fit_transform(X)
    fig = plt.figure()
    s = fig.add_subplot(1, 1, 1, xlabel='$x$', ylabel='$y$')
    s.scatter(X_pca[:, 0], X_pca[:, 1], c = X_color)
    plt.title(title + " - PCA")
    fig.savefig(data_generator_name  + '_pca.png')

    # Fit and plot tSNE
    X_tsne = TSNE(n_components=2).fit_transform(X)
    fig = plt.figure()
    s = fig.add_subplot(1, 1, 1, xlabel='$x$', ylabel='$y$')
    s.scatter(X_tsne[:, 0], X_tsne[:, 1], c = X_color)
    plt.title(title + " - tSNE")
    fig.savefig(data_generator_name + '_tsne.png')

produce_plots(generate_digits, "digits")
produce_plots(generate_hand, "hand")
produce_plots(generate_swiss_roll_data, "swiss_roll")
produce_plots(generate_gaussian_blobs, "blobs")
produce_plots(generate_gaussian_blobs_different, "blobs_diff")

Learning more

References

[1] van der Maaten, L.J.P.; Hinton, G.E. Visualizing High-Dimensional Data Using t-SNE. Journal of Machine Learning Research 9:2579-2605, 2008.[

[2] L.J.P. van der Maaten. Accelerating t-SNE using Tree-Based Algorithms. Journal of Machine Learning Research 15(Oct):3221-3245, 2014. https://lvdmaaten.github.io/publications/papers/JMLR_2014.pdf

The elements of statistical learning by Trevor Hastie, Robert Tibshirani, Jerome Friedman is a brilliant introduction to the topic and will help you have a better understanding of most of the PCA.

Python data science handbook

Websites

Distill article on tSNE