增量式分级判别回归树（IHDR）|翻译与笔记

为什么翻译这篇发表于2007年的老论文呢？主要有以下四方面原因： （1）这篇论文较清晰地给出增量式学习的概念。 （2）论文有一定的理论与应用方面的价值。 （3）论文的作者仍处于活跃状态，阅读前面的工作有助于理解最新的工作。 （4）经过一轮改换，感觉亲近些，读自己翻译的文本更容易理解些。 以下是译文，笔记会以斜体字显示。

---

[1] J. Weng, “Incremental Hierarchical Discriminant Regression,” p. 28, 2007.

摘要

该论文提出了Incremental Hierarchical Discriminant Regression (IHDR)。IHDR能够在线地给非常高维的回归或者决策空间增量式的构建一棵决策树或者回归树。其中的生物启发是大脑联想皮质的自动发育计算模型。这个联想模型是关于从上往下的感知输入到从下往上的运动投射。在IHDR树的每一个内部节点中，输出的信息直接被用来自动的导出局部子空间（这些子空间是由那些能按照输出的聚类对输入空间进行划分，使得各子空间区分度最高的特征）。方法中的一些参数都是自动的以及数据驱动的，这使得IHDR树能够自动的拟合未知的分布形状。一个嵌入到树中的分层概率分布模型，它被用来对那些在搜索中最不可能的情况进行剪枝。IHDR树动态地分配长期记忆来避免遗忘性问题（普遍存在于用于全局拟合学习的算法，例如神经网络）。在构建树的过程中，可能有大量的样本相差性非常大。对于这种情况，一个增量式更新的概率模型（样本大小相关的负对数似然SDNLL度量）被用来判断是大样本、小样本、不均衡的样本等情况。

keywords : online learning, incremental learning, cortical development, discrimenant analysis, local invariance, plasticity, decision trees, high dimensional data, classification, regression, and autonomus development.

引言

大脑皮质是执行数值输入到数据输出的回归模型。在一个合适的时间尺度上，神经元激活的频率已经被建模成一个由神经元轴突传递的主要的信息源。一块皮层区域会随着信号处理的经验逐步的发育（adapts and self-organizes）。人脑有四种皮质：初级的、联想（关联）的、多模态、以及旁系的。联想皮质是在初级皮质与运动区域之间。而本文构建的数学模型可以被认为是一种对大脑联想皮质的发育的一个粗糙的近似的计算模型。事实上，大脑对于人类来说，还有很多未知的地方。所以，本文的计算模型并不打算完全吻合生物学里的内容。

分类树（输出为类标签）与回归树（输出为数值向量）有两个目的：索引与预测。索引树，如K-D树和R树，已经被广泛地用于从数据集中索回已知树据，作为一种有效的数据存储结构，以降低搜索的时间复杂度（$\log{(n)}$）。预测树也被称为决策树被广泛地用于机器学习中，被用来生成一个基于树的预测规则集。因为构建的决策树通常是浅层的，所以时间复杂度对于决策树来说不是问题。并且，本文需要一个实时、在线、增量式学习的人工神经网络系统。所以，我们即考虑了时间效率与重要的泛化性能。而且，我们需要增量式的构建树，从而使得其能用于处理增量式到来的数据。例如实时输入到机器人系统的环境视频信息。

传统的分类或者回归树在每个内部节点利用单变量对其进行分裂的方式来构建树，如CART,C5.0，其他很多树方法。这意味着用于分割节点中输入空间的超平面与空间某一坐标轴正交。多变量线性分裂不需要与任何坐标轴正交，因此能够沿着决策边界对空间进行合适地划分，从而得到更好的预测与泛化效果。利用多变量进行构造的树被称为oblique trees。

已经存在的oblique trees有OC1,SHOSLIF。OC1利用递归地搜索一个分裂平面。SHOSLIF利用PCA与LDA直接计算分裂。SHOSLIF利用多变量非线性分裂（任意方向的曲面分割面）。LDA为什么重要，因为LDA利用了输出空间的信息来分裂输入空间。PCA仅仅用了输入空间的信息。这样做的结果是，输入空间中与输出无关的变量（例如噪声成分）也被PCA捕获。而包含判别分析技术，例如LDA，能够无视与输出无关的输入成分。

当输出是一个多维度的模拟信号时（例如实时的机器人运动控制信号），问题变得更难。类标签难以得到，因此LDA方法不能直接应用。Hierarchical Discriminant Regression (HDR) 在输入与输出空间中都执行聚类操作，输出空间的聚类提供虚拟的标签来引导输入空间的聚类。本文提出HDR的增量式版本，IHDR。

增量式回归

完全地增量式的学习是一个非常难的问题。完全的增量式学习意味着树必须对于每一个输入向量对树（模型）进行更新。把大脑假设成一个完全增量式的并不是不合理的，它必须对于每一个感知输入进行操作和更新。在增量式学习系统中，数据帧以序列的形式依次到来。每一帧$x(t)$数据向量都是变化世界在时间$t$的一个样本（快照）。因为，观测数据流可能延续的时间比较长，甚至是不会终止的，所以系统必须能够每次利用当前数据对系统进行更新一次。每一帧数据被用于更新之后，便被丢弃。序列性的学习的主要原因包括：(a)将所有到来的数据都存储是不可能的；(b)批处理（包括分块序列的学习）需要时间并引入时间延迟，因为从第一批数据到下一批数据到来之前的这中间的时间是可观的；(c)对树进行更新是非常快的，然而构造一棵完整的树的时间较长。因此，利用一帧数据对树进行更新所带来的延迟效应更短。

增量式学习性质的优势

(a)边构建边执行。IHDR能够在构建完成之前就能工作（用于预测）； (b)专注于难解决的学习部分。后来的训练样本是否被采集是依据当前树的性能，这使得更有效的训练样本能够被用来训练难以学习的部分； (c)训练样本的数量是动态确定的。给定一个分类，或者回归任务，我们很难确定该任务需要多少训练样本。而增量式学习使得机器人依据当前的系统性能来动态的确定所需的训练样本数据； (d)这些在不同时刻被接收到样本（存在于IHDR中的叶子结点中的micro-cluster），能够存储在树的不同位置，能够提升性能。

IHDR相比其它增量式学习方法，具有如下特性

• 具有系统化组织的长期记忆与短期记忆。长期记忆主要包括，浅层节点中的信息，叶子节点中的micro-clusters（不经常被访问）。短期记忆包括，叶子节点中的micro-clusters（经常被访问）。通过增量式平均来忽略细节。长期记忆使得IHDR避免了灾难性遗忘等问题（普遍存在于基于BP算的前向、反馈网络）。

• 参数是动态、自动被确定的，包括树的自由度。这些参数包括每个内部节点对应概率模型的均值与方差矩阵。树的自由度也不是固定的，意味着树的层数，以及每个节点具有的子节点个数都是根据学习的任务自动调整的。

• 由粗至精的多层分布估计。

• 快速的局部更新，不需要迭代。一些现有的方法，例如进化计算和模拟退火能够处理局部最小值问题，但是它们需要迭代计算。而这不适用于实时更新。

方法

A. 分类与回归问题的统一

根据输出的类型（类标或符号输出与数值输出），可以将判别分析任务分成两类：分类与回归。

分类任务可以定义如下： 给定训练样本集$L=\{(x_i, l_i)|i=1,2,..., n\}$,此处$x_i\in X$为输入（特征）向量，$l_i$为$x_i$对应的符号标签，分类任务是为了为任意未知类别的输入$x\in X$确定其所属的类别。

回归任务与分类任务相似，只是类标$l_i$替换成$y_i$，$y_i\in Y, i=1,2,...,n$。回归任务定义如下： 给定训练数据集$L'=\{(x_i, y_i)|i=1,2,...,n\}$，回归任务是为了对任意的输入向量$x\in X$估计其对应的输出向量$y(x)\in Y$。

回归任务非常常见。只要系统需要利用变化的值控制执行器，学习问题就是一个回归任务。例如：电机、航向、刹车、以及各种工业设置的机器。生物学中的大脑皮质也执行回归任务。

在分类问题中，类标并不能提供更多的信息来表明各类别之间的差异性，从而无法评价，某一类与其中一类比另一类更相似或更有区别。任意不同的两类都是完全不同的，不是一个平滑的关系，在类的边缘是突变的关系。而实际上，某一类别的确与其中一些类别更相近，而与另一些类更远。例如，把老虎判断为猫比判断为汽车更加让人容易接受。

我们可以将分类任务投射成一个回归任务，从而能够通过融合一些原始的类来方便的形成粗糙的类。这些粗糙的类被用来利用决策树实现由粗到精的分类与回归。本文中在后文就采取了这种做法。

此处，我们列出三种将分类任务投射成回归任务的方法。

1. 标准映射（Canonical mapping）。将$n$个类标映射到$n$维的输出空间。对于第$i$个类，对应的输出向量$y_i$为$n$维向量，其中第$i$个元素为1,其它所有的元素为0。该方法在增量式学习背景下，具有很大的限制，因为最大的类别数被限定为$n$。
2. 嵌入的代价矩阵（Embedding cost matrix）。如果代价矩阵$[c_{i,j}]$可用，则$n$个类标可以被嵌入进一个$(n-1)$维的输出空间中，利用向量$y_i$表示类别$i$，使得$||y_i-y_j||$与$c_{i,j}$接近。这个过程并不是一直可行，因为预定义的代价矩阵$[c_{i,j}]$并不经常容易得到。这里的类别数也是限定的，并不适用于增量式学习。
3. 类别的输入空间的均值。每个属于类别$l_i$的样本$(x_{i,j},l_i)$，都转换成$(x_{i,j},y_j)$，此处$y_i$为所有属于同一类别$l_i$的$x_{i,j}$的均值。在增量式学习中，该均值可以由amnesic average更新。该方法是比较优选的方案，因为在大部分学习任务中都满足，输入相似意味着相似的输出，这一假设。

B. 动机

在本文剩余部分，我们考虑一个普通的回归问题：增量式的近似一个映射:$h:X\rightarrow Y$

假设训练样本为$\{(x_i, y_i)|x_i\in X, y_i \in Y, i=1,2,...,n\}$。增量式近似的意义为，发展器d只利用前一刻的映射模型$h^{(i-1)}$以及当前的输入样本$(x_i,y_i)$来生成新的映射$h^{(i)}$：$h^{(i)}=d(h^{(i-1)}, x_i, y_i)$

IHDR树中总共有三种clusters（也可称为neurons）：x-clusters,y-clusters,and micro-clusters。

• x-clusters与y-clusters：为了刻画每个内部节点的判别特征子空间D的划分边界。y-clusters是在$Y$空间中表示，用来给x-clusters提供虚拟的标签。x-clusters是定义在$X$空间中，为了对沿着内部节点上的判别子空间D对输入空间进行划分。
• micro-cluster：是由$(x_i,y_i)$表示的，只存在于叶子节点上。当然，如果IHDR树允许删减操作，那么内部节点则需要micro-cluster，因为当某个叶节点被删除时，它的父节点能够直接充当叶节点。

D. 双聚类

IHDR的每一个内部节点都增量式的维护y-clusters和x-clusters。如下图所示，对于每个节点的每一个类别最多能分成q类(q是一个设定的最大值). 数学上的描述为，输入空间$X$的聚类是以输出空间$Y$作为条件的。所以，对于每一个x-cluster，都是随机变量$x\in X$关于随机变量$y\in Y$的条件概率分布，记为$p(x|Y\in c_i)$，其中$c_i$是第$i$个y-cluster，$i-=1,2,...,q$。

$q$个y-clusters基于$y$确定了每一个到来样本$(x,y)$的虚拟标签。虚拟标签被用来确定对于当前的$(x,y)$，应该利用它的$x$对哪一个x-cluster进行更新。每一个x-cluster近似在$X$空间中与之相关的样本总体。如果需要一个更加精细的近似，它可能会产生一个子节点。

增量式更新的步骤如下：

在每一个节点，当前样本$(x,y)$中的$y$，被用来找到离它最近的y-cluster（欧氏距离）并用它来更新这个y-cluster的中心。该最近的y-cluster指出了与$(x,y)$相关的x-cluster。然后$(x,y)$被用来更新x-cluster的统计学信息（均值向量与方差矩阵）。每一个x-cluster的统计学信息被用来估计当前样本$(x,y)$属于各x-cluster的概率（此处建模成多维高斯分布）。

此外，x-clusters的中心提供了许多重要信息，这些信息能够用来分裂子空间。因为这些x-clusters是由根据$Y$空间中的虚拟产生的。我们定义最大判别特征（MDF）空间$D$作为分开这些x-clusters的线性空间。q个x-clusters需要q-1个判别特征（撑起q-1维的判别空间D）

为什么称为最大判别空间？例如，假设我们的输入有两个特征，其中一个包含与输出相关的信息，另一个与输出无关（如下图所示）。因为两个x-clusters的中心，有一个特征的值基本一样，所以，该特征被忽略。另一个具有很大差异的特征对于分开这两个类很重要。从而达到将与输出无关的“噪声”特征屏蔽掉，减小对学习性能的影响。如下图所示的，最大判别特征，不一定是与某个特征所在的轴正交，它的方向是任意的。

E. 自适应局部准(quasi)不变性

输入样本$(x,y)$是被依序的接收。连续的向量$x(t)$与$x(t+1)$或许与一个跟踪对象的图片有关（例如，离开、沿着一个方向移动、旋转、变形、面部表情、光线变化等）.这种准不变性是自适应的，并且是局部的。

自适应意味着准不变性是从感知动作经验中推导出来的，不是手工设计的。局部性意味着准不变性能被应用于输入空间中的局部流形，而不是完整的输入空间。而准(quasi)意味着这个不变性是近似的，并不是完全无误差的。

到目前为止，还没有类似的方法能够有效的在增量式学习的前提下做到很好的局部准不变性。一般利用欧氏距离来测量两个向量的差异。

利用以上方法能够获取自适应准不变的最大判别特征子空间。

如上图的（a）所示，输入对应四种输出，用不同的符号表示。但是如果不考虑输出的信息，那么这些状态就程现出一团糟。看不到明显的规律。

而图（b）中，因为考虑了输出的信息，所以能够更合理的将这些样本分开。

读者补充：所有的这些都是建立在，相似的输入，具有相似的输出这一假设之上。要谨记！！！

F. 生物学观点

如下图所示为一个IHDR树。每一个节点都类似一个皮质区域。状态空间$X$是由d维的向量表示，它最初被粗糙的分成4份。每个部分都是根节点的子节点。每个子节点都具有它自己的小区域。这样从粗到精的划分，将输入空间逐渐划分成很小的区域。如果一个子结点接收到足够多的样本，它将又划分成q个子节点。在叶子结点中，micro-clusters以$(x_i, y_i)$的形式存在。

如果$x$给定，但是$y$未知，我们可以通过搜索IHDR树直到叶子节点，从而获取对应输出值作为当前状态的预测。

IHDR算法

下面介绍增量式构建IHDR树的算法。算法的训练数据是以按时间依次输入给机器人，将$t$时刻接收到的样本标记为$(x_t, y_t),t=1,2,...$。构建的树越深（层数越多），则x-clusters的方差越小。当节点中的样本数量太小，以致于不能得到一个对x-clusters的统计估计，则这个节点仍作为叶子叶点。

Procedure 1：IHDR Algorithm。

• 初始化根结点。
• For t=0,1,2,...do the following: 获取当前的样本$(x_t,y_t)$。 如果$y_t$已知，则调用update-tree(root, $x_t$, $y_t$)，否则，调用compute-response(root,$x_t$,$y$)，并且输出$y$。

A. Update tree

下面的是增量式更新树的算法。给定树的根结点与训练样本$(x,y)$，仅仅利用当前的样本$(x,y)$对树进行更新。每次调用update-tree函数，都可能对树进行增长。程序处理自适应的皮质可塑性（三层）。系统参数：q是每个内部节点的最大的子聚类个数。$b_s$是每个常量参数（自由度）所需要的样本数（本文取$b_s=20$）。

Procedure 2: Update-tree(root,$x$, $y$)

• 1）初始化。令A为被激活的节点（等待搜索），A被初始化为根节点。
• 2. 树搜索。当A不是叶子节点时，执行如下： a) 计算响应。对于输入$x$与A的每一个神经元，计算它们的响应。 b)竞争。对于A的所有q个子节点，拥有最大响应的子节点胜出。 c)A中胜出的子节点，被当作下一个激活的节点。
• 3. 更新可塑的内部节点。更新激活叶子节点的父节点A（注意：to avoid large-scale temporal space partition for constrained plasticity，叶子节点父节点之上的祖先节点不用更新）调用update-note(A,$x$,$y$)更新节点A。
• 4)增量式聚类。调用update-cluster-pair(A,$x$,$y$)。
• 5)对于叶子节点A，如果样本的数量$n$大于一个阈值（自动计算出来的$NSPP=2(n-q)/q^2>b_s$），则将A作为一个内部节点，并且生成q个叶子结点（通过将A的micro-clusters分配到他的新的子结点中）。

B. Update node

给定一个节点N以及一个训练样本$(x,y)$，它利用$(x,y)$更新节点N。这个程序处理皮质块（在树叶点内部）的自适应可塑性。

Procedure 3: Update-node(N, x,$$y)

1. 更新y-clusters。令c为N节点中的y-clusters集合。通过调用update-clusters($y$,$c$),该函数返回最近的$y_i$的索引。 2）更新与上一步y-cluster索引相关联的x-cluster。利用$x$更新x-cluster的均值，其中引入遗忘机制来加大新数据的作用。 3）利用最大的判别特征来更新节点N的子空间，因为第i个x-cluster已经被更新了。所以要对N节点的子空间进行重新划分。

C. Update clusters

给定样本$y$以及聚类集合$c=\{y_i|i=1,2,...,n\}$，程序更新$c$中的聚类。参数包括:q为$c$中聚类的最大个数；$\delta_y>0$为输出空间的分辨率。

Procedure 4: Update-clusters($y$, $c$)

1）如果$n<q$并且$||y-y_j||>\delta_y$（为了防止非常小的或者一样的样本产生多个聚类），增加$n=n+1$，设置新的聚类$y_n=y$，并将$y_n$添加到$c$中，返回。否则，执行以下操作： 2）按照以下度量找到最近邻的$y_j$：$j=arg\min_{1\leq i \leq n}\{||y_i-y||\}$。 3)更新最近聚类的确定的部分$p$（例如，$p=0.2$只更新聚类中20%），引入遗忘机制，并返回$j$。

D. Update-cluster-pair

给定样本$(x,y)$以及聚类对的集合$c=\{(x_i,y_i)|i=1,2,...,n\}$。该程序更新$c$中最匹配的聚类。该程序是对每一个输出神经元（micro-cluster）的局部操作。参数包括：$b_l>0$为叶子结点中micro-clusters的最大允许个数；$\delta_x>0$为输入空间的分辨率。

Procedure 5: Update-cluster-pair($x$,$y$,$c$)

1)如果$n<b_l$并且$||x-x_j||>\delta_X$，$n=n+1$，并设定一个新的聚类$(x_n,y_n)=(x,y)$，并将$(x_n,y_n)$增加进$c$，并返回。否则，执行以下操作： 2)利用如下度量找到最近的$x_j$：$j=arg\min_{1 \leq i \leq n}\{||x_i-x||\}$。 3）利用新样本$x$更新聚类$x_j$，并引入遗忘机制。 4）利用新样本$y$更新聚类$y_j$，并引入遗忘机制。 5）返回更新的$c$。

E. Compute response

当样本只有输入$x$，而不知道输出时。我们调用IHDR的compute-response程序来得到$x$的预测。

Procedure 6: Compute-response(root,$x$,$y$)

1. 先执行update-tree中同样的操作，找到最合适的叶子结点A。 2）计算输出$y$。令叶子结点中的micro-clusters为$c=\{(x_i,y_i)|i=1,2,...,n\}$。利用以下度量找到与$x$最合适的匹配：$j=arg\min_{1\leq i \leq n}\{||x-x_i||\}$令预测的输出$y=y_j$，返回。

F. 遗忘机制

遗忘机制的引入是由大脑神经元的可塑性启发的（需要对当前的经验作出自适应的调整）。它还受以下由数学语言描述的统计有效性：为了估计一个分布的均值向量$\theta$（例如观测值 $x_t,t=1,2,...,$的均值向量作为神经元突触的权值 向量），样本的均值是对一个大类时不变分布的最有效估计（例如，指数分布簇，像高斯分布）。通过定义最有效的估计$\Gamma$是所有可能估计中具有最少的期望误差，$E||\Gamma-\theta||^2$。然而，在实际情况中，观测分布并不是时不变的，此时就需要利用遗忘机制来得缓慢的改变分布，并同时保持估计的准最优有效性。

从增量式学习算法的观点看，初始的状态聚类的中心是由早期的数据决定的。当更多的数据可用时，这些中心移动到更加合适的位置。如果这些聚类的中心被用来判别每个聚类的边界，这些初始的数据将不能被正确的分类。换句话说，每个类别的中心包含的一些早期数据，可能此刻已经不属于这个类了。为了减少早期数据影响，遗忘机制能够被用来计算每一个聚类的中心。遗忘机制能更好的跟踪输入环境的动态性。

对于每个输入数据$x_1,x_2,...,x_t$的均值给出如下$\overline{x}^{(t)}=\frac{1}{t}\sum_{t=1}^t x_i=\sum_{i=1}^t \frac{1}{t}x_i$

对于上面的表述，每个$x_i$通过乘以一个权重$\frac{1}{t}$，然后加和。这被称为等权值平均。如果$x_i$增量式的到来，我们需要增量式的更新均值 ：$\overline{x}^{(t)}=\frac{(t-1)\overline{x}^{(t-1)}+x_t}{t}=\frac{t-1}{t}\overline{x}^{(t-1)}+\frac{1}{t}x_t$

而引入遗忘机制的平均如下：$\overline{x}^{(t)}=\frac{t-1-\mu}{t}\overline{x}^{(t-1)}+\frac{1+\mu}{t}x_t$。

当$\mu\geq0$，是一个遗忘参数。

但是随着数据不断接收，样本的个数不断增加，也即$t$不断增大。如果$\mu$是一个固定的值，那么在$t$很大时，它的作用几乎为零。因此，更加合适的做法是，$\mu$应该随着时间缓慢的增大，也即为时间的函数$\mu(t)$。可表示成一个分段函数，例如，当$t\leq t_1$时，不引入遗忘参数，即$\mu(t)=0$；当$t_1 < t \leq t_1$，我们使得$\mu$缓慢的从0变到$c$(例如，$c=1$)；当$t_2<t$时，我们令$\mu(t)$缓慢的变化，以增长率$1/m$增加。综上，可表示如下$\mu(t)=\left\{ \begin{matrix} 0 & t\leq t_1\\c(t-t_1)/(t_2-t_1) & t_1 < t \leq t_2\\ c+(t-t_1)/m & t_2< t\end{matrix}\right.$

因为$\lim_{t\rightarrow \infty}(1+\mu(t))/t=1/m$,当$t$增加到一个很大的值时，新数据的权重近似于在没有引入遗忘机制情况下的$m$个样本时的等权重均值计算。这样的增长$\mu(t)$使得均值能够跟踪上非平稳随机（均值随着时间缓慢地改变）输入$x_t$。

增量式更新样本的协方差矩阵的公式如下：$\Gamma_x^{(t)}=\frac{t-1-\mu(t)}{t}\Gamma_x^{(t-1)}+\frac{1+\mu(t)}{t}(x_t-\overline{x}^{(t)})(x_t-\overline{x}^{(t)})^T$

此处应注意，上面的协方差计算的表述中，假定样本的自由度为$t$，而不是批处理计算中的$t-1$。因为，即使有只有一个样本，相关的协方差矩阵也不应该被估计为零向量，因为$x_1$永远不会精确，这是由于物理中的测不准定理。例如，初始的方差矩阵$\Gamma_x^{(1)}$能被估计为$\delta^2 I$，而$\delta^2$为期望的数字噪声。

判别子空间与最优类别边界

每一个内部节点都维护一个最大判别特征（MDF）子空间，用它来提升模型的泛化能力。每个内部节点的MDF子空间是利用分配到该节点的样本来进行更新的。为了达到划分的目的，内部节点的每一个孩子节点都表示一个类。我们需要利用基于概率的最优类边界（基于MDF子空间）来对内部节点的输入空间进行划分。{输出空间--->MDF--->最优类边界--->划分输入空间}

A. 判别子空间

当输入的维度非常高时（一般最少几千维），为了计算效率，我们不应该在原始的输入空间$X$表示数据。更进一步，为了更好的表示特征，我们利用判别子空间$D$，也即那些与输出不相关的成分会被无视。

我们首先考虑x-clusters，每一个x-cluster是由它的统计学信息（均值与方差）刻画的。当输入状态的维度非常高时，方差是一个非常庞大的矩阵，对它进行维护是不切实际的。我们引进一个更加有效的方法：利用子空间来表示原始的状态空间。

每一个内部节点都保有q个x-clusters。这q个x-clusters的中心被表示为$C=\{c_1, c_2, ...,c_q|c_i\in X, i=1,2,...,q\}$

这些聚类的中心包含了子空间$D$在状态空间中的位置信息。$D$是一个判别空间的原因是，它的聚类形式是由输出空间$Y$决定的。

判别子空间$D$计算步骤如下：

假定第i个聚类中的样本个数为$n_i$，因些该内部节点总样本数为$n=\sum_{i=1}^q n_i$。

令$\overline{C}$为所有q个聚类的中心$\overline{C}=\frac{1}{n}\sum_{i=1}^q n_i c_i$。

从$\overline{C}$到聚类中心向量可以得到一列新的向量$s_1, s_2, ...$，如下图所示。我们称由$S=\{s_i|i=1,2,...,q\}$撑起的子空间为span(S)。

而span(S)的正交基$a_1, a_2,...,a_{q-1}$能够利用$s_1, s_2, ...,s_q$来构建（Gram-Schmidt Orthogonalization GSO 程序）：

在上面的程序中，需要注意的是，当第一步的分母为0,也即$s_1$为零时，被称为退化。与退化相关的$s_i$应该被删除掉。而GSO得到的基向量的个数是S中线性独立的元素个数。

假定，通过以上程序得到这些正交基，此刻到来一个状态向量$x\in X$，我们先计算它的点部分$s=x-\overline{C}$，然后计算$x$在判别空间上的投影 $f=M^Ts$，其中$M=[a_1, a_2, ..., a_{q-1}]$。我们称$f$为$x$在线性流形$S$上的判别特征。聚类的均值与方差都是基于判别子空间计算的。

B. 基于概率的测度

最大判别特征的子空间的一限制是，它的维度只允许为q-1。在子空间中，我们需要基于估计的概率分布自动地给每一个孩子节点确定一个最优的类边界。不同的概率分布模型得到不同的距离测度。

• negative-log-likelihood (NLL)： $G(x, c_i)=\frac{1}{2}(x-c_i)^T\Gamma_i^{-1}(x-c_i)+\frac{q-1}{2}ln(2\pi)+\frac{1}{2}ln(|\Gamma_i|)$

我们称之为状态$x$属性聚类i的Gaussian NLL。$c_i$与$\Gamma_i$为聚类样本的均值与方差矩阵。它们利用amnesic average进行更新。

• Mahalanobis NLL：

$M(x, c_i)=\frac{1}{2}(x-c_i)^T \Gamma^{-1}(x-c_i)+\frac{q-1}{2}ln(2\pi)+\frac{1}{2}ln(|F|)$

• Euclidean NLL: $E(x, c_i)=\frac{1}{2}(x-c_i)^T\rho^2 I^{-1}(x-c_i)+\frac{q-1}{2}ln(2\pi)+\frac{1}{2}ln(|\rho^2 I|)$

此处，$\Gamma$为q个聚类的平均方差矩阵:$\Gamma=\frac{1}{q-1}\sum_{i=1}^{q-1}\Gamma_i$

Euclidean NLL将所有的维度同等看待，当样本数很少时，Euclidean NLL是不错的选择。

Mahalanobis NLL利用类内点矩阵$\Gamma$。 所以，Mahalanobis NLL等于Euclidean NLL加上Fisher's LDA。每个维度的权值都不一样，参数的个数为$q(q-1)/2$，因此，Mahalanobis NLL需要更多的样本。

C. 各种测度之间的自动平滑转移

由于样本是增量式到来的，不同的距离测度的选择是依赖于每个节点中可用的样本个数。

我们将会用Euclidean NLL，当样本的个数很少时。当样本个数增加时，我们将利用Mahalanobis NLL。当聚类拥有丰富的观测样本时，我们将利用完全的Gaussian NLL。 NSPP的界限是为了限制NSPP的增长。

$\rho^2 I$的NSPP界限为$b_e=min\{(n-1)(q-1), n_s\}$

$\Gamma$的NSPP界限为$b_m=min\{max\{\frac{2(n-q)}{q}, 0\}, n_s\}$

$\Gamma_i$的NSPP界限为$b_g=\frac{2(n-q)}{q^2}$

size-dependent scatter matrix （SDSM）:$W_i=w_e\rho^2 I+w_m \Gamma+w_g\Gamma_i$

些处，$w_e = b_e/b, w_m=b_m/b, w_g=b_g/b$，$b$为归一化因子。

利用用size-dependent scatter matrix $W_i$，得到size-dependent negative log likelihood (SDNLL)为:

$L(x, c_i)=\frac{1}{2}(x-c_i)^TW_i^{-1}(x-c_i)+\frac{q-1}{2}ln(2\pi)+\frac{1}{2}ln(|W_i|)$

因为$b_e, b_m, b_g$会自动的变化，$(L(x, c_i))$会在这些合适的距离测度上平滑地切换。

D. 考虑计算效率

由于时实计算的需要，必须设计一个有效的非迭代的计算方案。

从$x$到每个x-cluster中心$c_i$的加权的平方距离矩阵如下：

$d^2(x, c_i)=(x-c_i)^TW_i^{-1}(x-c_i)$

所有节点中的状态距离的计算都是在$q-1$维的判别空间中进行的。所以SDSM为$(q-1)\times(q-1)$的对称方阵。我们只需要估计$q(q-1)/2$个参数。

为了更有效的计算上面的距离，利用Cholesky分解算法:对称方阵可以分解成$W=LL^T$

关于Cholesky分解可以参考博客矩阵分解 Cholesky分解

假设我们已经得到Cholesky分解，$W=LL^T$，且令$v=x-c_i$则$d^2(x, c_i)=v^TW_i^{-1}v=v^T(LL^T)^{-1}v=(L^{-1}v)^T(L^{-1}v)$

通过求解$Ly=v$，然后$y=L^{-1}v$以及$d^2(x,c_i)=(L^{-1}v)^T(L^{-1}v)=||y||^2$,求解$Ly=v$的时候可以用反代法，因为$L$为下三角阵。

遗留的问题

算法B update node的第（3）步中，说要对节点的最小判别子空间进行更新，是否需要重新计算最小判别子空间的基向量？ IHDR算法需要对激活叶子结点的父节点A的最小判别子空间进行更新。则，节点A中有两种空间：原始状态空间（一般维数较高，不易直接计算方差矩阵）X和最小判别子空间D。算法中只是说对节点A的最小判别子空间进行更新，直观的理解就是对基向量进行更新。但是，基向量一改变，D就变化了。之前在该空间中计算的均值与方差就不能用了。那么，在新的D中怎么恢复均值与方差信息？ 其中，均值可以通过同步的在X空间中维护一个均值来恢复。但是，方差呢？不可能通过在X空间中同步的维护一个方差信息吧？因为，最小判别子空间的初衷就是，避免在原始空间X中计算过大的方差矩阵。那么方差信息怎么恢复？ 如果利用它所属的叶子结点的x-y对来进行恢复，是否利用遗忘机制更新方差的步骤是多余的？因为，每到来一个新的数据，都会对D进行更新，方差都得重新计算。还没来得及迭代更新，方差又变了。

如果不需要重新计算基向量，那么最小判别子空间从一开始就定下来了，可以与增量式RBF网络结合。