发育网络(DN): 一个涌现的图灵机

如果一定得划分，翁教授的工作可以划分成两个时期：1）CCIPCA+IHDR；2）CCILCA+DN。IHDR作为一种偿试实现机器人大脑，显然是失败了。因为，前后两时期的工作没有太大的关联性【尽管如此，IHDR仍是一个很棒的增量式学习模型，值得深入学习与理解】。为什么IHDR不能作为机器人最终的发育模型（大脑）？翁教授的回答是：人的大脑没有一个中央的控制器，而IHDR有一个中央控制器调控所有的学习过程【IHDR是一种全局的学习模型，DN是一种局部的学习模型:::MAYBE】。结合背景，言外之意是发育网络（DN）没有中央调控。本文结合翁教授的讲座与相关论文介绍DN。首先系统的介绍发育网络的思想，然后介绍方法，接着利用一个发育网络的实现例子Where-what-network来进一步理解发育网络，最后对实验结果以及全文进行分析总结。本文例子python程序已经上传github。

文章目录

发育网络(DN): 一个涌现的图灵机
- 1. 大脑
- - 1.1 来自肌肉的上下文（Context from motor）
  - 1.2 大脑不是级联结构
  - 1.3 符号与联接的统一
  - 1.4 自主发育
- 2. 涌现通用图灵机（Emergent Universal Turing Machines）
- - 2.1 有限自动机（FA）
  - 2.2 图灵机
  - 2.3 通用图灵机
  - 2.4 涌现通用图灵机（Emergent UTM）
- 3. 发育网络（Delevelopment Network）
- - 3.1 从涌现UTM到发育网络
  - 3.2 发育网络的特点
  - 3.3 自主发育程序
  - 3.4 发育网络的例子
- 4. DN的一个例子：Where-what-network
- - 4.1 数据准备
  - - 4.1.1 构造过程说明
    - 4.1.2 样本构造程序
  - 4.1.3 构造的样本例子
  - 4.2 Where-What network介绍
  - 4.3 Where-What network训练
  - - 4.3.1 权值初始化
    - 4.3.2 响应计算
    - 4.3.3 top-k竞争响应
    - 4.3.4 参数更新
    - 4.3.5 DN利用
    - 4.3.5 DN神经元自适应增长
  - 4.4 结果
  - 4.5 程序[已上传github]
  - 4.6 分析与讨论
- 5. 结语

1. 大脑

1.1 来自肌肉的上下文（Context from motor）

翁教授认为 自主意图是从自上到下的上下文（top-down contexts）中学到的；进一步的，作为意图的上下文是来自肌肉（motor）。 下面是翁教授给的几个例子，方便理解上面那句话。

将一只刚出生的猫放在只有垂直边界线的环境中一段时间，然后将它放在正常环境中，我们发现它的V1区的神经细胞对水平边缘没有响应（具体表现就是：它不会盯着水平边缘看）。【Blakemore & Cooper, Nature 1970】
下图显示的是左视皮层的视神经细胞被右眼或左眼输入主导的情况。正常情况如下面的左图所示。相对来说，左视皮层神经细胞大多被右眼输入主导，也有不可忽略的一部分被左眼主导。现在把一个刚出生10天的小猫的右眼蒙上，经过21天后（也即出生后的第31天），情况如中间图所示。右眼的主导降低。再过6天，结果如右图所示。基本上全部由左眼主导了。【备注：正常情况下，左眼输入到右半脑，右眼输入到左半脑】

3. 将视觉神经切断，经过一段时间发现，视觉神经自动与听觉区域的神经元连接在了一起。也即，用听觉区域来“看”。【Sur, Angelucci and Sharm, Nature 1999】

1.2 大脑不是级联结构

深度学习利用深度神经网络作为学习模型，是一种级联结构。但是，大脑并不是级联的结构，各个神经元之间的突触连接方式错综复杂。

1.3 符号与联接的统一

大脑利用神经元相互连接实现知识的学习与存储，同时也具备逻辑推理能力。这说明大脑是符号主义与联接主义的统一体。

1.4 自主发育

人类从一枚受精卵发育成婴儿，从婴儿发育成成年人，这个过程是自主发育。很自然的想法是，如果给机器人造一个身体并给一个自主发育程序，机器人便会像人类一样通过自主的与环境以及其他智能体的交互进行发育，慢慢学到越来越多的知识，也越来越聪明。这也是翁教授研究自主心智发育理论的最核心驱动目标。

2. 涌现通用图灵机（Emergent Universal Turing Machines）

2.1 有限自动机（FA）

有限自动机或称为有限状态的机器，我们首先了解一下与FA非常相近的有限状态机(FSM)。有限状态机（FSM）由有限个状态及状态转换关系组成。如下图所示的有限状态机由6个状态组成，图中也明确给出状态到状态之间的转换条件。例如，当前状态为2，此时输入字符I，则转换到状态6。
然而，有限自动机由一个有限的内部状态集和一组控制规则组成，这些规则是用来控制在当前状态下读入输入符号后应转向什么状态。

【图灵机的控制就是一个有限自动机】

2.2 图灵机

图灵机是图灵在他24岁发表的论文《论可计算数及其在判定问题中的应用》中提出的一种抽象模型。该抽象模型很大程度上模仿了人类处理问题的过程：它拥有一个类似人眼和手的读写头能够读取信息以及输出信息；一条无限长的纸带，源源不断地提供信息以及供输出结果；一个类似于我们大脑的控制，能够根据问题的不同，进行不同的处理。

人类也可以抽象成图灵机。每一个会决策、会思考的人就可以被抽象的看成一个图灵机，每个人都有自己的操作系统。输入状态集合就是你所处的环境中能够看到、听到、闻到、感觉到的所有信息，可能的输出集合就是你的一言一行，以及你能够表达出来的表情动作。内部状态集合则要复杂得多。因为我们可以把任意一个神经细胞的状态组合看作是一个内部状态，那么所有可能的神经细胞的状态组合将是天文数字。这是人类的记忆。只要图灵机具有了内部状态，它就具有了记忆。

如果将图灵机与人进行类比，那么人脑就是图灵机的内部状态存储器与控制程序指令。

2.3 通用图灵机

我们可以将上面介绍的图灵机编码为字符串< M>。我们可以进一步构造出一个特殊的图灵机，它以任意图灵机M的编码< M>作为输入【图灵机以数据作为输入，通用图灵机以图灵机作为输入】，然后模拟M的运作，这样的图灵机称为通用图灵机。现代电子计算机其实就是这样一种通用图灵机的模拟，它能接受一段描述其他图灵机的程序，并运行程序实现该程序所描述的算法。

2.4 涌现通用图灵机（Emergent UTM）

上图为一个传统的有限自动机，可以表示成如下具有同样意义的符号查找表：

我们进一步编码：A-00;B-01;C-10;D-11。可得如下查找表。该查找表就是一个涌现有限自动机。

我们以Z作为状态，X作为输入。在初始状态00时，输入为01时，会到达下一个状态01，然后输入10时，又会跳到状态10，继续输入10，会跳到关态10，此时输入11，会跳入状态11。这已经非常接近本文的主角【发育网络（DN）】。

3. 发育网络（Delevelopment Network）

3.1 从涌现UTM到发育网络

还请紧盯着上面的查找表与转换例子。我们把状态Z当作人类的感知输入信息，而输入X当作人类的肌肉运动信息。我们将X与Z绑定成记忆图（发育网络）中的一个结点。Y是隐含在人类或机器人自身模型中的，不需要学习。

【备注：人类或者设计好的机器人，不需要学习他的转换关系（也即控制器）。因为，对于特定的肌肉动作，机器人或者人肯定按照自身的运动模型运行，从而自动的进入到下一步状态，在该状态感知新的信息，并输出绑定的肌肉动作信息。】

学习的过程可以描述为，对于新的感知输入信息Z【bottom】，假定机器人或者人类能够通过示教或者自主学习的方式得到当前合适的肌肉动作信息X【up】。我们将X与Z绑定在一起形成发育网络的一个记忆节点。

【备注：为了更好理解涌现图灵机与发育网络的关系，此处简化了发育网络的学习与利用过程，真实的发育网络的学习过程考虑了信息之间的相似性，是一个有效的①聚类过程 ：叶成分分析(LCA)，并按照②Hebbian 规则 来更新网络的参数。】

一个形象的例子如下图所示：【理解发育网络的重点是图灵机的控制规则是隐含在人类与机器人本体中的，不需要学习。人类与机器人本身就是一个模型，采取某一个动作，自然会显式的作用到环境】

3.2 发育网络的特点

翁教授认为机器人的自主发育程序应该具备以下八个特点：

接地地：感知物理世界，并作用到物理世界；
涌现地：内在的表述是自动生成的；
自然地：是自然编码，而不是手动编码；
增量地：不是批处理学习，没有数据集；
脑颅封闭地：人类不能改变颅内的参数；
注意地：对感知特征与电机特征具有注意机制；
有动机地：避免痛苦，寻找快东，好奇心等；
具有抽象能力地：能够从上下文中提取观念与规则；

还有一点是没有疑问的：发育网络统一了监督学习与无监督学习。

3.3 自主发育程序

通过以上介绍，我们可以大概知道发育网络的单个节点模式以及整个发育网络具有的特性。机器人需要自主发育程序自主在线构建发育网络：网络节点的增加、网络参数的更新等。翁教授提出的自主发育程序具有以下几个特点：

In-place learning （在位学习） :每个神经元通过自身内部的生理机制以及与其它神经元的交互来处理自身的学习。每个神经元都拥有同样一套发育程序，具有同等的发育能力。在位学习是以神经细胞为中心的学习。

【笔记：因此，只需要介绍清楚一个神经细胞的学习过程，就知道整个发育网络的学习，不管发育网络多么复杂。发育网络的每个神经元细胞的响应由三部分构成：a) 自底向上（上一层到当前层的输入）的响应yyy；b)自顶向下（下一层到当前层的输入）的响应aaa;c)同层间其它神经元的抑制响应hhh。总的响应为z=g(wby+wpa−whh)z=g(w_b y+w_p a -w_h h)z=g(wby+wpa−whh)，其中ggg为任意激活函数。下图绿色的圆与粗直线表示一个完整的神经元细胞，与圆相连的不带箭头实线为自底向上的输入，与粗直线相连的实线为自顶向下的输入，与粗直线相连的虚线为同层间神经元的抑制输入。】

Lateral inhibition （侧向抑制） :生物的神经网络中，例如人的视网膜中，存在着一种“侧向抑制”现象，即一个神经细胞兴奋后，通过它的分支会对周围其他神经细胞产生抑制。侧向抑制是同层的神经元竞争的一种机制。为了使侧向抑制更加有效的作用到发育网络中：a）强响应的神经元能够有效的抑制弱响应的神经元；b）弱响应的神经元不影响强响应的神经元。翁教授利用top-k竞争规则，也即同一层的神经元按响应程度从强到弱排序，排在前k的神经元有一个非零的权值更新，其他的神经元的权值不作改变。

【笔记：对于预先固定神经元数量的网络，竞争机制会使得之前所学到的模式被最近的模式替代，也即会产生遗忘。通过与翁教授的交流，也确认了这种遗忘的存在，并且翁教授认为这种遗忘与人脑的特性很像，是合理的。与翁教授的观点不同，本人坚持遗忘不是好的特性，我认为人具有的特性不一定都是好的。】

Hebbian learning（赫布学习） :赫布理论描述了突触可塑性的基本原理，即突触前神经元向突触后神经元的持续重复的刺激可以导致突触传递效能的增加。这一理论由唐纳德.赫布于1949年提出，又被称为赫布定律。赫布理论可以用于触释“联合学习”，在这种学习中，由对神经元的重复刺激，使得神经元之间的突触强度增加。这样的学习方法被称为赫布学习。赫布理论也成为非监督学习的生物学基础。

【笔记：假定当前竞争胜利的神经元权值为vjt−1v_j^{t-1}vjt−1，当前的输入为y(t)y(t)y(t)，则按以下公式对该胜利神经元进行更新：
vjt=αvjt−1+(1−α)y(t)(1)v_j^{t}=\alpha v_j^{t-1}+(1-\alpha) y(t) \tag{1}vjt=αvjt−1+(1−α)y(t)(1)
其中，1−α1-\alpha1−α为学习率，α\alphaα为遗忘率。】

Lobe components（叶成分） :叶成分分析将样本空间分为C个互相不重叠的区域，也即叶区域。顾名思义，就像树枝上的同一个芽点长出的几片树叶。一个二维空间的叶成分如下图所示：

【笔记：叶成分是由叶成分分析算法|详细推导与应用参见该博客得到的，叶成分分析是一个聚类算法，每个聚类中心由其个方向上的方向向量表示。说的再明白点，这个聚类算法中的距离是由两向量的cos(θ)cos(\theta)cos(θ)表示，其中θ\thetaθ表示两向量的夹角。我们知道，距离度量的方式对聚类算法效果影响非常大，某一距离度量可能只适用于某一些特定的数据。翁教授以叶成分分析作为距离度量，并且效果都还不错，主要原因是，它的输入为图片数据，维度非常高，在如此高维的输入空间，数据基本上是分布在一个超球面上的，所以用叶成分是合适的。但是，对于一些低维的，或者数据分布散布整个状态空间的，叶成分并不合适。如下图，是一个二维的数据集，我们无法用叶成分分析将它聚成合适的类。】

小结： 上面给出自主发育程序的主要特点，并对每个特点进行了详细的解释。侧向抑制、赫布学习以及叶成分并不是完全独立的个体，它们是你中有我，我中有你，互相依存的关系。单个神经元细胞的发育程序可以看成是叶成分分析，然而叶成分分析算法是由赫布规则以及侧向抑制胜利的聚类中心向量（神经元的权值参数）进行更新的。这是一套自成体系的算法，并且实验也验证了算法的有效性。本人对此算法抱有保持态度的地方主要有两个：1）遗忘是否合理？2）叶成分的聚类是否对所有数据分布都通用？关于聚类距离度量的选择，有一句话这么说：which one to use depends on what type of data we have and what our notion of similar is.当然，我自己也在偿试将其他距离度量方式融合到自主发育程序中。但是，当前的自主发育程序就是以叶成分分析为核心的，也即换了距离度量，就得重新推导一个新的自主发育程序。这个还是挺有挑战的。

3.4 发育网络的例子

Where-What Networks（WWNs）

4. DN的一个例子：Where-what-network

4.1 数据准备

学习时 ：Where-What-network利用图片数据（Input）以及对象的种类标签(What motor)和对象在图片的位置标签(Where motor)作为训练数据，训练得到发育网络。
利用时 ：对于训练中遇到的相似样本，发育网络能够得到对象的种类与所在图片的位置信息。

4.1.1 构造过程说明

首先介绍怎么构造Where-what-network训练所需要的数据。我们可以实景拍摄照片，然后加上位置与种类标签。但是，这个过程非常繁琐。本部分主要目的是为形象理解发育网络，真实应用性的说明居于其后。因此，本部分采用人工构造的数据来训练DN。人工构造的图片数据，可以自动的得到对象的种类与位置标签，省去了人工定标签的麻烦。下面简单给出人工构造训练数据的方法。

背景图片：为了模拟真实场景图片，一个对象一般都具有背景作为依托。例如，一只沙发上的猫，沙发就是猫的背景。本部分的背景大小为19x19，通过对上图左边的13张真实场景图片进行随机切割19x19的小块产生。
对象图片：本部分采用五种动物的大小为11x11的低分辨率图片作为对象。分别为1-猫、2-狗、3-大象、4-猪、5-卡车五种图片。
Input:随机选择一种对象图片，随机的放置于随机选择的某一背景图片的某一位置（此处有25个可选位置，上面中间图蓝色格子为可选的对象图片中心位置）。（如上图中间所示）
Where-What motors：将构造时的位置与种类作为Input的Where motor 和 What motor。

每次训练时，通过以上方式，随机得到一组训练样本：[Image,type,position],下面给出相关的样本构造的Python程序。

4.1.2 样本构造程序

1.对实景图片进行随机切割，得到19x19的背景图片。

def get_background(background_height, background_width):epsilon_error = 1.0/(256*sqrt(12))f = open('backgrounds/backgrounds_to_use.txt')splitData = f.read().split("\n")curr_num = np.random.randint(len(splitData) - 2) + 1img_file = 'backgrounds/'+splitData[curr_num]+'.'+splitData[0]curr_img = np.array(Image.open(img_file))rows, cols = curr_img.shapeul_row = np.random.randint(rows-background_height)ul_col = np.random.randint(cols-background_width)bg = curr_img[ul_row:ul_row + background_height, ul_col:ul_col + background_width]bg1 = copy.deepcopy(bg)samplemin = np.min(bg1)samplemax = np.max(bg1)bg = (bg - samplemin)/(samplemax-samplemin+epsilon_error)plt.imshow(bg)plt.show()ul_img = Image.fromarray(np.uint8(bg))ul_img.show()return bg

2.随机选择一个对象作为前景图片加入背景图片中。

def add_foreground(background, training_type):  # 01:cat; 02:dog; 03:elephant; 04:pig; 05:truckf = open('foregrounds/foregrounds_to_use.txt')splitData = f.read().split("\n")curr_num = np.random.randint(training_type) + 1img_type = curr_num - 1img_file = 'foregrounds/'+splitData[curr_num]+'.'+splitData[0] curr_img = Image.open(img_file, 'r')fg = np.array(curr_img.convert('L'))obj_width, obj_height = fg.shapefg1 = copy.deepcopy(fg)samplemin = np.min(fg1)samplemax = np.max(fg1)fg = (fg - samplemin)/(samplemax-samplemin)
#     plt.imshow(fg)
#     plt.show()p_r = np.random.randint(5)p_c = np.random.randint(5)r = (p_r)*2c = (p_c)*2position = p_r*5+p_cbackground[r:r+obj_width, c:c+obj_width] = fg
#     plt.imshow(background)
#     plt.show()return background, img_type, position

3.利用构造的输入图片，进一步得到 where and what motors，并封装成标准样本数据格式<x,y>。

def get_image(input_dim, z_neuron_num):"""1: cat; 2:dot; 3:elephant; 4:pig; 5:truckposition: (row-1)*col"""epsilon_error = 1/(256*sqrt(12))background_width = input_dim[0]background_height = input_dim[1]bg = get_background(background_height, background_height)training_type = 1if len(z_neuron_num) != 1:training_type = z_neuron_num[0]training_image, img_type, position = add_foreground(bg, training_type)true_z = []if len(z_neuron_num) != 1:true_z.append(img_type)true_z.append(position)else:true_z.append(position)true_z = np.array(true_z)return training_image, true_z

4.1.3 构造的样本例子

下图为随机选择的背景，随机选择的对象，随机安放的位置，得到最终的输出图片与Where and What motors。【备注：由于图片经过归一化，并用matplotlib的plt.imshow展示，所以显示与真实灰度图片有差异】

4.2 Where-What network介绍

本文介绍的Where-What network一个简单例子如下图所示【应该是最简单的形式了，但是麻雀虽小，该有的部分一样不少】。该网络由输入[Input]，发育网络的核心部分[DN]以及输出[Motors]三层组成。最核心的是红框标出的DN，该部分存储学习到的所有知识。并且，在本演示例子中，DN的网络神经元个数是自适应增长的，发育层的神经元个数等于25x(3^split_times)，其中，split_times为满足分裂的次数【有点像细胞的有丝分裂哈，与细胞只能一个分成两个不同，这里的一个节点可以分成我们设定的值，此处为3】。输入层的神经元个数为输入图片向量化后的长度19x19=361。What motor有5个神经元，Where motor有25个神经元。输入层与motor层的神经元个数都是直接由样本确定的。

4.3 Where-What network训练

DN的第一个特点便是in-place。所谓的in-place即，网络是局部的，单个神经元的更新只与它本身以及周围的神经元有关，并且每个神经元的更新是共用同一套发育程序【类似于人身体的细胞共用同一套DNA一样】。因此，我们只需要解释清楚DN中一个神经元的参数、输出计算以及权值更新即可。所有神经元的更新方法与步骤都是一样的。

4.3.1 权值初始化

DN神经元总共有五种参数：①bottom up weight;②top down weight;③lateral weight;④inhibit weight;⑤synapse factor。前三个参数的解释说明可见下图。我们先不考虑后两组参数，只利用前三个参数计算得到某一神经元的输出如下图所示【请仔细阅读下图】。先计算X的响应，Y的响应，然后将两者的响应加权求和得到临时的响应Temp Y，然后考虑侧向参数得到侧向神经元对该神经元的响应，并将其与Temp Y进行加权求和得到最终的神经元响应Y。

第四个参数：inhibit weight是同层抑制参数，也即当某一神经元被激活，只有附近的神经元会被相应比例激活，其他稍远的神经元则被抑制【对于当前样本，参数不被更新】。
第五个参数：synapse factor是注意力选择因子。对于不同的神经元，输入的不同维度的重要性不一样，因此，对于特定输入【X或Z】，每个Y中的神经元都附有一个因子向量，该因子向量与输入向量的维度一至，因子向量上每个元素表示对应输入向量元素的重要性，值越大，表示越重要。【这有点类似特征选择】

在本例子中，第四个参数用在top-k竞争函数中；第五个参数只用在了Y对输入X的注意力选择中，使发育网络能够自适应的选择感知数据的重要维度【特征选择】。

以下是本例中网络参数定义及初始的函数。

    def dn_create(self):self.y_lsn_flag = np.zeros((1, self.y_neuron_num))self.y_firing_age = np.zeros((1, self.y_neuron_num))self.y_inhibit_age = np.zeros((1, self.y_neuron_num))self.y_bottom_up_weight = np.ones((self.x_neuron_num, self.y_neuron_num))self.y_top_down_weight = []for i in range(self.z_area_num):self.y_top_down_weight.append(np.ones((self.z_neuron_num[i], self.y_neuron_num)))self.y_lateral_weight = np.zeros((self.y_neuron_num, self.y_neuron_num))self.y_inhibit_weight = np.ones((self.y_neuron_num, self.y_neuron_num))self.y_synapse_flag = 1"""1: only bottom-up2: bottom-up + top-down3: bottom-up + top-down + lateral4: bottom-up + top-down + lateral + inhibit"""self.y_synapse_coefficient = [0.8, 1.2, 5.0]self.y_synapse_age = 20self.y_bottom_up_synapse_diff = np.zeros(self.y_bottom_up_weight.shape)self.y_bottom_up_synapse_factor = np.ones(self.y_bottom_up_weight.shape)self.y_top_down_synapse_diff = []self.y_top_down_synapse_factor = []for i in range(self.z_area_num):self.y_top_down_synapse_diff.append(np.zeros(self.y_top_down_weight[i].shape))self.y_top_down_synapse_factor.append(np.ones(self.y_top_down_weight[i].shape))self.y_lateral_synapse_diff = np.zeros(self.y_lateral_weight.shape)self.y_lateral_synapse_factor = np.ones(self.y_lateral_weight.shape)self.y_inhibit_synapse_diff = np.zeros(self.y_inhibit_weight.shape)self.y_inhibit_synapse_factor = np.ones(self.y_inhibit_weight.shape)# z weightsself.z_bottom_up_weight = []self.z_firing_age = []for i in range(self.z_area_num):self.z_bottom_up_weight.append(np.zeros((self.y_neuron_num, self.z_neuron_num[i])))self.z_firing_age.append(np.zeros((1, self.z_neuron_num[i])))# responsesself.x_response = np.zeros((1, self.x_neuron_num))# pre lateral response is bottom up + top down, used to get lateral# pre response is bottom up + top down + lateralself.y_bottom_up_percent = 1/2self.y_top_down_percent = 1/2self.y_lateral_percent = 1/2self.y_bottom_up_response = np.zeros((1, self.y_neuron_num))self.y_top_down_response = np.zeros((self.z_area_num, self.y_neuron_num))self.y_pre_lateral_response = np.zeros((1, self.y_neuron_num))self.y_lateral_response = np.zeros((1, self.y_neuron_num))self.y_pre_response = np.zeros((1, self.y_neuron_num))self.y_response = np.zeros((1, self.y_neuron_num))self.z_response = []for i in range(self.z_area_num):self.z_response.append(np.zeros((1, self.z_neuron_num[i])))

4.3.2 响应计算

综合上图以及上节说明，一个完整的响应计算函数的输入包括：输入向量input_vec、权值向量weight_vec、注意力因子synapse_factor。首先，将注意力因子的每项元素与输入向量的对应元素相乘，再计算其归一化的单位向量，然后将归一化的融合注意力因子的输入向量与同样归一化的权值向量点乘得到最终的响应【两个归一化向量的点乘就是原向量的夹角余弦值cos(θ)cos(\theta)cos(θ)】。

def compute_response(input_vec, weight_vec, synapse_factor):"""input_vec is of shape 1x input_dimweight_vec is of shape input_dim x neuron_numsyanpse_factor is of shape input_dim x neuron_num"""_, neuron_num = weight_vec.shape_, input_dim = input_vec.shape# reshape input to neuron_num x input_dimtemp_input = np.tile(input_vec, (neuron_num, 1))temp_input = temp_input*synapse_factor.T# normalize inputtemp_input_norm = np.sqrt(np.sum(temp_input*temp_input, axis=1))temp_input_norm[temp_input_norm == 0] = 1temp_input = temp_input/np.tile(temp_input_norm.reshape(-1, 1), (1, input_dim))# normalize weightweight_vec_normalized = weight_vec*synapse_factor
#     plt.imshow(weight_vec_normalized)
#     plt.show()weight_vec_norm = np.sqrt(np.sum(weight_vec_normalized*weight_vec_normalized, axis=0))weight_vec_norm[weight_vec_norm == 0] = 1weight_vec_normalized = weight_vec_normalized/np.tile(weight_vec_norm, (input_dim, 1))output_vec = np.zeros((1, neuron_num))for i in range(neuron_num):output_vec[0, i] = np.dot(temp_input[i,:].reshape(1, -1), weight_vec_normalized[:, i].reshape(-1, 1))[0,0]return output_vec

4.3.3 top-k竞争响应

假定我们利用以上响应计算函数分别得到各输入的响应，并加权求和得到最终的神经元响应。
【函数相关代码如下】

        self.x_response = training_image.reshape(1, -1)for i in range(self.z_area_num):self.z_response[i] = np.zeros(self.z_response[i].shape)self.z_response[i][0,true_z[i]] = 1self.x_response = preprocess(self.x_response)# compute responseself.y_bottom_up_response = compute_response(self.x_response, self.y_bottom_up_weight, self.y_bottom_up_synapse_factor)for i in range(self.z_area_num):self.y_top_down_response[i] = compute_response(self.z_response[i], self.y_top_down_weight[i],self.y_top_down_synapse_factor[i])# top-down + bottom-up response    self.y_pre_lateral_response = (self.y_bottom_up_percent*self.y_bottom_up_response + self.y_top_down_percent*np.mean(self.y_top_down_response, axis=0).reshape(1,-1))/(self.y_bottom_up_percent + self.y_top_down_percent)# 侧边响应    self.y_lateral_response = compute_response(self.y_pre_lateral_response,self.y_lateral_weight,self.y_lateral_synapse_factor)self.y_pre_response = ((self.y_bottom_up_percent + self.y_top_down_percent)*self.y_pre_lateral_response + self.y_lateral_percent*self.y_lateral_response)self.y_response = top_k_competition(self.y_pre_response,self.y_top_down_response,self.y_inhibit_weight,self.y_inhibit_synapse_factor,self.y_top_k)

我们需要进一步通过top-k竞争计算激活的神经元，从而知道哪些神经元可以利用当前样本进行参数更新。top-k竞争计算就利用了inhibit weight。
程序如下：

def top_k_competition(response_input, top_down_response, inhibit_weight, inhibit_synapse_factor, top_k):"""TODO: there are two ways to do things1: if a neuron is within the synapse, then include that neuron in top-k2: if a neuron is within the synapse and the weight is > 0.5, theninclude that neuron in top-kthis version does things in the 1st wayresponse_input is of size 1xneuron_num"""response_output = np.zeros(response_input.shape)_, neuron_num = response_input.shapetop_down_flag = np.ones((1, neuron_num))for i in range(len(top_down_response)):top_down_flag = top_down_flag*top_down_response[i]for i in range(neuron_num):curr_response = response_input[0,i]curr_mask = (inhibit_synapse_factor[:,i] > 0)compare_response = response_input*curr_mask.T.reshape(1, -1)compare_response[0,i] = curr_responseneuron_id = np.argsort(-compare_response.reshape(-1))      for j in range(top_k):if len(top_down_response) != 0:if neuron_id[j] == i and top_down_flag[0,i]>0:response_output[0,i] = 1breakelif neuron_id[j] == i:response_output[0,i] = 1break            return response_output

4.3.4 参数更新

按上top-k竞争函数计算得到当前激活的神经元，就可根据赫布学习规则对各个权值进行更新。

此处需要说明两点：

本例子只对X进行了注意力选择，也即bottom_up_synapse_factor是有意义的，并且是在线更新的，通过synapse diff中间变量更新；输入Z没有注意力选择，synapse_factor所有元素都固定为1。

学习率lr是关于访问频次以及遗忘参数的函数，是随学习自适应调整的，不利用静态固定的学习率是为了合理分配历史数据与当前数据的权值，从而更好的跟踪输入环境的动态性。本例获取学习率的函数为(firing_age为该神经元的激活次数)：

程序如下：

        # hebbian learning and synapse maitenancefor i in range(self.y_neuron_num):if self.y_response[0,i] == 1: # firing neuron, currently set response to 1if self.y_lsn_flag[0,i] == 0:self.y_lsn_flag[0,i] = 1self.y_firing_age[0,i] = 0lr = get_learning_rate(self.y_firing_age[0,i]) # learning rate# bottom-up weight and synapse factorself.y_bottom_up_weight[:,i] = (1-lr)*self.y_bottom_up_weight[:,i] + lr*self.x_response.reshape(-1)self.y_bottom_up_synapse_diff[:,i] = ((1-lr)*self.y_bottom_up_synapse_diff[:,i]+lr*(np.abs(self.y_bottom_up_weight[:,i]-self.x_response.reshape(-1))))if self.y_synapse_flag>0 and self.y_firing_age[0,i] > self.y_synapse_age:self.y_bottom_up_synapse_factor[:,i] = get_synapse_factor(self.y_bottom_up_synapse_diff[:,i],self.y_bottom_up_synapse_factor[:,i],self.y_synapse_coefficient)# top-down weight and synapse factorfor j in range(self.z_area_num):self.y_top_down_weight[j][:,i] = (1-lr)*self.y_top_down_weight[j][:,i] + lr*self.z_response[j].reshape(-1)self.y_top_down_synapse_diff[j][:,i] = ((1-lr)*self.y_top_down_synapse_diff[j][:,i] + lr*np.abs(self.y_top_down_weight[j][:,i]-self.z_response[j].reshape(-1)))if (self.y_synapse_flag>1) and (self.y_firing_age[0,i]>self.y_synapse_age):self.y_top_down_synapse_factor[j][:,i] = get_synapse_factor(self.y_top_down_synapse_diff[j][:,i],self.y_top_down_synapse_factor[j][:,i],self.y_synapse_coefficient)# lateral weight and synapse factor# lateral exitation connection only exists within firing neuronsself.y_lateral_weight[:,i] = (1-lr)*self.y_lateral_weight[:,i]+lr*self.y_response.reshape(-1)self.y_lateral_synapse_diff[:,i] = ((1-lr)*self.y_lateral_synapse_diff[:,i] + lr*np.abs(self.y_lateral_weight[:,i] - self.y_response.reshape(-1)))if (self.y_synapse_flag > 2) and (self.y_firing_age[0,i]>self.y_synapse_age):self.y_lateral_synapse_factor[:,i] = get_synapse_factor(self.y_lateral_synapse_diff[:,i],self.y_lateral_synapse_factor[:,i],self.y_synapse_coefficient)self.y_firing_age[0,i] = self.y_firing_age[0,i] + 1elif self.y_lsn_flag[0,i] == 0: # initialization stage neuron is always updatinglr = get_learning_rate(self.y_firing_age[0,i])normed_input = self.x_response.reshape(-1, 1)*self.y_bottom_up_synapse_factor[:,i].reshape(-1, 1)self.y_bottom_up_weight[:,i] = (1-lr)*self.y_bottom_up_weight[:,i]+lr*normed_input.reshape(-1)self.y_bottom_up_weight[:,i] = self.y_bottom_up_weight[:,i]*self.y_bottom_up_synapse_factor[:,i]self.y_bottom_up_synapse_diff[:,i] = ((1-lr)*self.y_bottom_up_synapse_diff[:,i] +lr*np.abs(self.y_bottom_up_weight[:,i] - normed_input.reshape(-1)))if self.y_synapse_flag>0 and self.y_firing_age[0,i] > self.y_synapse_age:self.y_bottom_up_synapse_factor[:,i] = get_synapse_factor(self.y_bottom_up_synapse_diff[:,i],self.y_bottom_up_synapse_factor[:,i],self.y_synapse_coefficient)  # top-down weight and synapse factorfor j in range(self.z_area_num):self.y_top_down_weight[j][:,i] = (1-lr)*self.y_top_down_weight[j][:,i] + lr*self.z_response[j].reshape(-1)self.y_top_down_synapse_diff[j][:,i] = ((1-lr)*self.y_top_down_synapse_diff[j][:,i] + lr*np.abs(self.y_top_down_weight[j][:,i]-self.z_response[j].reshape(-1)))if (self.y_synapse_flag>1) and (self.y_firing_age[0,i]>self.y_synapse_age):self.y_top_down_synapse_factor[j][:,i] = get_synapse_factor(self.y_top_down_synapse_diff[j][:,i],self.y_top_down_synapse_factor[j][:,i],self.y_synapse_coefficient)# lateral weight and synapse factor# lateral exitation connection only exists within firing neuronsself.y_lateral_weight[:,i] = (1-lr)*self.y_lateral_weight[:,i]+lr*self.y_response.reshape(-1)self.y_lateral_synapse_diff[:,i] = ((1-lr)*self.y_lateral_synapse_diff[:,i] + lr*np.abs(self.y_lateral_weight[:,i] - self.y_response.reshape(-1)))if (self.y_synapse_flag > 2) and (self.y_firing_age[0,i]>self.y_synapse_age):self.y_lateral_synapse_factor[:,i] = get_synapse_factor(self.y_lateral_synapse_diff[:,i],self.y_lateral_synapse_factor[:,i],self.y_synapse_coefficient)   else:lr = get_learning_rate(self.y_inhibit_age[0,i])temp = np.zeros(self.y_inhibit_synapse_factor.shape)for j in range(self.y_neuron_num):temp[:,j] = self.y_pre_lateral_response.reshape(-1)*self.y_inhibit_synapse_factor[:,j]temp[:,j] = (temp[:,j] > self.y_pre_lateral_response[0,i])self.y_inhibit_weight[:,i] = (1-lr)*self.y_inhibit_weight[:,i]+lr*temp[:,i]self.y_inhibit_synapse_diff[:,i] = ((1-lr)*self.y_inhibit_synapse_diff[:,i] + lr*np.abs(self.y_inhibit_weight[:,i]-temp[:,i]))if (self.y_synapse_flag > 3) and (self.y_firing_age[0,i]>self.y_synapse_age):self.y_inhibit_synapse_factor[:,i] = get_synapse_factor(self.y_inhibit_synapse_diff[:,i],self.y_inhibit_synapse_factor[:,i],self.y_synapse_coefficient)self.y_inhibit_age[0,i] = self.y_inhibit_age[0,i] + 1## z neuron learningfor area_idx in range(self.z_area_num):for i in range(self.z_neuron_num[area_idx]):if self.z_response[area_idx][0,i] == 1:lr = get_learning_rate(self.z_firing_age[area_idx][0,i])self.z_bottom_up_weight[area_idx][:,i] = (1-lr)*self.z_bottom_up_weight[area_idx][:,i]+lr*self.y_response.reshape(-1)self.z_firing_age[area_idx][0,i] = self.z_firing_age[area_idx][0,i] + 1

4.3.5 DN利用

如果，我们觉得DN训练的可以，或者中途测试DN的学习效果，我们需要看看DN对于某一输入X，输出的where motor与what motor的准确与否。如果准确率很高，说明DN对于当前样本的学习效果很好，可以用来预测。

已知输入X，计算输出Z的程序如下：

    def dn_test(self, test_image):self.x_response = test_image.reshape(1, -1)self.x_response = preprocess(self.x_response)self.y_bottom_up_response = compute_response(self.x_response, self.y_bottom_up_weight, self.y_bottom_up_synapse_factor)self.y_pre_lateral_response = self.y_bottom_up_responseself.y_lateral_response = compute_response(self.y_pre_lateral_response,self.y_lateral_weight,self.y_lateral_synapse_factor)self.y_pre_response = ((self.y_bottom_up_percent*self.y_pre_lateral_response+self.y_lateral_percent*self.y_lateral_response)/(self.y_bottom_up_percent+self.y_lateral_percent))self.y_response = top_k_competition(self.y_pre_response, [], self.y_inhibit_weight, self.y_inhibit_synapse_factor, self.y_top_k)z_output = []for i in range(self.z_area_num):self.z_response[i] = compute_response(self.y_response,self.z_bottom_up_weight[i],np.ones(self.z_bottom_up_weight[i].shape))z_output_i = np.argmax(self.z_response[i])z_output.append(z_output_i)return np.array(z_output)

4.3.5 DN神经元自适应增长

本例中的发育程序会自适应的对发育网络的神经元个数进行增加，当某一神经元被连续激活超过某一次数时，程序中为20,则对神经元进行有丝分裂，一个神经元分裂成三个。
分裂函数如下：

def dn_split(dn, split_num, split_firing_age):input_dim = dn.x_neuron_numy_top_k = dn.y_top_kz_neuron_num = dn.z_neuron_numy_neuron_num = dn.y_neuron_num*split_numnew_to_old_index = np.zeros(y_neuron_num, dtype=np.int)for i in range(dn.y_neuron_num):start_ind = i*split_numend_ind = (i+1)*split_numnew_to_old_index[start_ind:end_ind] = inew_to_old_index = new_to_old_index.tolist()new_dn = DN(input_dim, y_neuron_num, y_top_k, z_neuron_num)for i in range(new_dn.y_neuron_num):j = new_to_old_index[i]new_dn.y_lsn_flag[0,i] = dn.y_lsn_flag[0,j]new_dn.y_firing_age[0,i] = split_firing_agenew_dn.y_inhibit_age[0,i] = split_firing_agenew_dn.y_bottom_up_weight[:,i] = (dn.y_bottom_up_weight[:,j] + generate_rand_mutate(dn.y_bottom_up_weight[:,j].shape))new_dn.y_bottom_up_weight[:,i] = (dn.y_bottom_up_weight[:,j])/np.max(new_dn.y_bottom_up_weight[:,j])for z_ind in range(new_dn.z_area_num):new_dn.y_top_down_weight[z_ind][:,i] = (dn.y_top_down_weight[z_ind][:,j]+generate_rand_mutate(dn.y_top_down_weight[z_ind][:,j].shape))new_dn.y_top_down_weight[z_ind][:,i] = new_dn.y_top_down_weight[z_ind][:,i]/np.max(new_dn.y_top_down_weight[z_ind][:,i])new_dn.y_lateral_weight[:,i] = dn.y_lateral_weight[new_to_old_index, j]new_dn.y_inhibit_weight[:,i] = (dn.y_inhibit_weight[new_to_old_index, j]+generate_rand_mutate(dn.y_inhibit_weight[new_to_old_index, j].shape))new_dn.y_inhibit_weight[:,i] = new_dn.y_inhibit_weight[:,i]/np.max(new_dn.y_inhibit_weight[:,i])new_dn.y_bottom_up_synapse_diff[:,i] = dn.y_bottom_up_synapse_diff[:,j]new_dn.y_bottom_up_synapse_factor[:,i] = np.ones(dn.y_bottom_up_synapse_factor[:,j].shape)for z_ind in range(new_dn.z_area_num):new_dn.y_top_down_synapse_diff[z_ind][:,i] = dn.y_top_down_synapse_diff[z_ind][:,j]new_dn.y_top_down_synapse_factor[z_ind][:,i] = np.ones(dn.y_top_down_synapse_factor[z_ind][:,j].shape)new_dn.y_lateral_synapse_diff[:,i] = dn.y_lateral_synapse_diff[new_to_old_index, j]new_dn.y_lateral_synapse_factor[:,i] = np.ones(dn.y_lateral_synapse_factor[new_to_old_index, j].shape)new_dn.y_inhibit_synapse_diff[:,i] = dn.y_inhibit_synapse_diff[new_to_old_index, j]new_dn.y_inhibit_synapse_factor[:,i] = np.ones(dn.y_inhibit_synapse_factor[new_to_old_index,j].shape)for z_ind in range(new_dn.z_area_num):new_dn.z_bottom_up_weight[z_ind][i,:] = dn.z_bottom_up_weight[z_ind][j,:]for z_ind in range(new_dn.z_area_num):new_dn.z_firing_age[z_ind] = np.ones(dn.z_firing_age[z_ind].shape)return new_dn

4.4 结果

每步的训练效果如下所示，我们发现中间有三次分裂操作，最终DN神经元个数为2533=225，最终的测试准确率为[0.99 0.962]。我们注意会发现，每当网络的性能遇到瓶颈时，发生分裂后，网络的性能重新被提升了。
500 training, current performance:
[0.226 0.054]
1000 training, current performance:
[0.278 0.068]
1500 training, current performance:
[0.362 0.164]
2000 training, current performance:
[0.41 0.17]
2500 training, current performance:
[0.4 0.172]
3000 training, current performance:
[0.442 0.182]
3500 training, current performance:
[0.476 0.212]
4000 training, current performance:
[0.48 0.226]
4500 training, current performance:
[0.486 0.206]
5000 training, current performance:
[0.466 0.242]
5500 training, current performance:
[0.462 0.204]
6000 training, current performance:
[0.48 0.232]
6500 training, current performance:
[0.436 0.174]
7000 training, current performance:
[0.522 0.256]
splitting at 7263
7500 training, current performance:
[0.444 0.292]
8000 training, current performance:
[0.776 0.558]
8500 training, current performance:
[0.8 0.558]
9000 training, current performance:
[0.78 0.59]
9500 training, current performance:
[0.8 0.574]
10000 training, current performance:
[0.784 0.6 ]
10500 training, current performance:
[0.794 0.586]
11000 training, current performance:
[0.82 0.596]
11500 training, current performance:
[0.824 0.614]
12000 training, current performance:
[0.826 0.602]
12500 training, current performance:
[0.818 0.58 ]
splitting at 12804
13000 training, current performance:
[0.688 0.444]
13500 training, current performance:
[0.802 0.696]
14000 training, current performance:
[0.922 0.882]
14500 training, current performance:
[0.934 0.908]
15000 training, current performance:
[0.92 0.9 ]
15500 training, current performance:
[0.924 0.876]
16000 training, current performance:
[0.96 0.922]
16500 training, current performance:
[0.974 0.928]
17000 training, current performance:
[0.988 0.954]
17500 training, current performance:
[0.97 0.956]
18000 training, current performance:
[0.984 0.95 ]
18500 training, current performance:
[0.992 0.984]
19000 training, current performance:
[0.988 0.956]
19500 training, current performance:
[0.988 0.97 ]
20000 training, current performance:
[0.996 0.97 ]
20500 training, current performance:
[0.986 0.958]
21000 training, current performance:
[0.99 0.956]
21500 training, current performance:
[0.994 0.972]
22000 training, current performance:
[0.994 0.978]
22500 training, current performance:
[0.992 0.982]
23000 training, current performance:
[0.992 0.978]
23500 training, current performance:
[0.992 0.966]
24000 training, current performance:
[0.99 0.958]
24500 training, current performance:
[0.99 0.962]
testing error is [0.996 0.975]
DN y_bottom_up_weight按15x15的矩阵点打印输出如下，每个矩阵点是一个19x19的小矩阵。

4.5 程序[已上传github]

声明：本例子的程序是参照翁教授的学生Zejia Zheng的matlab代码写的。
附上源码地址github/ZejiaZheng/multiresolutionDN[matlab]

本人通过将matlab代码改写成python代码的方式来加深入理解发育网络。本人是基于jupyter notebook交互式编译框架下写的python程序，现已上传到github。
github/huyanmingtoby/multiresolutionDN[python]

4.6 分析与讨论

本节（第4节）利用一个实现例子【Where-What-Network】来补充说明发育网络。利用构造的图片数据以及对象的位置与种类标答作为WWN的输入，对网络进行训练。在训练过程中，发育网络自己对发育层的神经元进行两次有丝分裂，最终形成225个神经元。每次分裂都能使测试的正确率大幅度提升。最后还将y_bottom_up_weight打印了出来，每个神经元与X输入相关的权值按19x19打印，发现它们是一系列与输入的图片非常相似的模式。对于所有对象可能出现的所有位置为5x25=125。最理想的神经元个数为125便能把所有可能的模式存储。初始化时DN的神经元个数为25个，显然不够；经过一次有丝分裂变成25x3=75个，也还不够；再经过一次有丝分裂75x3=225个，理论上是够了【实际实验结果也表明是够了，测试准确率非常高[0.99 0.962]】。但是，有一点需要说明，那就是背景是随机选的，对于同一个对象在同一个位置放入随机的背景中，那还是无数多种。最终发育网络的神经元在225时的效果就已经很好了，说明发育网络能够将抓住重要的信息（对象相关的像素点），忽略那些背景像素点，这是注意力选择起作用了。神经元注意力选择权值让其只对输入的对应区域有响应，从而避免了背景噪声的影响。

个人观点：我自己看论文，特别是机器学习方面的论文，看完感觉懂了，实则不然，很多重要的思想都掩藏在真实的实现过程中。如果真得觉得手头那篇论文非常有价值，最好能找到作者实现的程序读一读，并自己动手编程实现。很多难以言表却又重要的思想，作者为了可读性可能不会写进论文，但是肯定藏在程序里。如果觉得发育网络很有价值，建议完整程序git下来，自己看一遍，动手写一写，跑一跑程序。

5. 结语

本文在思想、方法、算法、程序等四个层面对发育网络进行了系统的介绍。发育网络具有in-place特点，每个神经元包含全部的发育程序，因此，我们着重对一个神经元的初始化、训练、利用进行说明。每个神经元的发育程序都是一个CCILCA，本文为了描述清楚DN的思想，将叶成分分析拆分开了，很难看到完成LCA的身影。但我认为这对于DN思想与方法的介绍是有益的。我将在另一篇博客将本文打碎的LCA还原成完整的算法。最后，利用发育网络的一个列子WWN加深DN的理解，并借此透露更多DN的细节。DN的发育层神经元个数自适应调整，最终的学习效果也是非常不错的。
另外，本文也掺入大量个人见解。不一定合适，欢迎交流。现将本人对DN的疑惑整理如下：

DN中的遗忘是否合理，人类具有的特质是否就是好的？

能否在不大改自主发育程序的前提下，用其他距离度量替代cos(θ)cos(\theta)cos(θ)？

DN归根结底还是以查找表的形式存储学习到的知识，也即<x, z>，人类大脑中的知识是否真的就是这种以一个个记忆状态存在？像深度学习、RBF网络会用非线性函数拟合x到z的关系，这显然是更加有效率的方式【先抛开人类增量式、序列式学习】，这种函数表征的方式能够对信息进行更有效的编码，从而使得同样的信息量，函数表示所需要的模型容量会更小。