信息熵应用随笔3：多维联合直方图的存储

站在计算机科学的角度，多维联合直方图的一个重要问题是，如何表示和存储其计算结果。

一，多变量联合直方图的多维数组存储法

一个简单的方法是将多变量联合直方图直接存储为多维数组。

• 一维变量V1的直方图可以存储成一维数组bin[v₁]=f(v₁), 其中，v₁是变量V₁上的bin下标，bin[v₁]用数组的方式保存了变量V₁上的bin位置信息，f(v₁)则是bin对应的频率。
• 二维变量V₁，V₂ 的联合直方图可以存储为二维数组bin[v₁][v₂] = f(v₁,v₂)。其中，v₁是变量V₁上的bin下标，v₂是变量V₂上的bin下标，
• n维变量V₁，V₂， ... , V_n的联合直方图可以存储为n维数组bin[v₁][v₂]...[v₂]=f(v₁, v₂, ..., v_n)。其中，v_i是变量V_i上的bin下标。

假设有N个数据维度，我们对所有维度都使用B个bin，如果频率以浮点数进行存储，那么对于一个多维联合直方图，若使用多维数组存储，我们就需要保存B^N个浮点数。显然，使用多维数组存储多维联合直方图，其内存需求会随着数据维度和bin数量的增加而指数型增加。这对存储空间的消耗也指数型增长的。这给我们处理多维联合直方图造成了很大的麻烦。

二，超稀疏多维数组压缩存储法

1. 多维直方图的压缩原理

为了减少存储多维直方图的内存需求，我们使用了多维直方图的稀疏性。给定一个具有P数据点的数据集，甭管这个数据集有多少个数据维度，画成的多变量直方图也最多只有P项具有非零频率，而所有其他项都是零。对于大规模数据集，如果我们将整个数据集划分为更小的块，并计算每个块的多变量直方图，则每个数据块的数据点数量甚至小于多变量直方图中的数据项数量。假设一个数据块的大小为8的3次方，每个数据点有3个属性（即数据维度），如果我们构造一个多变量直方图，每个维度有256个bin，那么至多有8^3/256^3≈0.0031%的多变量直方图数据项不是零。考虑到多维联合直方图的多维数组实际是超稀疏的，故可以采用数据空间变换来减小多元直方图的大小。

2. 超稀疏多维数组的空间变换方法

设V=V₀，V₁，…，V_n-1为一个多元直方图的变量集，记Bi为第i个变量的bin个数。那么，多维直方图可以表示为一个N维数组M，其大小为B0 * B1 * ... * Bn-1。为了减小多元直方图的大小，我们的目标是：用数据空间转换操作，将M0转换成另一个N维数组Mt, 其大小为P0 * P1 * ... * Pn-1，这里要求Pi << Bi 。

为了做到以上变换，需要构造一个字典来记录Mt的bin索引到M0的bin索引之间的1对1映射。让Di表示第i个变量构造的字典，给定bin索引j∈[0，Pi−1]，Di可以将j映射到一个唯一值k∈[0，B_i-1]，k即原始bin索引，我们将此步骤称为解码。

为了生成字典，对于每个变量Vi，让bin[v_i]表示其bin索引，bin[v_i]∈[0，Bi] (变量Vi的bin索引上下限是0到Bi)。而直方图数据项以原始的多维数组表示时，应为bin[v₁][v₂]...[v₂]=f(v₁, v₂, ..., v_n)。如果带有bin[v_i]索引的直方图数据项都等于0，则删除所有这些数据项。否则，我们向字典中插入一个新条目。当所有变量都经过空间转换后，将得到一个更小的多维数组，从而减少了存储开销。

3. 一个二维数组空间变换实例

图1，例子：一个二维直方图对应的二维数组

请看一个简单的例子来了解如何做变换的。有个二维直方图，包含变量V₀（对应图1的X轴）和V₁（对应图1的Y轴），每个变量各有15个bin。可以看到这个二维直方图保存为二维数组后，大多数数据项（方格）是空的（频率等于0）。为了降低数据存储开销，我们先计算V0变量的bin索引bin[v₀=0]对应的二维直方图数据项之和A(v₀=0)。参与这个计算的二维直方图数据项如上图蓝框所示。用公式表达即A(v₀=0) = ∑f(v₀=0，v1=i), i∈[0，14]。如果A(v₀=0) =0，则从二维数组中删除所有这些数据项。否则，我们向字典中插入一个新条目。依次类推，算完V₀变量的全部bin索引后，再算V₁变量的全部bin索引。最终，转换完成的二维数组如下图所示，只有5*5的大小，是原始二维数组的1/9。

图2，压缩后的二维数组

4. 字典表的构造

在数据空间转换之后，我们将压缩后的多维数组表示为索引和频率对的集合，其中每个键值对对应于多维数组中的一个非空条目。图2的例子就可以表示成：

• bin[v₀=8, v₁ =1] : 7
• bin[v₀=5, v₁ =4] : 6
• bin[v₀=9, v₁=4] : 3
• bin[v₀=7,v₁=6] : 6
• ......

可以看出此集合的键值对总数等于多维数组中非0项的个数。此集合可以无损地转换为原始的多维数组。这个集合就是我们的数据转换字典表。求字典表的计算过程非常简单，只要将原始的多维数组中非零项取出，组成如图3下所示的键值对即可。

图3，字典表中索引的表示方法

在字典表中，可以把指示多个变量的索引的数值放在一起，只要我们知道那几位表示的是哪个变量就行。还是之前的二维联合直方图的例子，其字典表如图3上所示。我们可以将索引值转换为二进制，对应于第一个变量的位被涂成红色，对应于第二个变量的位被涂成绿色，如图3下半部分所示。这样看起来更明显，也更方便计算：我们可以操纵这些位来操作不同变量的bin索引，从而计算边缘分布，条件分布。

 三、基于多维直方图压缩存储法的查询操作

1. 多维联合直方图的边缘化运算（Histogram Marginalization）

边缘分布（Marginal Distribution）指在概率论的多维随机变量中，只包含其中部分变量的概率分布，本质一种降维操作。考虑到多维联合直方图相当于多个变量的联合分布矩阵，多维联合直方图的边缘化运算就是计算其部分变量的边缘分布，也就是包含其部分变量的子集。为了构造这个新的直方图，我们将预先存储的多变量柱状图边缘化到用户没有选择的变量上。

假设原直方图涉及四个变量A、B、C和D，用户对B和D的柱状图感兴趣，想计算只包含B和D的子直方图。为此，我们将变量A和C上的多变量直方图边缘化，即：

P(B=b，D=d) = ∑_a∑_c P(A=a，B=b，C=c，D=d)。

基于上文描述的多维直方图压缩存储法，字典表索引中的不同位对应于不同变量的bin索引，因此可以使用按位与运算，然后使用频率和，轻松地执行柱状图边缘化操作。不需要对索引进行解码将转换后的bin索引映射到原始bin索引。

图4 求只包含变量B和D的子直方图

具体操作过程：包含四个变量A，B，C，D的多维直方图的字典表如图4(1)所示。每个变量对应两位的二进制索引。我们要求只包含B和D的子直方图。可以通过将用户未选择的变量对应的位设置为0，然后求具有相同索引的频率之和来完成边缘化。这一步可以用计算机的方法来做：

• 设置一个掩码，要取的变量bit设为1，未选的变量bit设为0。在这个例子里ABCD四个变量取B和D，则掩码为00110011，如图4(2)所示；
• 然后该掩码与各个变量的索引bit位按位与，结果如图4(3)所示；
• 将图4(3)的结果，将其中具有相同索引bit位对应的频率加起来，得到最终的边缘分布结果，如图4(4)所示。

2. 多维联合直方图的bin合并运算（Histogram Bin Merge）

多维联合直方图的bin合并，是相邻仓位中的频率沿某些维度进行组合，得到一个新的具有不同离散化水平的直方图。当某些分析任务需要不同分辨率的直方图时（分辨率可以直观地理解为bin的尺寸），使用此操作。

我们假设预先计算的多变量直方图是以最高分辨率（最小bin的尺寸）定义的。为了方便起见，在数据空间转换之前，每个变量的bin个数设置为2的幂。假设具有4个变量的高分辨率多维联合直方图在数据空间转换前的bin数为2^K1、2^K2、2^K3、2^K4，其中K1=K2=K3=K4=4，也就是bin数均为16。要求对这四个变量按照每组包含2¹、2²、2¹、2⁰个bin进行bin合并。最终计算出新的低分辨率直方图，每个变量应具有2^K1-1=8、2^K2-2=4、2^K3-1=8、2^Kn-0=16个bin。

图5，多维联合直方图的bin合并运算例子

对于这个要求，用户需要提供一个4维向量 L={1,2,1,0}, 代表每个维度下bin合并的指标。具体的bin合并计算，需要先解码bin索引，再进行操作。对于bin尺寸不变化的维度，也就是Li=0的维度，则不需要解码bin索引。在本例中只有A，B，C三个尺度需要解码，D不需要。例子直方图的字典表如图5(1)所示。解码bin索引后的结果如图5(2)所示。接下来的操作与多维联合直方图的边缘化运算一样，只是使用的位掩码不同。位掩码的计算方式如下图所示：

图6，掩码的计算方式

将掩码与解码后的bin索引按位与，结果如图5(4)所示。接下来将具有相同索引的频率加起来，得到最终结果，如图5(5)所示。

3. 条件直方图的运算（Conditional Histogram）

计算条件直方图就是计算条件分布，是计算条件熵的关键。

给定n个变量的变量集V=V₀，V₁，…，V_n-1，V有一个子集U ⊂ V， 和另一个子集W∈(V-U)。条件直方图的运算，是计算W中的变量在特定的bin范围下在U中的分布。当用户要研究一组变量U和另一组变量W之间的关系式，条件直方图是有用的。例如，我们可以计算P(U|W)的熵，它表示在特定值的范围内，给定W中的变量，U的不确定性。

计算条件柱状图的第一步是为w中的变量选择属于特定bin范围的所有条目。对于此选择步骤，我们需要对w中的变量进行解码，以获得其原始bin索引。然后，对于每个选定的条目，使用位掩码执行逻辑与操作，以将与U中的变量不对应的所有位设置为0。通过对具有相同索引的条目的频率求和来计算最终结果。

也是让我们看一个具体的例子：有一个具体的4维直方图，包含A，B，C，D，四个变量，每个变量有16个bin。其字典表如图7(1)所示。假设U包含两个变量B和C，W包含两个变量A和D。A的数值范围是(min_A, max_A)，D的数值范围是(min_D, max_D)，这条件直方图的查询操作定义为：

P(bin(B)=b, bin(C)=c |  min_A  ≤ bin(A) ≤ max_A , min_D  ≤ bin(D) ≤ max_D ) =  ∑_a=minA∑_d=minD P(bin(A)=a, bin(B)=b, bin(C)=c, bin(D)=d)

图7，条件直方图的运算示例

现在要求查询 P(A|11 ≤ bin(C) ≤ 15) 。我们首先将C的bin索引解码，获得原始的bin索引，然后选定11 ≤ bin(C) ≤ 15范围内的频率，如图7(2)绿色部分所示。接下来设置掩码，掩码只有要查询的A变量索引位是1，其他不查询的变量对应的索引位都是0，如图7(3)所示。接下来将图7(2)选定的绿色项与掩码按位与，并将相同bin索引的频率求和。结果如图7(5)所示。在数据空间转换完毕后，每个变量有4个bin。

时序直方图(Time Histogram)

对于时变多变量数据，我们将时间步数（time step）视为另一个变量，并将bin数设置为时间步数。查询的时序直方图有助于显示此选定数据块中随时间显示的功能。时序直方图查询与直方图边缘化查询类似，只是可变时间步长必须包含在边缘变量中。在查询过程中，用户指定用于计算时间柱状图的变量和数据块。变量时间步必须包含在这些选定变量中，然后查询的其余部分与柱状图边缘化查询相同。

 四、总结

以上介绍了一种多变量联合直方图的多维数组压缩存储法，以及使用这种存储方法如何进行的查询操作：包括多维联合直方图的边缘化运算，多维联合直方图的bin合并运算，条件直方图的运算，和时序直方图。可以看到文中所示的压缩存储法，不仅可以减小空间消耗，同时也支持多种运算。

回顾一下我们的分析路径：

• 如何在一个多维数据集中找到最具代表性的数据维度？
• 可以用信息论方法，选择具有最大信息增益的数据维度。
• 如何计算信息增益？需要计算联合熵和条件熵。
• 如何计算联合熵和条件熵？需要计算联合概率分布。
• 如何计算联合概率分布？可以用直方图做拟合。
• 多变量联合直方图如何表示和计算？可以用多维数组。
• 多维数组表示多变量直方图，其空间消耗呈指数型上升，不便于存储和计算，怎么办？可以用压缩的二维数组存储法。
• 有了压缩的二维数组存储法，如何计算联合熵和条件熵？下一节我们就具体讲这个。

参考文献：

Lu, K. and H.W. Shen. A compact multivariate histogram representation for query-driven visualization. 2015.