第5章基于矩阵机制的差分隐私
连续数据发布
5.1引言
随着数字技术的发展，数据越来越多地出现于现实生活当中。数据给人们的生活带来的好处不言而喻。人们不仅可以利用数据进行评估、分析和预测，还可以从中寻找有价值的结论，如“啤酒与尿布”的故事。然而，在享受数据带来的好处的同时，也应该注意到数据中包含的个人隐私信息可能存在泄露的风险。特别是当攻击者怀有恶意时，他就有可能利用已掌握的知识分析用户所发布的数据，并从中挖掘出数据所对应的用户的隐私信息。例如，只需根据4个时空点就能使95%的人泄露其位置信息［1］。因此，如何在发布数据的同时避免数据中包含的隐私被泄露是数据时代亟待解决的问题之一。针对这一问题，各种隐私保护模型被提出。其中，以提供严格数据保护为特点的差分隐私模型［24］得到了广泛的认可。该模型被提出后，人们基于该模型开展了很多研究工作。内容涉及直方图发布［58］、空间划分发布［9,10］、智能数据分析［11,12］等，有效克服了该模型基于k匿名［13］和划分［14］的隐私保护方法需要事先对攻击做出假设的不足。差分隐私数据发布研究的关键问题在于如何在保证差分隐私的前提下提高发布数据的可用性。
现有关于差分隐私的数据发布方法大多关注静态发布问题，而现实应用中更多情况下需要发布方法具有连续数据发布的能力。然而，研究表明，这些方法无法应用于连续数据发布问题。为此，本章对差分隐私下的连续数据发布问题展开研究。例如，某医疗数据库中记录了每个月的入院病人的信息，其中病人感染HIV的情况为敏感信息。表5.1展示了其中3个月的数据示例。同时，出于某研究目的，医院将按月统计并公布入院的HIV病人数。公布的数据形如表5.2所示，医院累计当前入院并且感染HIV的病人并于当月发布最新数据。与数据静态发布不同，在医院发布完每个月的统计信息后，该数据并非不再改变，而是在下个月将得到更新。更重要的是，在发布每一次数据的过程中，以后需要发布的数据是无法被预知的。该问题的核心是在满足差分隐私的条件下，寻找更精确、更高效的连续数据发布方法。表5.1病人感染HIV情况NameHIV+NameHIV+NameHIV+AliceYesAlanNoAndyNoBobNoBenYesBillNoCarolYesCariNoChenYes︙︙︙︙︙︙表5.2入院的HIV病人数MonthHIV+Sum of HIV+1531 5312392 92334261349︙︙︙以上连续发布问题的一种朴素解决方案［15］是直接将前一个月发布的HIV病人数与本月新增的HIV病人数相加，然后再添加噪声使其满足差分隐私。该方案导致每一次发布数据的噪声的均方误差线性累加，最终使发布的数据失去可用性。文献［15］针对该问题提出了一种基于二叉树的发布方法。然而，此方法仅仅引入二叉树模拟发布，并未对精确性提出有效优化方法。为此，本章以提高发布数据的精确性为主要目标，将矩阵机制引入差分隐私连续数据发布问题中，以期设计出高效的基于矩阵机制的差分隐私连续数据发布方法，可有效满足大规模连续数据发布的要求。
5.2基础知识与问题提出
在差分隐私的数据发布中，为提高数据发布的精度，Hay等人［7］和Xiao等人［8］分别提出基于一致性调节的区间树方法和小波变换方法，实现了较高精度的数据发布。然而，上述两种方法只适用于差分隐私下的数据静态发布，无法应用于差分隐私下的连续数据发布。Chan等人［15］提出了两种利用二叉树结构进行连续数据发布的方法。第一种方法是构建一棵叶节点数量为2m的完全二叉树，然后利用模拟二叉树统计发布的过程进行连续数据发布。第二种方法是第一种方法的改进版本，它试图通过调整二叉树各层节点的隐私预算分配来达到无限发布的效果。研究表明，虽然第一种方法相比于朴素方法，数据发布的精确性有显著的提升，然而该方法仅仅引入二叉树结构来模拟发布过程，并未做进一步改进，因此数据发布的精确性仍有较大的提升空间；第二种方法的隐私预算分配并不合理，导致发布数据的误差远大于第一种方法。
为了解决差分隐私下的线性查询问题，Li等人［16］提出了基于矩阵机制的批量查询方法。其基本思想是通过寻找策略矩阵对线性查询进行优化，进而提高发布数据的精确性。然而，该文献提出的矩阵机制仅能满足小规模数据集和查询负载的要求。此外，它还很容易产生次优化的查询策略，使得结果往往并不理想。为此，Yuan等人［17］利用了负载矩阵低秩的性质进行优化，提出低秩矩阵机制，在一定程度上改善了原有矩阵机制的不足，提升了数据发布效率与精确性。然而，该文献提出的优化查询使用半正定规划算法，同样只能适用于小规模数据集和查询负载的要求。由于本章研究的连续数据发布问题本质上也是线性查询问题，因此拟利用矩阵机制，结合连续数据发布问题本身具有的一些特性，设计出精确性更高的高效算法，使之具有大规模连续数据发布的能力。
矩阵机制是一种针对差分隐私下线性查询问题的优化方法。它通过将查询集Q转换成负载矩阵W，然后寻找最优策略矩阵M来实现差分隐私下线性查询的优化。其中，查询集Q是一组线性查询的集合，满足Q={q1,q2,…,qn}。每个线性查询表示如下: q=∑m1wixi其中,X=(x1,x2,…,xm)T为数据向量，wi为该查询在分量xi上的权重。负载矩阵W由每组线性查询的权重组成，并满足Q=WX=∑mj=1W1jxj,∑mj=1W2jxj,…,∑mj=1WnjxjT原始的矩阵机制［16］通过直接寻找策略矩阵M的方式求解问题，这种做法的效率和优化效果都不够理想。而低秩矩阵机制［17］则是采用分解负载矩阵的方法来寻找优化策略。在该机制下，它将W分解成两个矩阵B、M。其中M表示低秩矩阵下的策略矩阵，且满足W=BM。通过对中间结果MX添加噪声的方式减少误差。其形式化表示如下: A(W,X)=BMX+ΔMεn(5.1)文献［17］指出，式(5.1)的敏感度ΔM与策略矩阵的列范式相等，即ΔM=M1。而低秩矩阵的均方误差由下式求得: 2ε2trace(BTB)ΔM2(5.2)研究表明，差分隐私下的数据连续发布问题能够被转换成基于矩阵机制的优化问题。只需将每一次发布视为一个查询，然后将所有发布过程视为查询负载，并转换成相应的矩阵，利用矩阵机制进行求解。
5.3基于矩阵机制的差分隐私连续数据发布
考虑一个随着时间增长会不断产生记录并被添加进来的记录流。该记录流是数据发布的来源。记录流的每条记录都有需要保护的属性σ，满足σ∈{0,1},并假设记录集 i表示第i－1次和第i次发布之间记录流所中被添加的记录的集合。由 i可求出第i次发布的数据增量ai=|{σ|σ∈ i且σ=1}|。
定义5.1（连续数据发布）对于记录流，数据发布者随着记录流中的记录增长按照某种规则多次发布当前记录流中满足σ=1的记录数的行为即为连续数据发布。假设第i次发布的累计数据为si，那么si满足si=si－1+ai=∑ij=1aj(5.3)差分隐私下的连续数据发布行为即通过某种隐私算法A()发布添加了噪声的累计数据i，从而使数据连续发布的结果满足ε差分隐私。同时，在此基础上，本节还要求提出的算法能够精确而且高效地发布数据。
为了避免数据连续发布的敏感度过高而影响数据发布精度的问题，本节主要考虑了满足ε差分隐私的次数受限的数据连续发布算法。
定义5.2（次数受限的数据连续发布算法）如果隐私算法A()至多接受并发布N次满足ε差分隐私的统计数据，则称该算法为次数受限的数据连续发布算法。
上述问题能够转换成线性查询问题，而矩阵机制能够对线性查询问题进行优化。为了提出更精确的数据连续发布算法，本章将这一问题与矩阵机制相结合，并在此基础上寻找快速发布算法。而且，由于矩阵机制是成熟的且经过严格检验的满足ε差分隐私［17］的隐私保护机制。因此，基于矩阵机制的线性查询算法只要符合式(5.1)的形式就保证了该算法满足ε差分隐私。
利用矩阵机制优化前，需要将式(5.3)转换成负载矩阵W: W=100…
110…
111…
︙︙︙(5.4)可以看出，数据连续发布下的负载矩阵W为下三角矩阵，且满足Wij=1,i<j。
根据矩阵机制的特征及一般研究过程［16,17］。本章将通过以下步骤对该问题展开研究: 
(1) 寻找初始的策略矩阵M，将W分解为矩阵B、M，使矩阵机制能够较为精确地发布数据。
(2) 对策略矩阵M进行研究，寻找优化策略以优化数据发布的精确性。
(3) 综合分析矩阵B、M以及优化策略的性质，在保证数据发布的精确性得到优化的前提下，提出高效的优化算法。
5.4隐私连续数据发布算法〖*4/5〗5.4.1策略矩阵的构建由于数据随着发布过程动态产生，未来数据无法得知，因此只能根据当前以及过往的数据进行优化。这一特征反映到矩阵机制时，就要求矩阵机制所应用的策略矩阵为下三角矩阵。通过各种数据结构的对比研究，发现树状数组［19］更加适合构造基于矩阵机制的数据连续发布方法的策略矩阵。它能够自然并且快速求解数列的前k项和，符合连续数据发布的基本特征。同时，初步研究表明，将它与差分隐私结合而提出的隐私保护模型能够达到与基于二叉树的连续数据发布方法［15］相当的精确性。同时，深入研究表明，结合树状数组的差分隐私模型在实现的巧妙性以及进一步优化的潜力方面均比后者略胜一筹。
树状数组主要针对以下问题: 给定N个实数，记为a1,a2,…,an，要求快速求出前k项的和。记前k项的和为sk，则sk=∑kj=1aj。
针对该问题，利用树状数组提出如下解决方法。该方法计算了中间统计量ci。而ci由以下公式求得:ci=∑ij=i－lowbit(i)+1aj(5.5)其中，函数lowbit(x)表示将正整数x写成二进制形式后，将该二进制数值为1的最低位的数置为1，其余位均置为0。例如，当x=10时，其对应的二进制数为(1010)2。从低往高，20位为0，21位为1，那么，lowbit(x)的输出值就将21位置为1，其他位均置为0，即(0010)2=2。再将其转成十进制数，就得到lowbit(10)=2。
该函数的算法如下。算法5.1lowbit(x)
输入: 正整数x
输出: x对应的lowbit值
1. p←0，y←1;
2. while x mod 2=0
3. y←2y;x←x2;
4. wend
5. return y结合算法5.1以及式(5.5)，按照树状数组求解中间统计量的方式构造策略矩阵M,算法如下。算法5.2求解策略矩阵M
输入: 发布次数N
输出: 策略矩阵M
1. M←0N×N;//初始化为零矩阵
2. for p=1 to N do
3. pt←p;
4. while pt<N
5. Mpt,p←1;//更新矩阵元素
6. pt←pt+lowbit(pt);
7. wend
8. end for
9. return M下面讨论矩阵B的求解。由W=BM以及MN的可逆性可得B=WNM－1N。因此，只需根据树状数组的性质就能够快速地求解B。而树状数组的求和操作按照以下公式进行求解: 
si=ci+si－lowbit(i)(5.6)算法5.3求解矩阵B
输入: 发布次数N
输出: 矩阵B
1. B←0N×N;//初始化为零矩阵
2. for p=1 to N do 
3. pt←p;
4. while pt>0
5.Bpt,p←1//更新矩阵元素
6.pt←ptlowbit(pt);
7. wend
8. end for
9. return B上述算法结合低秩矩阵的表达式，即可得到基于树状数组的数据连续发布的表达式: A(W,D)=BMNX+ΔMNεn(5.7) 接下来分析该策略矩阵所产生的均方误差的情况。关于MN和B有如下两个相关定理。
定理5.1‖MN‖1=‖MN(:，1)‖1≥‖MN(:，j)‖1(j>1),且‖MN‖1= log2N +1。
证明: 通过算法5.2研究矩阵MN(:，1)的构造情况。可得当第一次迭代(p=1)时，有pt=p=1=20。设第t次迭代时有pt=2t－1，由更新表达式有pt=pt－lowbit(pt)=2t－1+2t－1=2t。而根据步骤6的判断条件，pt>N时，该列构造结束。因此有，2t－1≤Nt≤log2N+1。此时Mpt,p=1更新了 log2N +1次，因此‖MN(:，1)‖1= log2N +1(j>1)。
对于j>1的情况，假设第t次迭代时pt>2t－1且lowbit(pt)≥2t－1，由第一次迭代有pt=j>1可知满足该条件。根据lowbit函数的性质，有lowbit(pt′)=lowbit(pt+lowbit(pt))≥lowbit(2lowbit(pt))=2lowbit(pt)≥2t，从而pt′=pt+lowbit(pt)>2t。很显然，根据pt>2t－1可推得，j>1时的更新次数不大于j=1时。因此，‖MN(:，1)‖1≥‖MN(:，j)‖1(j>1),‖MN‖1=maxj‖MN(:，j)‖1=  log2N +1。
定理5.1得证。
定理5.2构造矩阵B的第p行的迭代次数为将p表示为二进制数(p)2时(p)2中包含的1的个数。
证明: 根据算法5.3的步骤6有操作pt=pt－lowbit(pt)，该操作的结果是将(pt)2中值为1的最低位置为0。而由步骤3,pt由p进行初始化。说明(p)2中包含多少个1，迭代次数就进行了多少次。
定理5.2得证。
由式(5.2)推出该策略矩阵所产生的均方误差为2ε2trace(BTB)Δ2MN=2ε2trace(BTB)( log2N +1)2同时，结合定理5.2可知，B的每一行至多有O(log2N)个元素为1，其余均为0。因此，B中有O(N log2N)个元素为1。即trace(BTB)的数值复杂度也为O(N log2 N)。又由定理5.1可知，MN的列范数 log2N +1的复杂度为O(log2N)。因而，根据式(5.2)可以求得总体的均方误差为O(N log32 N),而每条查询的均方误差则为O(log32N)。
综上所述，由树状数组构造的策略矩阵满足低敏感性以及均方误差复杂度低的特点，初步具备了在差分隐私下较为精确地进行连续数据发布的能力。然而，仅仅利用树状数组构造的策略矩阵并不能达到本章的精确性要求。接下来，将在此基础上寻找更精确的数据发布算法。
一般而言，可以用发布数据的一致性调节［7］、策略矩阵的权重系数调节等方法提高数据发布的精确性。然而，研究表明该问题已经满足线性一致性，具体分析如下。
定义5.3（线性一致性）对于负载矩阵W，记未加噪的查询结果为Y=WQ(D)，通过低秩矩阵机制获得的查询结果为Y′=A(W,D)。A(W,D)满足线性一致性当且仅当对任意可以表示成z=vY的行向量v都有vY′为定值，其中z为可以用Y线性表示的统计量。
定理5.3当式(5.1)中的矩阵M为行满秩矩阵时，低秩矩阵机制满足线性一致性。
证明: 当矩阵L为行满秩矩阵时，求得其右逆矩阵M+=MT(MMT)－1，满足MM+=I。
由于W=BM，因此B=WM+。
任取统计查询z，有多个vi满足z=viWX（其中vi之间互不相等）。
将其代入式(5.1)中，可得统计后的结果，令z′i表示由vi求得的查询结果:z′i=viBMX+ΔMεn=viWM+MX+ΔMεn=zM+MX+ΔMεn经过化简可以看出，z′i是与vi无关的噪声统计量。因此，z′i的值是相等的。
定理5.3得证。
而矩阵MN可知为可逆矩阵。结合定理5.3，可得出推论: 由树状数组构造基于矩阵机制的数据连续发布方法满足线性一致性，因此无法从线性一致性的角度提高数据发布的一致性。本章将在5.4.2节对策略矩阵的权重系数调节问题作进一步研究。
5.4.2查询均方误差的降低
5.4.1节分析了差分隐私下的连续数据发布的性质，并通过树状数组构造出基于矩阵机制的策略矩阵。本节在5.3节所描述的步骤的基础上进行优化。进一步研究矩阵MN，可发现该矩阵未饱和。以M3为例，表示如下: M3=100
110
001M′3=100
110
002计算可得‖M3‖1=2。若将M3的第三行乘以2，得到M′3，依旧满足‖M′3‖1=2，并不会影响整体的敏感度。同时，矩阵B也应转换为B′。B=100
010
011B′=100
010
010.5根据式(5.2)，可以直接求出转换前后两者之间的均方误差。转换前为32ε2，转换后为26ε2。经过转换，均方误差降低了,这说明直接由树状数组构造的矩阵MN优化得还不够彻底。经研究发现，可以通过在MN前面乘一个对角阵的方式提高精确性。
令ΣN=λ1
λ2
λ3表示N×N的系数对角阵,则可将式(5.7)拓展如下: A(W,D)=BΣ－1NΣNMNX+ΔΣNMNεn(5.8)式(5.8)即为添加系数对角阵后的隐私保护机制。当ΣN=IN时，该公式与式(5.7)等价。对应的均方误差公式如下: 2ε2trace(BTBΣ－2N)ΔΣNMN2(5.9)根据文献［17］的结论，令B′=αBΣ－1N，L′=α－1ΣNΜN，则有2ε2trace(B′TB′)Δ2L′=2ε2trace(BTBΣ－2N)ΔΣNMN2。因此，可将ΔΣNMN限制为|ΣNMN|1≤1，最小化trace(BTBΣ-2N)，则该优化问题可表示为如下形式: minΣNf(ΣN)=2ε2trace(BTBΣ－2N)
s.t.|ΣNMN|1≤1当上式取得最优解时，即等价于式(5.10)取得最优解。为简化推理过程，在式(5.10)中忽略了常数2ε2，实际计算时加上该常数即可。minΣNf(ΣN)=∑Ni=1B(:，i)TB(:，i)λ2i
s.t.MTNλ1
︙
λN≤IN×1
λi>0(5.10)其中B(:，i)表示矩阵B的第i列。
由分析可知，式(5.10)是一个线性约束下的凸优化问题。针对该问题，可直接采用SQP方法［20］求得最优解。然而SQP方法是一种时间复杂度很高的算法，无法满足大规模数据的要求。实验表明，对于一般计算机，该方法最多只能满足N<1000的求解规模。因此，需要进一步研究更快速的方法求得式(5.10)的最优解。
5.4.3最小误差的快速求解
由于SQP方法的时间复杂度很高，对于大规模数据，对角阵ΣN求解是无法完成的。因此，有必要对ΣN的求解进一步优化，并提出高效的解决方案。利用MN和B之间的特殊性质，本节提出一种高效的求解最小误差的算法——快速对角阵优化算法(Fast Diagonal Matrix Optimization Algorithm, FDA)。当N=2m－1时，该算法可以在O(log N)的时间复杂度下求解ΣN的任意系数值λi。该算法与未使用Σ前的方法有相当的求解效率，因此它保证了在不影响算法时间复杂度的前提下提高了隐私数据发布的精确性。该算法是基于以下定理提出的。
定理5.4令ΣN表示ΣN最优系数矩阵，则存在待定系数α使得Σ2m-1和Σ2m－1－1之间满足以下递推关系: Σ2m－1=αΣ2m－1
1－α
Σ2m－1(5.11)证明: 对于矩阵B2m－1进一步分析,可发现它满足以下特性: 令a<2m－1，b=2m－1+a。将其写成二进制形式可以描述为: a=(xm－2xm－3…x0)2和b=(1xm－2xm－3…x0)2。根据算法5.2，不难发现B2m－1(a,t)=1(t>0)当且仅当B2m－1(b,2m－1+t)=1。同时，由于b的2m-1位为1，因此B2m－1(b,2m－1)=1。
通过以上分析，可将B2m-1写成如下形式: B2m－1=B2m－1－1O（2m－1－1）×1O(2m－1－1)×(2m－1－1)
O1×(2m－1－1)1O1×(2m－1－1)
O(2m－1－1)×(2m－1－1)I(2m－1－1)×1B2m－1－1(5.12)通过式(5.12)得B2m－1(:，b)TB2m－1(:，b)=O(2m－1)×1
B2m－1－1(:，b－2m－1)TO(2m－1)×1
B2m－1－1(:，b－2m－1)
=B2m－1－1(:，a)TB2m－1－1(:，a)下面分析M2m-1与M2m-1-1之间的关系。
根据算法5.4，有M2m－1(t,a)=1(1≤t≤2m－1－1)当且仅当M2m-1(2m-1+t,b)=1，满足1≤t≤2m-1-1，M2m-1(2m-1,t)=1。
因此，将M2m-1与M2m－1－1写成如下递推关系: M2m－1=M2m－1－1O(2m－1－1)×1O(2m－1－1)×(2m－1－1)
I1×(2m－1－1)1O1×(2m－1－1)
O(2m－1－1)×(2m－1－1)O(2m－1－1)×1M2m－1－1(5.13)令RN表示ΣN的对角线元素组成的列向量， RN=(λ1λ2…λN)T。当N=2m-1时，可将R2m-1拆分成3个部分:R2m－1=(R(1)T2m－1λ2m－1R(2)T2m－1)T(5.14)其中,R(1)2m－1=(λ1λ2…λ2m－1－1)T，R(2)2m-1=(λ2m-1+1λ2m-1+2…λ2m-1)T。
对于式(5.10)，也可将其拆分为以下3个子部分: f(R2m－1)=∑2m－1－1i=1BT2m－1(:，i)B2m－1(:，i)λ2i①
+BT2m－1(:，2m－1)B2m－1(:，2m－1)λ2m－12②
+∑2m－1－1i=1BT2m－1(:，i)B2m－1(:，i)λ2m－1+i2③令f(i)()分别表示这3个子部分,从而将上述表达式转换为3个子部分的和: f(R2m－1)=f(1)(R(1)2m-1)+f(2)(λ2m-1)+f(3)(R(2)2m-1)而对于其限制条件，有MT2m-1R2m-1≤I（2m-1）×1。依照式(5.13)和式（5.14)展开得MT2m－1－1R(1)2m－1+λ2m－1I（2m－1－1）×1
λ2m－1
MT2m－1－1R(2)2m-1≤I(2m－1)×1(5.15)由式(5.15)可将限制条件分解成以下3个子条件: 
MT2m－1－1R(1)2m－1≤(1－λ2m－1)I（2m－1）×1
λ2m－1≤1
MT2m－1－1R(2)2m-1≤I(2m－1－1)×1①
②
③
通过以上3个子限制条件，可知子条件①受限于子条件②中的λ2m-1的取值。因此，先假设λ2m－1为待定系数，令λ2m－1=1－α(0<α<1)。
式(5.10)取最优时，子部分①满足: f(R2m－1)=minΣ2m－1f(R2m－1)f(1)(R(1)2m－1)=minΣ2m－1f(1)(R(1)2m－1)
MT2m-1R2m－1≤I(2m－1)×1MT2m－1－1R(1)2m－1≤(1－λ2m－1)I(2m－1)×1由于MT2m－1－1R(1)2m－1≤αI(2m－1)×1MT2m－1－11αR(1)2m－1≤I(2m－1)×1因此，令μi=1αλi,QN=1αRN =(μ1μ2…μN)，并将其代入式(5.10)的子部分①后有f(1)(R(1)2m－1)=∑2m－1－1i=1BT2m－1(:，i)B2m－1(:，i)λi2=1α2∑2m－1－1i=1BT2m－1(:，i)B2m－1(:，i)(μi)2
=1α2f(1)(Q(1)2m－1)通过以上分析，可将式（5.10）的子部分①的问题描述如下: opt: minQ(1)2m－11α2f(1)(Q(1)2m－1)opt: minQ(1)2m－1f(1)(Q(1)2m－1)
s.t.MT2m－1－1(Q(1)2m－1)≤I(2m－1)×1
μi>0将Q(1)2m－1用R2m－1－1代入，则问题等价于求解R2m－1－1，即Σ2m－1－1。由此可得Q(1)2m－1=R2m－1－1 =1αR(1)2m－1，即R(1)2m－1=αR2m－1－1。
而式（5.10）的子部分③可看成是α＝1的特殊情况。因此，可以得出R(2)2m－1=R2m－1－1。
综上所述，式(5.11)成立。
定理5.4得证。
由定理5.4，可得形如式(5.10)的ΣN递推关系，从而可由Σ2m－1－1的结果来求解关于Σ2m－1的最优结果。
假设已经求得N=2m－1－1下的最小均方误差errm－1=minΣ2m－1－1f(Σ2m－1－1)及Σ2m－1－1。将上述问题转化为关于α的最优化问题: opt:h(α)=minαerrm－1α2+2m－1(1－α)2+errm－1