ml-sdk/use_case_examples/titanic/KaggleTitanic.ipynb

Privacy-Preserving Machine Learning on Titanic

This notebook introduces a Privacy-Preserving Machine Learning (PPML) solution to the Kaggle Titanic competition using the TorusML open-source framework. Its main ambition is to show that Fully Homomorphic Encryption (FHE) can be used for protecting data when using a Machine Learning model to predict outcomes without degrading its performance. In this example, a XGBoost classifier model will be considered as it achieves near state-of-the-art accuracy.

It also took some ideas from several upvoted public notebooks, including this one from Manav Sehgal and this one from LD Freeman.

Imports

In [1]:

import time

import numpy as np
import pandas as pd
from sklearn.datasets import fetch_openml
from sklearn.model_selection import GridSearchCV, ShuffleSplit, train_test_split
from xgboost import XGBClassifier

from concrete.ml.sklearn import XGBClassifier as ConcreteXGBClassifier

No CUDA runtime is found, using CUDA_HOME='/usr/local/cuda'

In [2]:

# pylint: disable-next=no-member
all_data = fetch_openml(data_id=40945, as_frame=True, cache=True).frame

all_data = all_data.convert_dtypes()
all_data["survived"] = all_data["survived"].astype(int)
print(all_data.columns)

Index(['pclass', 'survived', 'name', 'sex', 'age', 'sibsp', 'parch', 'ticket',
       'fare', 'cabin', 'embarked', 'boat', 'body', 'home.dest'],
      dtype='object')

Preprocessing the Data

Be sure to launch the download_data.sh script in order to have local versions of the data-set.

In [3]:

train_data, test_data = train_test_split(all_data, test_size=len(all_data.index) - 891)
datasets = [train_data, test_data]

Let's take a closer look at the train data:

In [4]:

train_data.head()

Out [4]:

	pclass	survived	name	sex	age	sibsp	parch	ticket	fare	cabin	embarked	boat	body	home.dest
995	3	0.0	Markoff, Mr. Marin	male	35.0	0	0	349213	7.8958	<NA>	C	<NA>	<NA>	<NA>
1147	3	0.0	Riihivouri, Miss. Susanna Juhantytar 'Sanni'	female	22.0	0	0	3101295	39.6875	<NA>	S	<NA>	<NA>	<NA>
485	2	0.0	Levy, Mr. Rene Jacques	male	36.0	0	0	SC/Paris 2163	12.875	D	C	<NA>	<NA>	Montreal, PQ
1265	3	0.0	Van Impe, Miss. Catharina	female	10.0	0	2	345773	24.15	<NA>	S	<NA>	<NA>	<NA>
632	3	0.0	Andersson, Mrs. Anders Johan (Alfrida Konstant...	female	39.0	1	5	347082	31.275	<NA>	S	<NA>	<NA>	Sweden Winnipeg, MN

We can observe:

• the target variable: Survived
• some numerical variables: Pclass, SbSp, Parch, Fare
• some categorical (non-numerical) variables: Name, Sex, Ticket, Cabin, Embarked

Missing Values

Then, we can notice that some values are missing for the Cabin variable. We must therefore investigate a bit more about this by printing the total amounts of missing values for each variables:

In [6]:

print("-" * 3, "Train data", "-" * 3)
print(train_data.isnull().sum(), "\n")
print("-" * 3, "Test data", "-" * 3)
print(test_data.isnull().sum())

--- Train data ---
pclass         0
survived       0
name           0
sex            0
age          173
sibsp          0
parch          0
ticket         0
fare           0
cabin        691
embarked       1
boat         556
body         806
home.dest    378
dtype: int64 

--- Test data ---
pclass         0
survived       0
name           0
sex            0
age           90
sibsp          0
parch          0
ticket         0
fare           1
cabin        323
embarked       1
boat         267
body         382
home.dest    186
dtype: int64

Only four variables are incomplete: Cabin, Age, Embarked and Fare. However, the Cabin one seems to be missing quite more data than the others:

In [7]:

for incomp_var in train_data.columns:
    missing_val = pd.concat(datasets)[incomp_var].isnull().sum()
    if missing_val > 0 and incomp_var != "survived":
        total_val = pd.concat(datasets).shape[0]
        print(f"Percentage of missing values in {incomp_var}: {missing_val/total_val*100:.1f}%")

Percentage of missing values in age: 20.1%
Percentage of missing values in fare: 0.1%
Percentage of missing values in cabin: 77.5%
Percentage of missing values in embarked: 0.2%
Percentage of missing values in boat: 62.9%
Percentage of missing values in body: 90.8%
Percentage of missing values in home.dest: 43.1%

Since the Cabin variable misses more than 2/3 of its values, it might not be relevant to keep it:

In [8]:

drop_column = ["cabin", "boat", "body", "home.dest", "ticket"]

for dataset in datasets:
    # Remove irrelevant variables
    dataset.drop(drop_column, axis=1, inplace=True)

For the other ones, we can replace the missing values using:

• the median value for Age and Fare
• the most common value for Embarked

In [9]:

for dataset in datasets:
    # Complete missing Age values with median
    dataset.age.fillna(dataset.age.median(), inplace=True)

    # Complete missing Embarked values with the most common one
    dataset.embarked.fillna(dataset.embarked.mode()[0], inplace=True)

    # Complete missing Fare values with median
    dataset.fare.fillna(dataset.fare.median(), inplace=True)

Feature Engineering

We can manually extract and create new features in order to help the model interpret some behaviors and correctly predict an outcome. Those new features are:

• FamilySize: The size of the family the individual was traveling with, with 1 being someone that traveled alone.
• IsAlone: A boolean variable stating if the individual was traveling alone (1) or not (0). This might help the model to emphasize on this idea of traveling with relatives or not.
• Title: The individual's title (Mr, Mrs, ...), often indicating a certain social status.
• Farebin and AgeBin: Binned version of the Fare and Age variables. It groups values together, generally reducing the impact of minor observation errors.

In [10]:

def get_bin_labels(bin_name, number_of_bins):
    labels = []
    for i in range(number_of_bins):
        labels.append(bin_name + f"_{i}")
    return labels


for dataset in datasets:
    # Emphasize on relatives
    dataset["FamilySize"] = dataset.sibsp + dataset.parch + 1

    dataset["IsAlone"] = 1
    dataset.IsAlone[dataset.FamilySize > 1] = 0

    # Consider an individual's social status
    dataset["Title"] = dataset.name.str.extract(r" ([A-Za-z]+)\.", expand=False)

    # Group fares and ages in "bins" or "buckets"
    dataset["FareBin"] = pd.qcut(dataset.fare, 4, labels=get_bin_labels("FareBin", 4))
    dataset["AgeBin"] = pd.cut(dataset.age.astype(int), 5, labels=get_bin_labels("AgeBin", 5))

    # Remove now-irrelevant variables
    drop_column = ["name", "sibsp", "parch", "fare", "age"]
    dataset.drop(drop_column, axis=1, inplace=True)

Let's have a look on the titles' distribution:

In [11]:

data = pd.concat(datasets)
titles = data.Title.value_counts()
print(titles)

Title
Mr          757
Miss        260
Mrs         197
Master       61
Dr            8
Rev           8
Col           4
Ms            2
Major         2
Mlle          2
Lady          1
Mme           1
Countess      1
Dona          1
Sir           1
Jonkheer      1
Capt          1
Don           1
Name: count, dtype: Int64

We can clearly observe that only a few titles represent most of the individuals. In order to prevent the model from becoming overly specific, we decide to group all the "uncommon" titles together:

In [12]:

uncommon_titles = titles[titles < 10].index

for dataset in datasets:
    dataset.Title.replace(uncommon_titles, "Rare", inplace=True)

Dummification

Finally, we can "dummify" the remaining categorical variables. Dummification is a common technique of transforming categorical (non-numerical) data into numerical data without having to map values or consider any order between each of them. The idea is to take all the different values found in a variable and create a new associated binary variable.

For example, the "Embarked" variable has three categorical values: "S", "C" and "Q". Dummifying the data will create three new variables, "Embarked_S", "Embarked_C" and "Embarked_Q", and set the value of "Embarked_S" (resp. "Embarked_C" and "Embarked_Q") to 1 for each data point initially labeled with "S" (resp. "C" and "Q") in the "Embarked" variable, else 0.

In [13]:

categorical_features = train_data.select_dtypes(exclude=np.number).columns
x_train = pd.get_dummies(train_data, prefix=categorical_features)
x_test = pd.get_dummies(test_data, prefix=categorical_features)

We then split the target variable from the others.

In [14]:

target = "survived"
x_train = x_train.drop(columns=[target])
y_train = train_data[target]

x_test = x_test.drop(columns=[target])
y_test = test_data[target]

Training

Training with XGBoost

Let's first train a classifier model using XGBoost. Since several parameters have to be fixed beforehand, we use scikit-learn's GridSearchCV method to perform cross validation in order to maximize our chance to find the best ones.

In [15]:

# Instantiate the Cross-Validation generator
cv = ShuffleSplit(n_splits=5, test_size=0.3, random_state=0)

# Set the parameters to tune.
# Those ranges are voluntarily small in order to keep the FHE execution time per inference
# relatively low. In fact, we found out that, in this particular Titanic example, models with
# larger numbers of estimators or maximum depth don't score a much better accuracy.
param_grid = {
    "max_depth": list(range(1, 5)),
    "n_estimators": list(range(1, 5)),
    "learning_rate": [0.01, 0.1, 1],
}

# Instantiate and fit the model through grid-search cross-validation
time_begin = time.time()
model = GridSearchCV(XGBClassifier(), param_grid, cv=cv, scoring="roc_auc")
model.fit(x_train, y_train)
cv_xgb_duration = time.time() - time_begin

print(f"Best hyper-parameters found in {cv_xgb_duration:.2f}s :", model.best_params_)

Best hyper-parameters found in 7.98s : {'learning_rate': 1, 'max_depth': 2, 'n_estimators': 4}

Training with TorusML

Now, let's do the same with TorusML's XGBClassifier method.

In order to do so, we need to specify the number of bits over which inputs, outputs and weights will be quantized. This value can influence the precision of the model as well as its inference running time, and therefore can lead the grid-search cross-validation to find a different set of parameters. In our case, setting this value to 2 bits outputs an excellent accuracy score while running faster.

In [16]:

# The TorusML model needs an additional parameter used for quantization
param_grid["n_bits"] = [2]

# Instantiate and fit the model through grid-search cross-validation
time_begin = time.time()
concrete_model = GridSearchCV(ConcreteXGBClassifier(), param_grid, cv=cv, scoring="roc_auc")
concrete_model.fit(x_train, y_train)
cv_concrete_duration = time.time() - time_begin

print(f"Best hyper-parameters found in {cv_concrete_duration:.2f}s :", concrete_model.best_params_)

Best hyper-parameters found in 40.99s : {'learning_rate': 0.1, 'max_depth': 4, 'n_bits': 2, 'n_estimators': 4}

Predicting the Outcomes

Computing the predictions in FHE on the complete test set of 418 data points using the above hyper-parameters can take up to 1 minute, using a c5.4xlarge AWS instance.

In [17]:

# Compute the predictions in clear using XGBoost
clear_predictions = model.predict(x_test)

# Compute the predictions in clear using TorusML
clear_quantized_predictions = concrete_model.predict(x_test)

# Compile the TorusML model on a subset
fhe_circuit = concrete_model.best_estimator_.compile(x_train.head(100))

# Generate the keys
# This step is not absolutely necessary, as keygen() is called, when necessary,
# within the predict method.
# However, it is useful to run it beforehand in order to be able to
# measure the prediction executing time separately from the key generation one
time_begin = time.time()
fhe_circuit.keygen()
key_generation_duration = time.time() - time_begin

# Compute the predictions in FHE using TorusML
time_begin = time.time()
fhe_predictions = concrete_model.best_estimator_.predict(x_test, fhe="execute")
prediction_duration = time.time() - time_begin

print(f"Key generation time: {key_generation_duration:.2f}s")
print(f"Total execution time for {len(clear_predictions)} inferences: {prediction_duration:.2f}s")
print(f"Execution time per inference in FHE: {prediction_duration / len(clear_predictions):.2f}s")

Key generation time: 0.94s
Total execution time for 418 inferences: 65.49s
Execution time per inference in FHE: 0.16s

FHE computations are expected to be exact. This means that the model executed in FHE results in the same predictions as the TorusML one, which is executed in clear and only considers quantization.

In [18]:

number_of_equal_preds = np.sum(fhe_predictions == clear_quantized_predictions)
pred_similarity = number_of_equal_preds / len(clear_predictions) * 100
print(
    "Prediction similarity between both TorusML models (quantized clear and FHE): "
    f"{pred_similarity:.2f}%",
)

accuracy_fhe = np.mean(fhe_predictions == y_test) * 100
print(
    "Accuracy of prediction in FHE on the test set " f"{accuracy_fhe:.2f}%",
)

Prediction similarity between both TorusML models (quantized clear and FHE): 100.00%
Accuracy of prediction in FHE on the test set 78.23%

However, as seen previously, the grid-search cross-validation was done separately between the XGBoost model and the TorusML one. For this reason, the two models do not share the same set of hyper-parameters, making their decision boundaries different.