concrete/benchmarks/glm.py

from copy import deepcopy
from typing import Any, Dict

import numpy as np
import progress
from common import BENCHMARK_CONFIGURATION
from sklearn.compose import ColumnTransformer
from sklearn.datasets import fetch_openml
from sklearn.decomposition import PCA
from sklearn.linear_model import PoissonRegressor
from sklearn.metrics import mean_poisson_deviance
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline, make_pipeline
from sklearn.preprocessing import (
    FunctionTransformer,
    KBinsDiscretizer,
    OneHotEncoder,
    StandardScaler,
)
from tqdm import tqdm

from concrete.quantization import QuantizedArray, QuantizedLinear, QuantizedModule
from concrete.quantization.quantized_activations import QuantizedActivation


class QuantizedExp(QuantizedActivation):
    """
    Quantized Exponential function

    This class will build a quantized lookup table for the exp function
    applied to input calibration data
    """

    def calibrate(self, x: np.ndarray):
        self.q_out = QuantizedArray(self.n_bits, np.exp(x))

    def __call__(self, q_input: QuantizedArray) -> QuantizedArray:
        """Process the forward pass of the exponential.

        Args:
            q_input (QuantizedArray): Quantized input.

        Returns:
            q_out (QuantizedArray): Quantized output.
        """

        quant_exp = np.exp(self.dequant_input(q_input))

        q_out = self.quant_output(quant_exp)
        return q_out


class QuantizedGLM(QuantizedModule):
    """
    Quantized Generalized Linear Model

    Building on top of QuantizedModule, this class will chain together a linear transformation
    and an inverse-link function
    """

    def __init__(self, n_bits, sklearn_model, calibration_data) -> None:
        self.n_bits = n_bits

        # We need to calibrate to a sufficiently low number of bits
        # so that the output of the Linear layer (w . x + b)
        # does not exceed 7 bits
        self.q_calibration_data = QuantizedArray(self.n_bits, calibration_data)

        # Quantize the weights and create the quantized linear layer
        q_weights = QuantizedArray(self.n_bits, np.expand_dims(sklearn_model.coef_, 1))
        q_bias = QuantizedArray(self.n_bits, sklearn_model.intercept_)
        q_layer = QuantizedLinear(self.n_bits, q_weights, q_bias)

        # Store quantized layers
        quant_layers_dict: Dict[str, Any] = {}

        # Calibrate the linear layer and obtain calibration_data for the next layers
        calibration_data = self._calibrate_and_store_layers_activation(
            "linear", q_layer, calibration_data, quant_layers_dict
        )

        # Add the inverse-link for inference.
        # This function needs to be quantized since it's computed in FHE.
        # However, we can use 7 bits of output since, in this case,
        # the result of the inverse-link is not processed by any further layers
        # Seven bits is the maximum precision but this could be lowered to improve speed
        # at the possible expense of higher deviance of the regressor
        q_exp = QuantizedExp(n_bits=7)

        # Now calibrate the inverse-link function with the linear layer's output data
        calibration_data = self._calibrate_and_store_layers_activation(
            "invlink", q_exp, calibration_data, quant_layers_dict
        )

        # Finally construct out Module using the quantized layers
        super().__init__(quant_layers_dict)

    def _calibrate_and_store_layers_activation(
        self, name, q_function, calibration_data, quant_layers_dict
    ):
        # Calibrate the output of the layer
        q_function.calibrate(calibration_data)
        # Store the learned quantized layer
        quant_layers_dict[name] = q_function
        # Create new calibration data (output of the previous layer)
        q_calibration_data = QuantizedArray(self.n_bits, calibration_data)
        # Dequantize to have the value in clear and ready for next calibration
        return q_function(q_calibration_data).dequant()

    def quantize_input(self, x):
        q_input_arr = deepcopy(self.q_calibration_data)
        q_input_arr.update_values(x)
        return q_input_arr


def score_estimator(y_pred, y_gt, gt_weight):
    """Score an estimator on the test set."""

    y_pred = np.squeeze(y_pred)
    # Ignore non-positive predictions, as they are invalid for
    # the Poisson deviance. We want to issue a warning if for some reason
    # (e.g. FHE noise, bad quantization, user error), the regressor predictions are negative

    # Find all strictly positive values
    mask = y_pred > 0
    # If any non-positive values are found, issue a warning
    if (~mask).any():
        n_masked, n_samples = (~mask).sum(), mask.shape[0]
        print(
            "WARNING: Estimator yields invalid, non-positive predictions "
            f" for {n_masked} samples out of {n_samples}. These predictions "
            "are ignored when computing the Poisson deviance."
        )

    # Compute the Poisson Deviance for all valid values
    dev = mean_poisson_deviance(
        y_gt[mask],
        y_pred[mask],
        sample_weight=gt_weight[mask],
    )
    print(f"mean Poisson deviance: {dev}")
    return dev


def score_sklearn_estimator(estimator, df_test):
    """A wrapper to score a sklearn pipeline on a dataframe"""
    return score_estimator(estimator.predict(df_test), df_test["Frequency"], df_test["Exposure"])


def score_concrete_glm_estimator(poisson_glm_pca, q_glm, df_test):
    """A wrapper to score QuantizedGLM on a dataframe, transforming the dataframe using
    a sklearn pipeline
    """
    test_data = poisson_glm_pca["pca"].transform(poisson_glm_pca["preprocessor"].transform(df_test))
    q_test_data = q_glm.quantize_input(test_data)
    y_pred = q_glm.forward_and_dequant(q_test_data)
    return score_estimator(y_pred, df_test["Frequency"], df_test["Exposure"])


@progress.track([{"id": "glm", "name": "Generalized Linear Model", "parameters": {}}])
def main():
    """
    This is our main benchmark function. It gets a dataset, trains a GLM model,
    then trains a GLM model on PCA reduced features, a QuantizedGLM model
    and finally compiles the QuantizedGLM to FHE. All models are evaluated and poisson deviance
    is computed to determine the increase in deviance from quantization and to verify
    that the FHE compiled model acheives the same deviance as the quantized model in the 'clear'
    """

    df, _ = fetch_openml(
        data_id=41214, as_frame=True, cache=True, data_home="~/.cache/sklean", return_X_y=True
    )
    df = df.head(50000)

    df["Frequency"] = df["ClaimNb"] / df["Exposure"]

    log_scale_transformer = make_pipeline(
        FunctionTransformer(np.log, validate=False), StandardScaler()
    )

    linear_model_preprocessor = ColumnTransformer(
        [
            ("passthrough_numeric", "passthrough", ["BonusMalus"]),
            ("binned_numeric", KBinsDiscretizer(n_bins=10), ["VehAge", "DrivAge"]),
            ("log_scaled_numeric", log_scale_transformer, ["Density"]),
            (
                "onehot_categorical",
                OneHotEncoder(sparse=False),
                ["VehBrand", "VehPower", "VehGas", "Region", "Area"],
            ),
        ],
        remainder="drop",
    )

    df_train, df_test = train_test_split(df, test_size=0.2, random_state=0)
    df_calib, df_test = train_test_split(df_test, test_size=100, random_state=0)

    poisson_glm = Pipeline(
        [
            ("preprocessor", linear_model_preprocessor),
            ("regressor", PoissonRegressor(alpha=1e-12, max_iter=300)),
        ]
    )

    poisson_glm_pca = Pipeline(
        [
            ("preprocessor", linear_model_preprocessor),
            ("pca", PCA(n_components=15, whiten=True)),
            ("regressor", PoissonRegressor(alpha=1e-12, max_iter=300)),
        ]
    )

    poisson_glm.fit(df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"])

    poisson_glm_pca.fit(
        df_train, df_train["Frequency"], regressor__sample_weight=df_train["Exposure"]
    )

    # Let's check what prediction performance we lose due to PCA
    print("PoissonRegressor evaluation:")
    _ = score_sklearn_estimator(poisson_glm, df_test)
    print("PoissonRegressor+PCA evaluation:")
    _ = score_sklearn_estimator(poisson_glm_pca, df_test)

    # Now, get calibration data from the held out set
    calib_data = poisson_glm_pca["pca"].transform(
        poisson_glm_pca["preprocessor"].transform(df_calib)
    )

    # Let's see how performance decreases with bit-depth.
    # This is just a test of our quantized model, not in FHE
    for n_bits in [28, 16, 6, 5, 4, 3, 2]:
        q_glm = QuantizedGLM(n_bits, poisson_glm_pca["regressor"], calib_data)
        print(f"{n_bits}b Quantized PoissonRegressor evaluation:")
        score_concrete_glm_estimator(poisson_glm_pca, q_glm, df_test)

    q_glm = QuantizedGLM(2, poisson_glm_pca["regressor"], calib_data)
    dev_pca_quantized = score_concrete_glm_estimator(poisson_glm_pca, q_glm, df_test)
    test_data = poisson_glm_pca["pca"].transform(poisson_glm_pca["preprocessor"].transform(df_test))
    q_test_data = q_glm.quantize_input(test_data)

    engine = q_glm.compile(
        q_test_data,
        BENCHMARK_CONFIGURATION,
        show_mlir=False,
    )

    y_pred_fhe = np.zeros((test_data.shape[0],), np.float32)
    for i, test_sample in enumerate(tqdm(q_test_data.qvalues)):
        with progress.measure(id="evaluation-time-ms", label="Evaluation Time (ms)"):
            q_sample = np.expand_dims(test_sample, 1).transpose([1, 0]).astype(np.uint8)
            q_pred_fhe = engine.run(q_sample)
            y_pred_fhe[i] = q_glm.dequantize_output(q_pred_fhe)

    dev_pca_quantized_fhe = score_estimator(y_pred_fhe, df_test["Frequency"], df_test["Exposure"])

    if dev_pca_quantized_fhe > 0.001:
        difference = abs(dev_pca_quantized - dev_pca_quantized_fhe) * 100 / dev_pca_quantized_fhe
    else:
        difference = 0

    print(f"Quantized deviance: {dev_pca_quantized}")
    progress.measure(
        id="non-homomorphic-loss",
        label="Non Homomorphic Loss",
        value=dev_pca_quantized,
    )

    print(f"FHE Quantized deviance: {dev_pca_quantized_fhe}")
    progress.measure(
        id="homomorphic-loss",
        label="Homomorphic Loss",
        value=dev_pca_quantized_fhe,
    )

    print(f"Percentage difference: {difference}%")
    progress.measure(
        id="relative-loss-difference-percent",
        label="Relative Loss Difference (%)",
        value=difference,
        alert=(">", 7.5),
    )