基于全同态加密的逻辑回归心脏病预测示例详解

1 背景介绍

文章基于TenSEAL开源项目的Tutorial 1内容进行，可以直接打开最后参考1链接查看原文。主要内容是用逻辑回归算法，根据患者的身体数据，预测他未来10年是否有心脏病风险。
预测目标：10 年内患冠心病（CHD）的风险；
输出结果：是（1）或否（0）→ 二分类问题；
使用模型：逻辑回归（Logistic Regression）；
1.1 数据集

该数据集可在Kaggle网站上公开获取，来源于美国马萨诸塞州弗雷明汉镇居民的一项持续进行的心血管研究。研究分类目标是预测患者在未来10年内是否存在冠心病（CHD）风险。数据集包含患者信息，共有4000多条记录和16个属性。

（1）人口统计学（Demographic）
性别（Sex）：男性或女性（分类变量）
年龄（Age）：患者年龄（连续变量——虽然记录为整数，但年龄本质上是连续的）
教育（Education）：教育类型分为4类
（2）行为因素（Behavioral）
当前是否吸烟（Current Smoker）：当前是否为吸烟者（分类变量）
每日吸烟量（Cigs Per Day）：每天平均吸烟数量（可视为连续变量，因为可以是任意数量）
（3）医学史（Medical - history）
是否服用降压药（BP Meds）：是否服用降压药（分类变量）
是否曾中风（Prevalent Stroke）：是否有中风史（分类变量）
是否高血压（Prevalent Hyp）：是否患有高血压（分类变量）
是否糖尿病（Diabetes）：是否患有糖尿病（分类变量）
（4）当前医学指标（Medical - current）
总胆固醇（Tot Chol）：连续变量
收缩压（Sys BP）：连续变量
舒张压（Dia BP）：连续变量
体重指数（BMI）：连续变量
心率（Heart Rate）：连续变量（虽然本质离散，但通常视为连续变量）
血糖（Glucose）：连续变量
（5）预测变量（目标变量）
10 年冠心病风险（CHD）：二分类变量（1 表示“是”，0 表示“否”）。
1.2 数据处理

在原文中提到的P-value（P值）指的是显著性水平，显著性水平是假设检验中的一个概念，是指当原假设为正确时人们却把它拒绝了的概率或风险，这里它用于判断某个特征（如“抽烟”）对结果的影响是否是真实存在的，还是仅仅因为随机巧合。显著性水平是假设检验中的一个概念，是指当原假设为正确时人们却把它拒绝了的概率或风险。它是公认的小概率事件的概率值，必须在每一次统计检验之前确定，通常取α=0.05或α=0.01。在该实例中P值量化了“虚假关联”的风险，即P值α风险高，这个特征可能是“骗子”，可以把它丢弃。文中纸用向后剔除法（Backward Elimination）剔除P值>0.05的特征，最后只保留P值小于0.05的有效特征，如性别：男性患病几率比女性高78.8%，年龄：年龄每增加1岁，风险增加 7%等。 
2 理论基础

2.1 线性回归（Linear Regression)

1 概念

首先看下线性、非线性和回归的概念：
线性：两个变量之间的关系是一次函数关系的——图象是直线，叫做线性。
非线性：两个变量之间的关系不是一次函数关系的——图象不是直线，叫做非线性。
回归：人们在测量事物的时候因为客观条件所限，求得的都是测量值，而不是事物真实的值，为了能够得到真实值，无限次的进行测量，最后通过这些测量数据计算回归到真实值，这就是回归的由来。
所以线性回归就是：用一条直线（或超平面）去拟合数据，通过已知的输入 X，预测一个连续的输出 Y。
2 解决的问题

对大量的观测数据进行处理，从而得到比较符合事物内部规律的数学表达式。也就是说寻找到数据与数据之间的规律所在，从而就可以模拟出结果，也就是对结果进行预测。解决的就是通过已知的数据得到未知的结果，例如：对房价的预测、判断信用评价、电影票房预估等。
3 模型及公式

最简单的线性回归可以用公式：y = wx + b表示，这里y是预测结果，x是输入特征，w是权重（斜率），b是偏置（截距），这就是一条直线，线性回归的任务就是：找到最好的w和b，让这条直线最贴合数据。
在此基础上扩展到多个特征时，有公式：

这时y仍为预测值，x1~xn是n个特征，w1~wn是对用的系数（权重），b仍为截距，该公式还可以表示为更简洁的向量形式：

此时就可通过损失函数来得到预测的损失，如基于MSE（Mean Squared Error）均方误差：

基于RMSE（Root Mean Squared Error）均方根误差：

基于MAE（Mean Absolute Error）平均绝对值误差：

以上m是样本数量，yi是实际值，^yi是预测值，由线性回归公式计算得到，之后即可由梯度下降法，通过对损失函数对各个参数求偏导计算梯度，并沿梯度方向反方向迭代更新参数使得Loss损失最小。
2.2 逻辑回归（Logistic Regression）

逻辑回归是一种统计中的回归分析方法，用于根据一组自变量预测分类因变量的结果。逻辑回归主要用于预测，同时也可以用于计算某事件发生的概率。Logistic回归的因变量可以是二分类的，也可以是多分类的，但是二分类的更为常用，也更加容易解释。逻辑回归常用于数据挖掘，疾病自动诊断，经济预测等领域。逻辑回归和线性回归关系密切，简单来说：逻辑回归 = 线性回归 + Sigmoid激活函数，逻辑回归的核心计算部分，本质就是线性回归；只是在输出层套了一层Sigmoid，把线性输出压缩到0~1的概率区间，从而实现二分类。
Sigmoid函数：

其图形如下：

意义是将连续的(-∞, +∞)映射到[0, 1]，这里[0, 1]正好可以对应概率p取值范围，即p=sigmoid(x) = 1/(1+e-x)，将之前预测结果表示为logit(p) = WTX + b，则有：

同时可知：

即logit(p)正是线性回归的线性部分，同时它是概率的对数几率，把0~1的概率线性映射到整个实数域，它正是线性模型和概率之间的桥梁。在进行逻辑回归时，不直接预测概率，而是预测概率的logit值，再由logit值通过Sigmoid函数还原概率，完美解决了“线性模型预测概率”的问题。对应到心脏病预测实例，通过以下公式完成预测：

2.3 Sigmoid函数近似表示

在全同态加密时，不能简单地在加密数据上计算sigmoid，需要使用低次多项式来近似该，而且受限于同态加密乘法运算深度限制，次数越低越好，所以目标是执行尽可能少的乘法，以便能够使用较小的参数，从而优化计算。参考3中的文章给出两个接近Sigmoid函数的近似函数，这里选取在[-5, 5]范围内更接近Sigmoid函数的σ(x) = 0.5 + 0.197x - 0.004x3。

2.4 Z-Score标准化

目的是将不同量纲（单位）的特征——比如“年龄”（20-80岁）和“胆固醇”（150-400mg/dL）——缩放到同一个“起跑线”上，即将特征取值标准化为均值为0、标准差为1的分布，对应的数学公式是：

xi：特征的取值
μ：特征取值的平均值，特征值减去它能让数据的中心点移动到0；
σ：特征取值的标准差，除以它能让数据的波动范围（缩放比例）变为1；
执行该操作出于以下原因：
（1）梯度爆炸与难以收敛
在训练逻辑回归或神经网络时，如果“胆固醇”数值是 300，而“是否抽烟”是1，模型会给胆固醇分配极小的权重，给抽烟分配极大的权重。这会导致梯度下降时像在“深谷”中反复震荡，很难找到最优解。
（2）特征权重不公平
如果一个特征的数值范围是0-1000，另一个是0-1，模型会下意识地认为数值大的特征更重要。标准化确保了每个特征对预测结果的贡献是基于其变化规律，而不是数值大小。
3 源码分析

3.1 Setup

程序最一开始需要import所有依赖的相关模块，请确保各模块已经按照到系统中：

1. import torch
2. import tenseal as ts
3. import pandas as pd
4. import random
5. from time import time
7. # those are optional and are not necessary for training
8. import numpy as np
9. import matplotlib.pyplot as plt

复制代码

接下来程序会对之前Kaggle网站上下载的程序进行处理，如删除无效数据及无关特征列，并按相同的数量随机采样患病和不患病数据集，这样做的目的是解决类不平衡问题（Class Imbalance）。因为在心脏病预测（Framingham）数据集中，患病（1）的人数远少于未患病（0）的人数，如果不处理，模型会倾向于预测所有人都不患病。程序还提供random_data()函数，该函数生成随机的、线性可分的点，对于那些只想看看事情是如何运作的人来说，可以使用它来代替Kaggle的数据集。这部分程序完整源码如下：

1. torch.random.manual_seed(73)
2. random.seed(73)
4. # 将原始数据进行随机洗牌后，按7:3的比例分割，分别做为训练集和测试集
5. def split_train_test(x, y, test_ratio=0.3):
6.     idxs = [i for i in range(len(x))]
7.     random.shuffle(idxs)
8.     # delimiter between test and train data
9.     delim = int(len(x) * test_ratio)
10.     test_idxs, train_idxs = idxs[:delim], idxs[delim:]
11.     return x[train_idxs], y[train_idxs], x[test_idxs], y[test_idxs]
13. def heart_disease_data():
14.     data = pd.read_csv("./data/framingham.csv")
15.     # drop rows with missing values
16.     data = data.dropna()
17.     # drop some features
18.     data = data.drop(columns=["education", "currentSmoker", "BPMeds", "diabetes", "diaBP", "BMI"])
19.     print("save cleaned data to data/framingham_cleaned.csv")
20.     data.to_csv("./data/framingham_cleaned.csv", index=False)
21.     # 根据TenYearCHD（十年内心脏病风险，0或1）这一列，将原始数据集拆分为两个子集
22.     grouped = data.groupby('TenYearCHD')
23.     #data = grouped.apply(lambda x: x.sample(grouped.size().min(), random_state=73).reset_index(drop=True))
24.     # 1. 执行分组并采样，使用数量少的那个类别个数进行采样，每个组（0 和 1）都随机抽取等量的样本
25.     data = grouped.apply(lambda x: x.sample(grouped.size().min(), random_state=73), include_groups=False)
26.     #data.to_csv("./data/twogroupdata1.csv")
27.     # 2. 核心步骤：恢复被排除的分组列
28.     # 因为 include_groups=False 把它变成了索引，我们需要把它变回列，使数据结构回归到普通的表格形式
29.     data = data.reset_index(level=0).reset_index(drop=True)
30.     #data.to_csv("./data/twogroupdata2.csv")
32.     # extract labels
33.     y = torch.tensor(data["TenYearCHD"].values).float().unsqueeze(1)
34.     # 丢弃TenYearCHD列
35.     data = data.drop(columns="TenYearCHD")
36.     # standardize data，标准化确保了每个特征对预测结果的贡献是基于其变化规律，而不是数值大小
37.     data = (data - data.mean()) / data.std()
38.     #data.to_csv("./data/last.csv")
39.     #print("have save data/last.csv")
40.     x = torch.tensor(data.values).float()
41.     return split_train_test(x, y)
44. def random_data(m=1024, n=2):
45.     # data separable by the line `y = x`
46.     x_train = torch.randn(m, n)
47.     x_test = torch.randn(m // 2, n)
48.     y_train = (x_train[:, 0] >= x_train[:, 1]).float().unsqueeze(0).t()
49.     y_test = (x_test[:, 0] >= x_test[:, 1]).float().unsqueeze(0).t()
50.     return x_train, y_train, x_test, y_test
53. # You can use whatever data you want without modification to the tutorial
54. # x_train, y_train, x_test, y_test = random_data()
55. x_train, y_train, x_test, y_test = heart_disease_data()
57. print("############# Data summary #############")
58. print(f"x_train has shape: {x_train.shape}")
59. print(f"y_train has shape: {y_train.shape}")
60. print(f"x_test has shape: {x_test.shape}")
61. print(f"y_test has shape: {y_test.shape}")
62. print("#######################################")

复制代码

程序运行后输出如下：

3.2 训练逻辑回归模型

将首先训练一个逻辑回归模型（没有任何加密），它可以被视为一个具有单个节点的单层神经网络，后续将使用此模型作为与加密训练和评估进行比较的手段。

1. class LR(torch.nn.Module):
2.     def __init__(self, n_features):
3.         super(LR, self).__init__()
4.         self.lr = torch.nn.Linear(n_features, 1)
5.         
6.     def forward(self, x):
7.         out = torch.sigmoid(self.lr(x))
8.         return out
10. n_features = x_train.shape[1]
11. model = LR(n_features)
12. # use gradient descent with a learning_rate=1
13. optim = torch.optim.SGD(model.parameters(), lr=1)
14. # use Binary Cross Entropy Loss
15. criterion = torch.nn.BCELoss()
17. # define the number of epochs for both plain and encrypted training
18. EPOCHS = 5
20. def train(model, optim, criterion, x, y, epochs=EPOCHS):
21.     for e in range(1, epochs + 1):
22.         optim.zero_grad()
23.         out = model(x)
24.         loss = criterion(out, y)
25.         loss.backward()
26.         optim.step()
27.         print(f"Loss at epoch {e}: {loss.data}")
28.     return model
30. model = train(model, optim, criterion, x_train, y_train)
32. def accuracy(model, x, y):
33.     out = model(x)
34.     correct = torch.abs(y - out) < 0.5
35.     return correct.float().mean()
37. plain_accuracy = accuracy(model, x_test, y_test)
38. print(f"Accuracy on plain test_set: {plain_accuracy}")

复制代码

该段程序输出如下：

正如原文所说，高精度不是该程序的目标，这里只是想看看加密数据的训练不会影响最终结果，所以将比较加密数据的准确性和在这里得到的plain_accuracy。
3.3 基于明文参数的加密数据评估

在这一部分中，将只关注在加密测试集上使用明文参数（可选加密参数）评估逻辑回归模型。首先创建一个类似PyTorch的LR模型，可以评估加密数据：

1. class EncryptedLR:
2.     def __init__(self, torch_lr):
3.         # TenSEAL processes lists and not torch tensors,
4.         # so we take out the parameters from the PyTorch model
5.         self.weight = torch_lr.lr.weight.data.tolist()[0]
6.         self.bias = torch_lr.lr.bias.data.tolist()
7.         
8.     def forward(self, enc_x):
9.         # We don't need to perform sigmoid as this model
10.         # will only be used for evaluation, and the label
11.         # can be deduced without applying sigmoid
12.         enc_out = enc_x.dot(self.weight) + self.bias
13.         return enc_out
14.    
15.     def __call__(self, *args, **kwargs):
16.         return self.forward(*args, **kwargs)
17.         
18.     ################################################
19.     ## You can use the functions below to perform ##
20.     ## the evaluation with an encrypted model     ##
21.     ################################################
22.    
23.     def encrypt(self, context):
24.         self.weight = ts.ckks_vector(context, self.weight)
25.         self.bias = ts.ckks_vector(context, self.bias)
26.         
27.     def decrypt(self, context):
28.         self.weight = self.weight.decrypt()
29.         self.bias = self.bias.decrypt()
31. eelr = EncryptedLR(model)

复制代码

之后创建一个TenSEAL Context，用于指定要使用的方案和参数。在这里，选择小而安全的参数，允许进行一次乘法。这足以评估逻辑回归模型，之后会发现，在对加密数据进行训练时，会需要更大的参数。

1. # parameters
2. poly_mod_degree = 4096
3. coeff_mod_bit_sizes = [40, 20, 40]
4. # create TenSEALContext
5. ctx_eval = ts.context(ts.SCHEME_TYPE.CKKS, poly_mod_degree, -1, coeff_mod_bit_sizes)
6. # scale of ciphertext to use
7. ctx_eval.global_scale = 2 ** 20
8. # this key is needed for doing dot-product operations
9. ctx_eval.generate_galois_keys()
11. t_start = time()
12. enc_x_test = [ts.ckks_vector(ctx_eval, x.tolist()) for x in x_test]
13. t_end = time()
14. print(f"Encryption of the test-set took {int(t_end - t_start)} seconds")

复制代码

在代码中，会在评估前加密测试数据集。接下来在构建EncryptedLR类时，不会在线性层的加密输出上计算sigmoid函数，仅仅是因为它不是必需的，在加密数据上计算sigmic会增加计算时间并需要更大的加密参数，但是，在之后的加密训练中会使用sigmoid。当前，直接进行加密测试集上的评估，并将其准确性与普通测试集进行比较。

1. def encrypted_evaluation(model, enc_x_test, y_test):
2.     t_start = time()
3.    
4.     correct = 0
5.     for enc_x, y in zip(enc_x_test, y_test):
6.         # encrypted evaluation
7.         enc_out = model(enc_x)
8.         # plain comparison
9.         out = enc_out.decrypt()
10.         out = torch.tensor(out)
11.         out = torch.sigmoid(out)
12.         if torch.abs(out - y) < 0.5:
13.             correct += 1
14.    
15.     t_end = time()
16.     print(f"Evaluated test_set of {len(x_test)} entries in {int(t_end - t_start)} seconds")
17.     print(f"Accuracy: {correct}/{len(x_test)} = {correct / len(x_test)}")
18.     return correct / len(x_test)
19.    
21. encrypted_accuracy = encrypted_evaluation(eelr, enc_x_test, y_test)
22. diff_accuracy = plain_accuracy - encrypted_accuracy
23. print(f"Difference between plain and encrypted accuracies: {diff_accuracy}")
24. if diff_accuracy < 0:
25.     print("Oh! We got a better accuracy on the encrypted test-set! The noise was on our side...")

复制代码

程序运行结果如下：

不仅比明文直接评估的精度有所下降，而且对比原文章加密评估的精度0.6736526946107785，还要低一些只有0.6167664670658682，这里不知是什么原因。
3.4 基于加密数据训练的加密逻辑回归模型

在这一部分中，将重新定义一个类似PyTorch的模型，该模型既可以向前传播加密数据，也可以反向传播以更新权重，从而在加密数据上训练加密的逻辑回归模型，以下是关于训练的更多细节。
1 损失函数

这里使用带有正则化的二元交叉熵损失函数，y(i)是第i个预期标签，^y(i)是是逻辑回归模型的第i个输出，θ是n维权重向量，损失函数如下：

上面公式中m是样本数量，n是特征数量，损失可以分成两部分，前半部分是交叉熵损失：

当真实值y=0和y=1时，以上公式可分别简化为：

以真实值y=0为例，预测^y越接近0，则损失值越小越接近零，反之预测越接近1，损失值越大，同样对于真实值y=1时，有相同的规律，而且通过响应函数的图形也能直观的看到该规律：

公式中的后半部分是正则化项（Regularization）：

θj：权重向量中的第j个值，λ：正则化系数，该部分意义在于惩罚过大的权重，以保证模型不能过度依赖某一个特征，必须保持权重相对“温和”。
2 参数更新

为了进行参数更新，使用如下规则，这里x(i)是第i个输入数据：

然而，由于同态加密约束，这里使用α=1以减少乘法，并使用λ/m=0.05，从而得出以下更新规则：

3 同态加密参数


从输入数据到参数更新，密文需要深度为6的乘法运算，1用于点积运算，2用于sigmoid算法近似计算，3用于反向传播阶段（其中1个隐藏在backward函数中的self._delta_w += enc_x * out_minus_y运算中，该运算将1维度的向量与n维度的向量相乘，需要掩码提取第一个槽位的值并依次复制到n个槽位的其他位置）。对于大约20位的缩放，我们需要6个与缩放具有相同比特大小的系数模数，加上最后一个需要更多比特的系数，我们已经超出了4096个多项式模数（如果我们考虑128位的安全性，这需要系数模数的总比特数