精华内容
下载资源
问答
  • 主要功能有 SimFix.me SimGlobal.m testf.me testg.m 是测试特定值的文件 项目相关文档请咨询 http://www.dicecca.net/giovanni/biblo/calcolo/ 项目由乔瓦尼·迪切卡弗吉尼亚·贝利诺
  • 多层感知机解决了之前无法模拟异或逻辑的缺陷,同时更多的层数也让网络更能够刻画现实世界中的复杂情形。...但是随着神经网络层数的加深,优化函数越来越容易陷入局部最优解(即过拟合,在训练样本上有很好的拟...

    多层感知机解决了之前无法模拟异或逻辑的缺陷,同时更多的层数也让网络更能够刻画现实世界中的复杂情形。理论上而言,参数越多的模型复杂度越高,“容量”也就越大,也就意味着它能完成更复杂的学习任务。多层感知机给我们带来的启示是,神经网络的层数直接决定了它对现实的刻画能力——利用每层更少的神经元拟合更加复杂的函数。但是随着神经网络层数的加深,优化函数越来越容易陷入局部最优解(即过拟合,在训练样本上有很好的拟合效果,但是在测试集上效果很差),并且这个“陷阱”越来越偏离真正的全局最优。利用有限数据训练的深层网络,性能还不如较浅层网络。同时,另一个不可忽略的问题是随着网络层数增加,“梯度消失”(或者说是梯度发散diverge)现象更加严重。具体来说,我们常常使用sigmoid作为神经元的输入输出函数。对于幅度为1的信号,在BP反向传播梯度时,每传递一层,梯度衰减为原来的0.25。层数一多,梯度指数衰减后低层基本上接受不到有效的训练信号。那么深度学习中是如何解决局部极值及梯度消失问题的呢?

    根据我个人的理解,第一个阶段就是2006年Hinton提出的逐层预训练方法,为了解决深层神经网络的训练问题,一种有效的手段是采取无监督逐层训练(unsupervised layer-wise training),其基本思想是每次训练一层隐节点,训练时将上一层隐节点的输出作为输入,而本层隐节点的输出作为下一层隐节点的输入,这被称之为“预训练”(pre-training);在预训练完成后,再对整个网络进行“微调”(fine-tunning)训练。比如Hinton在深度信念网络(Deep Belief Networks,简称DBN)中,每层都是一个RBM,即整个网络可以被视为是若干个RBM堆叠而成。在使用无监督训练时,首先训练第一层,这是关于训练样本的RBM模型,可按标准的RBM进行训练;然后,将第一层预训练号的隐节点视为第二层的输入节点,对第二层进行预训练;... 各层预训练完成后,再利用BP算法对整个网络进行训练。虽然也解决了一些问题,但并没有特别火。

    事实上,“预训练+微调”的训练方式可被视为是将大量参数分组,对每组先找到局部看起来较好的设置,然后再基于这些局部较优的结果联合起来进行全局寻优。这样就在利用了模型大量参数所提供的自由度的同时,有效地节省了训练开销。
           第二个阶段开始的标志就是2012年IMAGENET比赛中,CNN以压倒性优势取得胜利,自此开始深度学习才真正引人关注起来。虽然都叫深度学习,但其侧重点完全不同,通过一些手段,比如relu, dropout等小技巧,第二波深度学习算法已经基本抛弃了预训练的做法。caffe也是第二波深度学习的产物,所以里面除了自编码网络基本没有逐层预训练这种东西。

    引入RELU代替sigmoid激活函数。为什么引入Relu呢? 第一,采用sigmoid等函数,算激活函数时(指数运算),计算量大,反向传播求误差梯度时,求导涉及除法,计算量相对大,而采用Relu激活函数,整个过程的计算量节省很多。第二,对于深层网络,sigmoid函数反向传播时,很容易就会出现梯度消失的情况(在sigmoid接近饱和区时,变换太缓慢,导数趋于0,这种情况会造成信息丢失),从而无法完成深层网络的训练。第三,Relu会使一部分神经元的输出为0,这样就造成了网络的稀疏性,并且减少了参数的相互依存关系,缓解了过拟合问题的发生。
           当然现在也有一些对relu的改进,比如prelu,random relu等,在不同的数据集上会有一些训练速度上或者准确率上的改进,具体的大家可以找相关的paper看。多加一句,现在主流的做法,会在做完relu之后,加一步batch normalization,尽可能保证每一层网络的输入具有相同的分布。

    传统Sigmoid系激活函数

    深度学习解决局部极值和梯度消失问题方法简析

       稀疏性激活函数ReLU

    深度学习解决局部极值和梯度消失问题方法简析

    Dropout

        训练神经网络模型时,如果训练样本较少,为了防止模型过拟合,Dropout可以作为一种trikc供选择。Dropouthintion最近2年提出的,源于其文章Improving neural networks by preventing co-adaptation of feature detectors.中文大意为:通过阻止特征检测器的共同作用来提高神经网络的性能。

           Dropout是指在模型训练时随机让网络某些隐含层节点的权重不工作,不工作的那些节点可以暂时认为不是网络结构的一部分,但是它的权重得保留下来(只是暂时不更新而已),因为下次样本输入时它可能又得工作了。

      按照hinton的文章,他使用Dropout时训练阶段和测试阶段做了如下操作:

      在样本的训练阶段,在没有采用pre-training的网络时(Dropout当然可以结合pre-training一起使用),hintion并不是像通常那样对权值采用L2范数惩罚,而是对每个隐含节点的权值L2范数设置一个上限bound,当训练过程中如果该节点不满足bound约束,则用该bound值对权值进行一个规范化操作(即同时除以该L2范数值),说是这样可以让权值更新初始的时候有个大的学习率供衰减,并且可以搜索更多的权值空间。

      在模型的测试阶段,使用”mean network(均值网络)”来得到隐含层的输出,其实就是在网络前向传播到输出层前时隐含层节点的输出值都要减半(如果dropout的比例为50%),其理由文章说了一些,可以去查看。

      关于Dropout,文章中没有给出任何数学解释,Hintion的直观解释和理由如下:

      1. 由于每次用输入网络的样本进行权值更新时,隐含节点都是以一定概率随机出现,因此不能保证每2个隐含节点每次都同时出现,这样权值的更新不再依赖于有固定关系隐含节点的共同作用,阻止了某些特征仅仅在其它特定特征下才有效果的情况。

      2. 可以将dropout看作是模型平均的一种。对于每次输入到网络中的样本(可能是一个样本,也可能是一个batch的样本),其对应的网络结构都是不同的,但所有的这些不同的网络结构又同时share隐含节点的权值。这样不同的样本就对应不同的模型,是bagging的一种极端情况。个人感觉这个解释稍微靠谱些,和baggingboosting理论有点像,但又不完全相同。

      3. native bayesdropout的一个特例。Native bayes有个错误的前提,即假设各个特征之间相互独立,这样在训练样本比较少的情况下,单独对每个特征进行学习,测试时将所有的特征都相乘,且在实际应用时效果还不错。而Droput每次不是训练一个特征,而是一部分隐含层特征。

      4. 还有一个比较有意思的解释是,Dropout类似于性别在生物进化中的角色,物种为了使适应不断变化的环境,性别的出现有效的阻止了过拟合,即避免环境改变时物种可能面临的灭亡。

      文章最后当然是show了一大把的实验来说明dropout可以阻止过拟合。这些实验都是些常见的benchmark,比如Mnist, Timit, Reuters, CIFAR-10, ImageNet.

       后来又衍生出Droput的改进版maxoutDropConnect

          第三个阶段可能就算是最近的高速公路网络(highway network)和深度残差学习(deep residual learning)进一步避免了梯度消失,网络层数达到了前所未有的一百多层(深度残差学习:152层)!

          对于deep learning越深越难训的问题,微软残差学习这个工作的直观想法应该就是说既然太深了难训,那我就把前面的层,夸几层直接接到后面去(短接),会不会既保证了深度,又让前面的好训一些呢?其实这个并不是特别新的想法,这个训得这么好,应该也是经过了很辛苦和仔细的调参、初始化等工作。换一个角度看,这也是一种把高阶特征和低阶特征再做融合,从而得到更好的效果的思路。其实细看会发现,它是highway networks的一个特例。 最主要的贡献,应该是把这种很深的highway networks在主流的benchmark上做了一遍,并且实践证明:1.这么深的网络效果比浅的好多了;2 highway networks能够在很深的网络训练中还能保持梯度稳定,不易消失或爆炸。相信不久后会有一堆工作出现类似的网络结构。residual block 其实就是 transform function T 为两个卷积和一个relu, transform gate是一个data independent 值为0.5 的一个 highway layer...残差只是 kaiming 的换了种说法的马甲..

          下面是复制的highway networkdeep residual learning论文解析,由于本人是初学者,很多地方还不是很理解,anyway,先mark一下吧,方便今后研读吧。。。

    Deep Residual Learning for Image Recognition

      首先第一篇肯定要谈,开创了152层 deep network,成功摘取 ILSVRC2015 全(主要)类别第一名 的工作,来自 MSRA 的 Deep Residual Network 。虽然 residual network 的名字很新,但其实可以把它看成是Highway Networks 的特例 Highway Networks 的介绍见下)。尽管如此,作者在解释 residual 的 motivation 时还是非常充分合理,让人不得不佩服其在背后自己的思考。

    深度学习解决局部极值和梯度消失问题方法简析

      作者提出 residual network 的 motivation 其实依然是 information flow 在 deep network 中受阻。大家都可以想象到这个问题。但是这篇工作中,作者是如何验证这个问题的呢?他们用了两种角度,第一种是在 Introduction 里就提到的:当他们用 deep network 和可类比的 shallower network 横向对比(如上图),发现 deep network 的 training error 总是比 shallower network 要高。可是这在理论上是不对的。因为如果一个 shallow network 可以被训练优化求解到某一个很好的解(subspace of that of deep network),那么它对应的 deep 版本的 network 至少也可以,而不是更差。但事实并不是如此。这种 deep network 反而优化变差的事情,被作者称为 “degration problem”。第二种 角度是在实验部分提到的,从理论角度分析这种 degration 不应该是 gradient vanishing 那种问题,而是一种真正的优化困难。于是,为了解决这个 degration problem,作者提出了 residual network,意思是说,如果我们直接去逼近一个 desired underlying mapping (function) 不好逼近(优化困难,information flow 搞不定),我们去让一个 x 和 mapping 之间的 residual 逼近 0 应该相对容易。这就是作者给 highway networks 找的 residual 的解释。

    深度学习解决局部极值和梯度消失问题方法简析

      那么,在实际上,residual network block(上图)就相当于是 Highway network block 中 transform gate 和 transform function 被特例后的结果——因为是特例了,也自然 reduce 了 parameter,也算是这篇工作中作者一个卖点。现在,问题就是,它为啥是 Highway Networks 的特例呢?介绍完 Highway Networks 就明白了。

    Training Very Deep Networks

      好了,现在让我们来再介绍一遍 Highway Networks 这个工作。这篇论文前身是《Highway Networks》,发表于 ICML workshop 上。最初放在 arXiv 上,现在已经被 NIPS'15 接收。这个工作纯被 LSTM 的 gate 启发,既然 gate 也是为了解决 Information flow,有没有其他方式去解决?更直观一点的,不通过 gradient 的?既然 information 像被阻隔了一样,我们就暴力让它通过一下,给它们来个特权——在某些 gate 上,你就不用接受审查transform)了,直接通过吧。这像高速公路一样——于是就有了这个名字,Highway NetworksHW-Nets)。

       To overcome this, we take inspiration from Long Short Term Memory (LSTM) recurrent networks. We propose to modify the architecture of very deep feedforward networks such that information flow across layers becomes much easier. This is accomplished through an LSTM-inspired adaptive gating mechanism that allows for paths along which information can flow across many layers without attenuation. We call such paths information highways. They yield highway networks, as opposed to traditional ‘plain’ networks.

      加粗了 adaptive,这就是这个 mechanism 的重点  也是之所以说 Deep Residual Network 是它的一个特例的原因所在  在文章中,公式(2-3)就是他们的机制。

      公式(3)是公式(2)变形。 核心是 transform function H 和 transform gate T 。这样,对于当前的这个input,在这个 layer 里,公式(3)决定多大程度是让它去进行 nonlinear transform(隐层),还是多大程度上让它保留原貌 直接传递给下一层,直通无阻。

    深度学习解决局部极值和梯度消失问题方法简析

      那么,在这里,其实我们也可以把这个公式(3)像 residual block 一样,拆成两个 component,一个是 H,一个是 x。如果是 x,就变成了 residual block 中 identity mapping 的 skip connection。这是一个 intuition 的解释。那么再具体等价一下,Deep Residual Network 那篇工作,讲到,自己比 Highway Networks 的一个优势是,parameter 少,并且自己的“gate”永远不 close。这两个优势,既都对,也都不对。 关于第一点 ,这是事实,而这恰恰是把 Highway Networks 的一个优势给抹掉了。在 Highway Networks 文章的 4.1 部分,有讨论自己这种 adaptive mechansim for information flow learning 的优点。也就是说,如果说 Highway Networks 是 data-independent 的,那么 Deep Residual Network 就是 data-dependent 的,不 flexible,更可能不最优Deep Residual Network 就是 Highway Networks 中提到的 hard-wired shortcut connections), 它只是一种 p=0.5 的二项分布的期望(np  关于第二点 ,就不一定是事实了。因为实际上,虽然公式(3)看起来,transform gate 有可能完全 close,然而 transform function 却永远不可能为 0 或者 1,它的取值范围只是 (0,1) 开区间,所以在这点上它和 Deep Residual Network 是一样的。

      对比说完,继续介绍 Highway NetworksHighway Networks 更巧妙的是,它在公式(3)这种 layerwise 的基础上,进一步拆分, 把每一个 layer 的处理从 layerwise 变成了 blockwise ,也就是我们可以对于 input 中的某一部分(某些维度)选择是 transform 还是 carry。这个 blockwise 的实现也在实验中被验证特别有用,大家可以去看下,我就不重复粘贴了。

      最后说一点关于CNNRNN的个人感觉吧, CNN卷积神经网络是一种特殊的深层的神经网络模型,它的特殊性体现在两个方面,一方面它的神经元间的连接是非全连接的,另一方面同一层中某些神经元之间的连接的权重是共享的(即相同的)。它的非全连接和权值共享的网络结构使之更类似于生物神经网络,降低了网络模型的复杂度(对于很难学习的深层结构来说,这是非常重要的),减少了权值的数量。上述就是所谓的局部感受野和权值共享。对卷积和池化操作,感觉可以使CNN具有原生态的稀疏性,具有一定的泛化能力和鲁棒性。RNN可以看成一个在时间上传递的神经网络,它的深度是时间的长度!正如我们上面所说,“梯度消失现象又要出现了,只不过这次发生在时间轴上对于t时刻来说,它产生的梯度在时间轴上向历史传播几层之后就消失了,根本就无法影响太遥远的过去。因此,之前说所有历史共同作用只是理想的情况,在实际中,这种影响也就只能维持若干个时间戳。为了解决时间上的梯度消失,机器学习领域发展出了长短时记忆单元LSTM,通过门的开关实现时间上记忆功能,并防止梯度消失

          以上是我最近阅读深度学习博客总结的深度学习解决局部极值和梯度消失问题的方法,本人初学者,理解总结的不一定准确,如有错误以后再纠正吧。。。

    展开全文
  • 本文提出了一种利用不断调大学习率的方法试图跳出SGD优化过程中的局部极小值或者鞍点的方法,并提出了两种具体的方法:随机漫步起跳法历史最优起跳法,实验证明相对常规优化方法有一定性能提升。

    /* 版权声明:可以任意转载,转载时请标明文章原始出处和作者信息 .*/

                                                  

                                                                                             黄通文   张俊林( 2016年12月)

                                                                


    注:这篇文章的思路以及本文是我们在2016年底左右做自动发现探索网络结构过程中做的,当时做完发现ICLR 2017有类似论文提出相似的方法,所以没有做更多实验并就把这篇东西搁置了。今天偶然翻到,感觉还是有一定价值,所以现在发出来。


    对于神经网络目标函数优化来说,在参数寻优过程中最常用的算法是梯度下降,目前也出现了很多基于SGD基础上的改进算法,比如带动量的SGD,Adadelta,Adagrad,Adam,RMSProp等梯度下降改进方法,这些算法大多都是针对基础更新公式进行改进:

       


          其中,w是训练过程中的当前参数值,alpha是学习率,grad(w)是此刻的梯度方向。目前的改进算法要么修正学习率alpha要么修正梯度计算方法grad(w),要么两者都修正,比如Adagrad/Adam动态调整学习率,Adam动态调整梯度的更新量或者有些算法加入动量部分。


     一般深度神经网络由于很深的深度以及非线性函数这两个主要因素,导致目标的损失函数是非凸函数,该损失函数曲面中包含许多的驻点,驻点可能是以下三类关键点中的一种:鞍点,局部极小值,全局极小值,其中全局极小值是我们期望训练算法能够找到的,而其它两类是期望能够避免的。一个示意性的损失函数曲面如下图所示:



    不少学习算法在训练过程中,随着误差的减少,迭代次数的增加,步长变化越来越小,训练误差越来越小直到变为零,局部探索到一个驻点,如果这个点是误差曲面的鞍点,即不是最大值也不是最小值的平滑曲面,那么一般结果表现为性能比较差;如果这个驻点是局部极小值,那么表现为性能较好,但不是全局最优值。目前大多数训练算法在碰到无论是鞍点还是局部极值点的时候,因为此刻学习率已经变得非常小,所以会陷入非全局最优值不能自拔。这其实就是目前采用SGD类似的一阶导数方法训练机器学习系统面临的一个很重要的问题。


    目前来看,没有方法能够保证找到非凸函数优化过程中的全局最小值,也就是说,不论你怎么调参将性能调到你能做到的最好,仍然大概率获得了优化过程中的一个鞍点或者局部极小值,那么一个很自然的想法就是:能否在优化过程的末尾,就是已经大概率陷入某个局部极小值或者鞍点的时候,放大学习率,强行让训练算法跳出当前的局部最小值或者鞍点,所谓“世界这么大,我想去看看”,继续寻找可能效果会更好的极值点,这样多次跳出优化过程中的局部最小值,然后比较其中性能最好的极小值,并以那个时刻的参数作为模型参数。当然,这仍然没有办法保证找到全局最小值,只是说多比较一些极值点,矮子里面拔将军的意思。这个思路是去年我们在利用人工生命和遗传算法探索自动生成神经网络结构的不断优化网络性能过程中产生的,当时是看到了SNAPSHOT ENSEMBLES : TRAIN 1,GET M FOR FREE这篇论文受到启发动了这个念头,于是单独测试了一下这个思路,结果表明一般情况下,这种方法是能够有效提升优化效果的,不过提升幅度不算太大(后来看到了这篇论文:SGDR: STOCHASTIC GRADIENT DESCENT RESTARTS 发现基本思路是有点类似的)。


    一.两种跳出驻点的方法


    上面说了基本思路,其实思路很朴实也很简单:常见的训练方法在刚开始的时候学习率设置的比较大,这样可以大幅度快步向极值点迈进,加快收敛速度,随着越来越接近极值点,逐步将学习率调小,这样避免步子太大跳过极值点或者造成优化过程震荡结果不稳定。当优化到性能无法再提升的时候,此时大概率陷入了鞍点或者局部极小值,表现为梯度接近于0而且一般此时学习率已经调整到非常小的程度了。跳出驻点的思路就是说:此时可以重新把已经调整的很小的学习率数值放大,强行逼迫优化算法跳出此刻找到的鞍点或者极值点,在上述示意的损失函数曲面上等于当前在某个山谷的谷底(或者某个平台上,即鞍点),此时迈出一大步跳到比较远的地方,然后继续常规的优化过程,就是不断调小学习率,期望找到附近的另外一个山谷谷底,如此往复若干次,然后找出这些谷底哪个谷底更深,然后采取这个最深的谷底对应的参数作为模型参数。


    所以这里的关键是在优化进行不下去的时候重新调大学习率跳出当前的山谷,根据起跳点的不同,可以有两种不同的方法:


      方法1:随机漫步起跳法

    这个方法比较简单,就是说当前优化走到哪个山谷谷底,就从当前位置起跳,直接调大学习率即可,因为这样类似于走到哪里走不下去就从哪里开始跳到附近的远处,走到哪算哪,所以称之为“随机漫步起跳法”。


      方法2:历史最优起跳法

    上面说过,在这种优化方法过程中,会不断进入某个谷底,那么假设在训练集上优化历史上已经经过K个谷底,在验证集上其对应参数的模型性能是可以知道的,分别为P1,P2……Pk。随机漫步起跳法就是从Pk作为起跳点进行起跳。而历史最优起跳法则从P1,P2……Pk中性能最优的那个山谷谷底开始起跳,假设Pi=Max(P1,P2……Pk),那么从第i个山谷开始起跳,如果跳出后发现第j个山谷的性能Pj要优于Pi,则后面从Pj的位置开始起跳,当从Pi跳过K次都没有发现比Pi性能更好的山谷,此时可以终止优化算法。这个有点像是Pi性能的“一山望着一山高”,然后不断从历史最高峰跳到更高的高山上,可以保证性能应该越来越好。从道理上讲感觉这种历史最优起跳法要好于随机漫步起跳法,但是实验结果表明两种不同的起跳方法性能是差不多的,有时候这个稍好点,有时候那个稍好点,没有哪个方法确定性的比另外一个好。


    二.相关的实验及其结果


    为了检验上述跳出驻点的思路是否有效,我们选择了若干深层CNN模型来测试,这里列出其中某个模型的结果,这个模型结构可能看上去有点奇怪,前面说过这个思路是我们去年在做网络最优结构自动探索任务中产生的,所以这两个比较奇怪的结构是从自动生成的网络结构里面随机选出的,性能比较好。分类数据是cifar10图片分类数据,选择带动量的SGD算法作为基准方法来进行相关实验的验证。


         2.1模型结构及参数配置

    CNN模型结构如下图所示:


    其模型结构为:model structure:('init-29', [('conv', [('nb_filter', 16)]), ('conv',[('nb_filter', 256)]), ('maxp', [('pool_size', 2)]), ('conv', [('nb_filter',256)]), ('conv', [('nb_filter', 256)]), ('conv', [('nb_filter', 256)]),('maxp', [('pool_size', 2)]), ('conv', [('nb_filter', 256)])])")


    后面部分接上fc(512,relu)+dropout(0.2)+fc(10,softmax)的分类层;


    学习算法设置的参数为: 选择使用Nesterov的带动量的随机梯度下降法,初始学习率为0.01,衰减因子为1e-6,momentum因子为0.9。

     

    2.2 Nesterov+动量SGD(固定学习率)实验效果

    在上节所述的模型设置下,从结果看,最好的测试精度是0.9060,模型的误差和精度都变化地比较平滑


    2.3学习率自适应减少进行训练的实验结果

    其它配置类似于2.1.2所述,唯一的不同在于:使用常规策略对学习率进行缩小的改变。改变学习率的策略是如果误差10次不变化,那么学习率缩小为原来的0.1倍。 总迭代次数设置为100次。 训练误差、训练准确率、测试误差、测试准确率的变化曲线为:


     

    从训练结果看,第47次学习率变化,性能上有一次大幅度的提高:


    从训练和测试的结果看,测试精度最高的是在第80轮,达到了0.9144。此时的训练误差0.3921,训练准确率0.9926,测试误差0.3957,测试准确率0.9144。


    由此可见,增加学习率的动态改变,在学习速度和质量上要比原始的学习率恒定的训练方法要好,这里性能上增加近一个百分点。

     

    2.4随机漫步起跳法实验结果

     

    从上面的两组实验结果过程的最后阶段可以看出,误差和准确率变化较小,很可能优化过程陷入一个局部最优值,可能是一个极值,也可能是一个鞍点,不管怎样称呼它,它就是一个驻点。(从结果看,应该是一个局部最优值)。在其它实验条件与2.1.3节所述相同情况下,采取上文所述的随机漫步起跳法强迫优化算法跳出当前局部最优,即在学习率比较低的时候,进行扩大学习率的方式(直接设置为0.1,后面再常规的逐步减少),让其在误差曲面的区域中去寻找更多的驻点。 具体做法是:如果当前发现经过三次学习率变小没有带来准确率的提升,那么在当前训练的点,将学习率恢复到初始的学习率即0.01。 迭代的效果如下训练误差、训练准确率、测试误差、测试准确率的变化曲线如下图为:  


    从训练过程看,从当前迭代的当前位置出发,采用学习率缩小和扩大相结合的方式,目前最高能在第153轮的准确率提高到0.9236,比之前最好的0.9144有一点提升,说明这样的学习率变大策略也是一种可行的思路。


    另外,从图中看,由于我们变大学习率许多次,每次学习率变大都会呈现一个波动形式,故结果表现出锯齿状的波动趋势。


    2.5历史最优起跳法实验结果

       我们做了类似的实验,发现历史最优起跳法与随机漫步起跳法性能相当,有时这个好点有时那个好点,没有发现性能明显差异。

     

    我们随机抽了其它几个最优网络结构自动探索出的深度神经网络结构,表现与上述网络结构类似,而如果基础网络结构本身性能越差,这个差别越大。


    所以可以得出的结论是:优化过程中采用常规动态学习率调整效果要优于固定学习率的方法,而随机漫步起跳法效果略优于常规的动态学习率调整方法,但是历史最优起跳法相对常规动态学习率方法没有发现改善。

     

    相关文献

    1. SGDR: STOCHASTIC GRADIENT DESCENT RESTARTS .ICLR-2017
    2. SNAPSHOT ENSEMBLES : TRAIN 1,GET M FOR FREE

    展开全文
  • 《图像处理实例》 之 局部极值提取

    千次阅读 2017-08-23 02:38:00
    局部极值提取算法 这是http://www.imagepy.org/的作者原创,我只是对其理解之后改进说明,欢迎大家使用这个小软件! 有需要C++的朋友,可有留邮箱,经测试改进的版本! 算法原理:(按照程序思路来,...

     


     

    局部极值提取算法

    这是http://www.imagepy.org/的作者原创,我只是对其理解之后改进和说明,欢迎大家使用这个小软件!

    有需要C++的朋友,可有留邮箱,经测试和改进的版本! 

     


     

    算法原理:(按照程序思路来,当然其中很多可以改进的地方)

     

        第一步:    

            先进行距离变换(这类问题都是要通过几何操作进行,距离变换就是几何距离)

     

        第二步:

            先找全局可能最大值,然后对图像进行标记

            全局可能最大值是基础,3X3核在图像上滑动去找全局可能最大值。标记是为了掩膜操作更方便,你看opencv很多函数都有掩膜,PS上面也有那个白色的膜。     

            图片结果是:背景和最大值标记3,前景标记1,图像边界标记0(为了防止越界)
           貌似下面的图片显示有问题,matploylib的显示有时候确实有问题,Debug的时候就和显示不一样,没关系这一步比较简单。

     

        第三步:

            对全局可能最大值进行排除,这一步是程序核心!这里以求最大值为例,最小值相反。

            •   首先对找到的最大值进行排序--从大到小
            •   按照从大到小的顺序去寻找每个最大值的领地
            •   领地按照“目标范围”和“目标深度”两个原则,具体表现形式为:“不重合”、“不包括”、“领主大于等于领地”             

              下面以图说明两个大原则:

                目标范围:

                    下图A1这个最大最他所占领的范围由用户输入,假设半径为y,则 x1 = x2 = y ,x1和x2的范围就是A1的范围,P点虽然很低,但是不属于A1的范围,所以不去占领。

     

                     目标深度:

                     
                    下图A1的深度也又用户输入,假设输入为y,则 y1 = y,高度确定之后,那么P点虽然也在山脊上,但是不会被A1所占领。

                    

     

              下面以图为例进一步说明三个具体形式:

              注释:三个具体形式是建立在两个原则之上的,所以理解上面是结局算法的基础。

                不重合形式:

                    下图A1和A2有领土纠纷,但是A1峰值高所以先进行领地划分,当A2进行领地划分的时候一旦接触A1的领地就立马被消灭,最后只剩A1的存在。

                不包括形式:

                           下图A1先进行领土划分,过程中直接把A2消灭了,不存在A2的领土划分

                领主大于等于领地形式,也可以叫同等峰值先来先得形式:

                      由于A1先进行领土划分,所以后面的都被消灭了,读到这里你会发现这个程序的BUG,后面再详细谈论这一点

     

     

                     目标深度原则

                    要确定和你预定的深度是不是符合,比如你想要山峰之间的深度都是在10cm以上,那么有些不符合的就得去除,去除原则肯定矮的GG

    注释:大佬的程序没有几何距离,还有部分Bug,部分经过修正了。

     

     

     

     

     程序有待改进:

             这部分可以通过处理---“同样的峰值且山脉相连”的一类极值点,对其平均、覆盖范围最多等操作取一个中间的值!

             由于python用的不熟练,而且现在时间太晚了脑子转不动,以后用C++写的时候再改进这一点。

     

     

    2017.10.17更新:

     

          1.对上面第二步:可能最大值进行补充

            请看下面这个图,中间位置肯定不是最大值区域,而且通过上面的第三步过滤也是没办法过滤的,具体看上面原理就知道了,所以需要在第二步进行过滤!

            本程序利用如下方法:

              A.每一层的最大值都大于外面一层,sum[0]>=sum[1] && sum[1]>sum[2] && sum[2]>sum[3],这里使用到9X9的内核

              B.记录每一层最大值出现的位置X0和X1,然后最大值的距离(abs(x0.x - x1.x) <= 2 && abs(x0.y - x1.y) <= 2)小于等于2,这里不严谨,因为每层最大值不止一个。

              C.记录最大值出现的数量n_count >= 3

     

     1        float sum[4] = { 0 };//存储最大值进行对比、在这里说不清楚看博客2017.10.17更新
     2             sum[0] = img[3][j];
     3             Point x0 = Point(0, 0);
     4             Point x1 = Point(0, 0);
     5             uchar n_count = 0;
     6             for (int m = 2; m < 5; m++)
     7             {
     8                 for (int n = -1; n < 2; n++)
     9                 {
    10                     if (m == 3 && n == 0) continue;
    11                     sum[1] = sum[1] < img[m][j + n] ? img[m][j + n] : sum[1];
    12                     x0 = sum[0] == img[m][j + n] ? Point(m, n) : x0;
    13                     n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
    14                     //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
    15                 }
    16             }
    17             for (int m = 1; m < 6; m++)
    18             {
    19                 for (int n = -2; n < 3; n++)
    20                 {
    21                     if (2 <= m && m <= 4 && -1 <= n && n <= 1) continue;
    22                     sum[2] = sum[2] < img[m][j + n] ? img[m][j + n] : sum[2];
    23                     x1 = sum[0] == img[m][j + n] ? Point(m, n) : x1;
    24                     n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
    25                     //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
    26                 }
    27             }
    28             for (int m = 0; m < 7; m++)
    29             {
    30                 for (int n = -3; n < 4; n++)
    31                 {
    32                     sum[3] = sum[3] < img[m][j + n] && !(1 <= m && m <= 5 && -2 <=n && n <= 2) ? img[m][j + n] : sum[3];
    33                     //flag = img[3][j+0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
    34                 }
    35             }
    36             x1 = (x1.x == 0 && x1.y == 0) || n_count >= 3 ? x0 : x1;//判断是否存在5X5的最大值(和目标点相同)
    37             int tmp = sum[0] >= sum[1] && sum[1] >= sum[2] && sum[2] >= sum[3] && (abs(x0.x - x1.x) <= 2 && abs(x0.y - x1.y) <= 2)
    38                 ? 0 : FillBlock(src, mask_tmp, Point(j, i));//tmp没意义,就是为了调用后面的函数而已
    39         }

     

     

     

           2.对第三步进行补充

            A.搜索过程遇到边界,那就把这个最大值点去除if (img[i][fill_point[count].x + j] == 2 || msk[i][fill_point[count].x + j] == 0) max_point[k].data = 1;

          3.效果图

            注释:满足一个条件就判定为最大值!

               1 Find_Max(img, mask,0, 5000); 

                   1 Find_Max(img, mask,5000, 0); 

                   1 Find_Max(img, mask,5000, 5000); 

                   1 Find_Max(img, mask,10, 20); 

                   1 Find_Max(img, mask,10, 50); 

     

     

      1 import scipy.ndimage as ndimg
      2 import numpy as np
      3 from numba import jit
      4 import cv2
      5 
      6 def neighbors(shape):
      7     dim = len(shape)
      8     block = np.ones([3] * dim)
      9     block[tuple([1] * dim)] = 0
     10     idx = np.where(block > 0)
     11     idx = np.array(idx, dtype=np.uint8).T
     12     idx = np.array(idx - [1] * dim)
     13     acc = np.cumprod((1,) + shape[::-1][:-1])
     14     return np.dot(idx, acc[::-1])
     15 
     16 
     17 @jit  # trans index to r, c...
     18 
     19 def idx2rc(idx, acc):
     20     rst = np.zeros((len(idx), len(acc)), dtype=np.int16)
     21     for i in range(len(idx)):
     22         for j in range(len(acc)):
     23             rst[i, j] = idx[i] // acc[j]
     24             idx[i] -= rst[i, j] * acc[j]
     25     return rst
     26 
     27 
     28 #@jit  # fill a node (may be two or more points)
     29 
     30 def fill(img, msk, p, nbs, buf):
     31     msk[p] = 3
     32     buf[0] = p
     33     back = img[p]
     34     cur = 0
     35     s = 1
     36     while cur < s:
     37         p = buf[cur]
     38         for dp in nbs:
     39             cp = p + dp
     40             if img[cp] == back and msk[cp] == 1:
     41                 msk[cp] = 3
     42                 buf[s] = cp
     43                 s += 1
     44                 if s == len(buf):
     45                     buf[:s - cur] = buf[cur:]
     46                     s -= cur
     47                     cur = 0
     48         cur += 1
     49     #msk[p] = 3
     50 
     51 
     52 #@jit  # my mark
     53 
     54 def mark(img, msk, buf, mode):  # mark the array use (0, 1, 2)
     55     omark = msk     
     56     nbs = neighbors(img.shape)
     57     idx = np.zeros(1024 * 128, dtype=np.int64)
     58     img = img.ravel()  # 降维
     59     msk = msk.ravel()  # 降维
     60     s = 0
     61     for p in range(len(img)):
     62         if msk[p] != 1: continue  
     63         flag = False             
     64         for dp in nbs:
     65             if mode and img[p + dp] > img[p]: 
     66                 flag = True
     67                 break
     68             elif not mode and img[p + dp] < img[p]:
     69                 flag = True
     70                 break
     71         
     72         if flag : continue
     73         else    : fill(img, msk, p, nbs, buf)
     74         idx[s] = p
     75         s += 1
     76         if s == len(idx): break
     77     plt.imshow(omark, cmap='gray')
     78     return idx[:s].copy()
     79 
     80 
     81 
     82 def filter(img, msk, idx, bur, tor, mode):
     83     omark = msk  
     84     nbs = neighbors(img.shape)
     85     acc = np.cumprod((1,) + img.shape[::-1][:-1])[::-1]
     86     img = img.ravel()
     87     msk = msk.ravel()
     88 
     89     arg = np.argsort(img[idx])[::-1 if mode else 1] 
     90 
     91     for i in arg:
     92         if msk[idx[i]] != 3:    
     93             idx[i] = 0
     94             continue
     95         cur = 0
     96         s = 1
     97         bur[0] = idx[i] 
     98         while cur < s:
     99             p = bur[cur]
    100             if msk[p] == 2:     
    101                 idx[i] = 0
    102                 break
    103 
    104             for dp in nbs:
    105                 cp = p + dp
    106                 if msk[cp] == 0 or cp == idx[i] or msk[cp] == 4: continue
    107                 if mode and img[cp] < img[idx[i]] - tor: continue
    108                 if not mode and img[cp] > img[idx[i]] + tor: continue
    109                 bur[s] = cp
    110                 s += 1
    111                 if s == 1024 * 128:
    112                     cut = cur // 2
    113                     msk[bur[:cut]] = 2
    114                     bur[:s - cut] = bur[cut:]
    115                     cur -= cut
    116                     s -= cut
    117 
    118                 if msk[cp] != 2: msk[cp] = 4    
    119             cur += 1
    120         msk[bur[:s]] = 2    
    121         #plt.imshow(omark, cmap='gray')
    122 
    123     return idx2rc(idx[idx > 0], acc)
    124 
    125 
    126 def find_maximum(img, tor, mode=True):
    127     msk = np.zeros_like(img, dtype=np.uint8)  
    128     msk[tuple([slice(1, -1)] * img.ndim)] = 1  
    129     buf = np.zeros(1024 * 128, dtype=np.int64)
    130     omark = msk
    131     idx = mark(img, msk, buf, mode)
    132     plt.imshow(msk, cmap='gray')
    133     idx = filter(img, msk, idx, buf, tor, mode)
    134     return idx
    135 
    136 
    137 if __name__ == '__main__':
    138     from scipy.misc import imread
    139     from scipy.ndimage import gaussian_filter
    140     from time import time
    141     import matplotlib.pyplot as plt
    142 
    143     img = cv2.imread('123.png')
    144     img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    145     ret2, img = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)
    146     img[:] = ndimg.distance_transform_edt(img)
    147     plt.imshow(img, cmap='gray')
    148     pts = find_maximum(img, 20, True)
    149     start = time()
    150     pts = find_maximum(img, 10, True)
    151     print(time() - start)
    152     plt.imshow(img, cmap='gray')
    153     plt.plot(pts[:, 1], pts[:, 0], 'y.')
    154     plt.show()

     

     C++版本

     老版本、不稳定,可以看思路

      1 //---fill black value
      2 int FillBlock(Mat src, Mat &mask, Point center)
      3 {
      4     uchar back = src.at<uchar>(center.y, center.x);
      5     vector<Point> fill_point;
      6     int count = 0, count_mount = 1;
      7     fill_point.push_back(center);
      8     while (count < count_mount)
      9     {
     10         vector<uchar*> img;
     11         vector<uchar*> msk;
     12         for (int i = -1; i < 2; i++)
     13         {
     14             img.push_back(src.ptr<uchar>(fill_point[count].y + i));
     15             msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
     16         }
     17         for (size_t i = 0; i < 3; i++)
     18         {
     19             for (int j = -1; j < 2; j++)
     20             {
     21                 if (img[i][fill_point[count].x + j] == back && !(j == 0 && i == 1) && msk[i][fill_point[count].x + j] == 255)
     22                 {
     23                     fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
     24                     msk[i][fill_point[count].x + j] = 1;
     25                 }
     26             }
     27         }
     28         msk[1][fill_point[count].x] = 1;
     29         count_mount = fill_point.size() - 1;
     30         fill_point.erase(fill_point.begin());
     31     }
     32     return 0;
     33 }
     34 //---cal mask
     35 //---@_src        
     36 //---@mask        
     37 void MaskImage(InputArray _src, Mat &mask)
     38 {
     39     Mat src = _src.getMat(),mask_tmp = Mat::zeros(src.size(), CV_8UC1);
     40     mask_tmp.setTo(255);
     41     Mat rows = Mat::zeros(Size(src.cols, 1), CV_8UC1), cols = Mat::zeros(Size(1, src.rows), CV_8UC1);
     42     Mat src_rows_beg = mask_tmp.row(0);
     43     Mat src_rows_end = mask_tmp.row(src.rows - 1);
     44     Mat src_cols_beg = mask_tmp.col(0);
     45     Mat src_cols_end = mask_tmp.col(src.cols - 1);
     46     rows.copyTo(src_rows_beg); rows.copyTo(src_rows_end);
     47     cols.copyTo(src_cols_beg); cols.copyTo(src_cols_end);
     48     for (size_t i = 1; i < src.rows-1; i++)
     49     {
     50         uchar *img0 = src.ptr<uchar>(i - 1);
     51         uchar *img  = src.ptr<uchar>(i);
     52         uchar *img1 = src.ptr<uchar>(i + 1);
     53         uchar *msk  = mask_tmp.ptr<uchar>(i);    
     54         for (size_t j = 1; j < src.cols-1; j++)
     55         {
     56             bool flag = false;
     57             //msk[j] = img[j] == 0 ? 0 : msk[j];
     58             if (msk[j] != 255) continue;
     59             flag = (img[j] < img[j - 1] || img[j] < img[j + 1]
     60                 || img[j] < img0[j] || img[j] < img0[j - 1]
     61                 || img[j] < img0[j + 1] || img[j] < img1[j]
     62                 || img[j] < img1[j - 1] || img[j] < img1[j + 1])
     63                 ? true : false;
     64             int tmp = flag == true ? FillBlock(src, mask_tmp, Point(j, i)) : 0;
     65         }
     66     }
     67     mask = mask_tmp.clone();
     68 }
     69 //---filter parts max value
     70 //---@
     71 //---@
     72 //---@gap        
     73 //---@radius    
     74 
     75 vector<Point> Find_Max(InputArray _src, Mat&mask,int gap,int radius)
     76 {
     77     Mat src = _src.getMat();
     78     
     79     typedef struct MyStruct
     80     {
     81         Point position;
     82         float data;
     83     }MyStruct;
     84 
     85     MaskImage(src, mask);
     86     vector<MyStruct> max_point;
     87     for (size_t i = 0; i < src.rows; i++)
     88     {
     89         uchar *img = src.ptr<uchar>(i);
     90         uchar *msk = mask.ptr<uchar>(i);
     91         for (size_t j = 0; j < src.cols; j++)
     92         {
     93             if (msk[j] != 255) continue;
     94             MyStruct my_struct;
     95             my_struct.data = img[j];
     96             my_struct.position = Point(j, i);
     97             max_point.push_back(my_struct);
     98         }
     99     }
    100     for (size_t i = 0; i < max_point.size(); i++)
    101     {
    102         for (size_t j = i; j < max_point.size(); j++)
    103         {
    104             MyStruct temp;
    105             if (max_point[i].data <= max_point[j].data)
    106             {
    107                 if (max_point[j].data == 0) continue;
    108                 temp = max_point[j];
    109                 max_point[j] = max_point[i];
    110                 max_point[i] = temp;
    111             }
    112         }
    113     }
    114     //---find max
    115     for (size_t k = 0; k < max_point.size(); k++)//---
    116     {
    117         uchar back = src.at<uchar>(max_point[k].position.y, max_point[k].position.x);
    118         vector<Point> fill_point;
    119         int count = 0, count_mount = 1;
    120         fill_point.push_back(max_point[k].position);
    121         
    122         while (count < count_mount &&  max_point[k].data != 1)
    123         {
    124             vector<uchar*> img;
    125             vector<uchar*> msk;
    126             for (int i = -1; i < 2; i++)
    127             {
    128                 img.push_back(src.ptr<uchar>(fill_point[count].y + i));
    129                 msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
    130             }
    131             for (int i = 0; i < 3; i++)
    132             {
    133                 for (int j = -1; j < 2; j++)
    134                 {
    135                     //---
    136                     uchar x = pow((max_point[k].position.x - fill_point[count].x + j), 2); //(max_point[k].position.x - img[i][fill_point[count].x + j])*(max_point[k].position.x - img[i][fill_point[count].x + j]);
    137                     uchar y = pow((max_point[k].position.y - (fill_point[count].y + i - 1)) , 2); // (max_point[k].position.y - img[i][fill_point[count].y + j])*(max_point[k].position.y - img[i][fill_point[count].x + j]);
    138                     uchar distance = sqrt(x + y);
    139                     if (img[i][fill_point[count].x + j] <= img[1][fill_point[count].x] - gap 
    140                         || msk[i][fill_point[count].x + j] == 3
    141                         || msk[i][fill_point[count].x + j] == 0
    142                         || (j == 0 && i == 1)
    143                         || distance >= radius) continue;
    144                     if (img[i][fill_point[count].x + j] == 2) max_point[k].data = 1;
    145                     msk[i][fill_point[count].x + j] = 3;
    146                     fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
    147                     count_mount++;
    148                 }
    149             }
    150             count++;
    151         }    
    152         if (max_point[k].data == 1)
    153         {
    154             for (size_t i = 0; i < fill_point.size(); i++)
    155             {
    156                 mask.at<uchar>(fill_point[i]) = 1;
    157             }
    158         }
    159         else
    160         {
    161             for (size_t i = 0; i < fill_point.size(); i++)
    162             {
    163                 mask.at<uchar>(fill_point[i]) = 2;
    164             }
    165             max_point[k].data = 255;
    166             mask.at<uchar>(max_point[k].position) = 255;
    167         }
    168     }
    169 }
    View Code

     

     

    C++版本:

    2017.10.17更新

      1 int FillBlock(Mat src, Mat &mask, Point center)
      2 {
      3     uchar back = src.at<uchar>(center.y, center.x);
      4     vector<Point> fill_point;
      5     int count = 0, count_mount = 1;
      6     fill_point.push_back(center);
      7     while (count < count_mount)
      8     {
      9         vector<uchar*> img;
     10         vector<uchar*> msk;
     11         for (int i = -1; i < 2; i++)
     12         {
     13             img.push_back(src.ptr<uchar>(fill_point[count].y + i));
     14             msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
     15         }
     16         for (size_t i = 0; i < 3; i++)
     17         {
     18             for (int j = -1; j < 2; j++)
     19             {
     20                 if (img[i][fill_point[count].x + j] == back && !(j == 0 && i == 1) && msk[i][fill_point[count].x + j] == 255)
     21                 {
     22                     fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
     23                     msk[i][fill_point[count].x + j] = 1;
     24                 }
     25             }
     26         }
     27         msk[1][fill_point[count].x] = 1;
     28         count_mount = fill_point.size() - 1;
     29         fill_point.erase(fill_point.begin());
     30     }
     31     return 0;
     32 }
     33 
     34 void MaskImage(InputArray _src, Mat &mask)
     35 {
     36     Mat src = _src.getMat(),mask_tmp = Mat::zeros(src.size(), CV_8UC1);
     37     mask_tmp.setTo(255);
     38     Mat rows = Mat::zeros(Size(src.cols, 1), CV_8UC1), cols = Mat::zeros(Size(1, src.rows), CV_8UC1);
     39     Mat src_rows_beg = mask_tmp.row(0);
     40     Mat src_rows_end = mask_tmp.row(src.rows - 1);
     41     Mat src_cols_beg = mask_tmp.col(0);
     42     Mat src_cols_end = mask_tmp.col(src.cols - 1);
     43     rows.copyTo(src_rows_beg); rows.copyTo(src_rows_end);
     44     cols.copyTo(src_cols_beg); cols.copyTo(src_cols_end);
     45     
     46     for (size_t i = 3; i < src.rows-3; i++)
     47     {
     48         vector<uchar*> img;
     49         uchar* msk = mask_tmp.ptr(i);
     50         uchar* img1 = src.ptr(i);
     51         for (int k = -3; k < 4; k++)
     52         {
     53             img.push_back(src.ptr<uchar>(k + i));
     54         }
     55         for (size_t j = 3; j < src.cols-3; j++)
     56         {
     57             bool flag = false;
     58             if (msk[j] != 255) continue;
     59             float sum[4] = { 0 };
     60             sum[0] = img[3][j];
     61             Point x0 = Point(0, 0);
     62             Point x1 = Point(0, 0);
     63             uchar n_count = 0;
     64             for (int m = 2; m < 5; m++)
     65             {
     66                 for (int n = -1; n < 2; n++)
     67                 {
     68                     if (m == 3 && n == 0) continue;
     69                     sum[1] = sum[1] < img[m][j + n] ? img[m][j + n] : sum[1];
     70                     x0 = sum[0] == img[m][j + n] ? Point(m, n) : x0;
     71                     n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
     72                     //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
     73                 }
     74             }
     75             for (int m = 1; m < 6; m++)
     76             {
     77                 for (int n = -2; n < 3; n++)
     78                 {
     79                     if (2 <= m && m <= 4 && -1 <= n && n <= 1) continue;
     80                     sum[2] = sum[2] < img[m][j + n] ? img[m][j + n] : sum[2];
     81                     x1 = sum[0] == img[m][j + n] ? Point(m, n) : x1;
     82                     n_count = sum[0] == img[m][j + n] ? n_count+1 : n_count;
     83                     //flag = img[3][j + 0] < img[m][j + n] ? true : flag;//如果目标像素小于周围任何一个像素就说明这个一定不是最大值
     84                 }
     85             }
     86             for (int m = 0; m < 7; m++)
     87             {
     88                 for (int n = -3; n < 4; n++)
     89                 {
     90                     sum[3] = sum[3] < img[m][j + n] && !(1 <= m && m <= 5 && -2 <=n && n <= 2) ? img[m][j + n] : sum[3];
     91                     //flag = img[3][j+0] < img[m][j + n] ? true : flag;
     92                 }
     93             }
     94             x1 = (x1.x == 0 && x1.y == 0) || n_count >= 3 ? x0 : x1;
     95             int tmp = sum[0] >= sum[1] && sum[1] >= sum[2] && sum[2] >= sum[3] && (abs(x0.x - x1.x) <= 2 && abs(x0.y - x1.y) <= 2)
     96                 ? 0 : FillBlock(src, mask_tmp, Point(j, i));
     97         }
     98     }
     99     mask = mask_tmp.clone();
    100 }
    101 
    102 vector<Point> Find_Max(InputArray _src, Mat&mask,int gap,int radius)
    103 {
    104     Mat src = _src.getMat();
    105     
    106     typedef struct MyStruct
    107     {
    108         Point position;
    109         float data;
    110     }MyStruct;
    111 
    112     MaskImage(src, mask);
    113     vector<MyStruct> max_point;
    114     for (size_t i = 0; i < src.rows; i++)
    115     {
    116         uchar *img = src.ptr<uchar>(i);
    117         uchar *msk = mask.ptr<uchar>(i);
    118         for (size_t j = 0; j < src.cols; j++)
    119         {
    120             if (msk[j] != 255) continue;
    121             MyStruct my_struct;
    122             my_struct.data = img[j];
    123             my_struct.position = Point(j, i);
    124             max_point.push_back(my_struct);
    125         }
    126     }
    127     for (size_t i = 0; i < max_point.size(); i++)
    128     {
    129         for (size_t j = i; j < max_point.size(); j++)
    130         {
    131             MyStruct temp;
    132             if (max_point[i].data <= max_point[j].data)
    133             {
    134                 if (max_point[j].data == 0) continue;
    135                 temp = max_point[j];
    136                 max_point[j] = max_point[i];
    137                 max_point[i] = temp;
    138             }
    139         }
    140     }
    141 
    142     for (size_t k = 0; k < max_point.size(); k++)//---
    143     {
    144         uchar back = src.at<uchar>(max_point[k].position.y, max_point[k].position.x);
    145         vector<Point> fill_point;
    146         int count = 0, count_mount = 1;
    147         fill_point.push_back(max_point[k].position);
    148         
    149         while (count < count_mount &&  max_point[k].data != 1)
    150         {
    151             vector<uchar*> img;
    152             vector<uchar*> msk;
    153             for (int i = -1; i < 2; i++)
    154             {
    155                 img.push_back(src.ptr<uchar>(fill_point[count].y + i));
    156                 msk.push_back(mask.ptr<uchar>(fill_point[count].y + i));
    157             }
    158             for (int i = 0; i < 3; i++)
    159             {
    160                 for (int j = -1; j < 2; j++)
    161                 {
    162                     //---
    163                     double x = pow((max_point[k].position.x - fill_point[count].x + j), 2); //(max_point[k].position.x - img[i][fill_point[count].x + j])*(max_point[k].position.x - img[i][fill_point[count].x + j]);
    164                     double y = pow((max_point[k].position.y - (fill_point[count].y + i - 1)), 2); // (max_point[k].position.y - img[i][fill_point[count].y + j])*(max_point[k].position.y - img[i][fill_point[count].x + j]);
    165                     int distance = sqrt(x + y);
    166                     if (img[i][fill_point[count].x + j] <= img[0][fill_point[count].x] - gap 
    167                         || msk[i][fill_point[count].x + j] == 3
    168                         //|| msk[i][fill_point[count].x + j] == 0 
    169                         || (j == 0 && i == 1)
    170                         || distance >= radius) continue;
    171                     if (img[i][fill_point[count].x + j] == 2 || msk[i][fill_point[count].x + j] == 0) max_point[k].data = 1;
    172                     msk[i][fill_point[count].x + j] = 3;
    173                     fill_point.push_back(Point(fill_point[count].x + j, fill_point[count].y + i - 1));
    174                     count_mount++;
    175                 }
    176             }
    177             count++;
    178         }    
    179         if (max_point[k].data == 1)
    180         {
    181             for (size_t i = 0; i < fill_point.size(); i++)
    182             {
    183                 mask.at<uchar>(fill_point[i]) = 1;
    184             }
    185         }
    186         else
    187         {
    188             for (size_t i = 0; i < fill_point.size(); i++)
    189             {
    190                 mask.at<uchar>(fill_point[i]) = 2;
    191             }
    192             max_point[k].data = 255;
    193             mask.at<uchar>(max_point[k].position) = 255;
    194         }
    195     }
    196     vector<Point> print_wjy;
    197     for (size_t i = 0; i < mask.rows; i++)
    198     {
    199         uchar *msk = mask.ptr<uchar>(i);
    200         for (size_t j = 0; j < mask.cols; j++)
    201         {
    202             if (msk[j] == 255)
    203                 print_wjy.push_back(Point(j, i));
    204         }
    205 
    206     }
    207     return print_wjy;
    208 }

     

    转载于:https://www.cnblogs.com/wjy-lulu/p/7416135.html

    展开全文
  • Given symmetric matrices B,D ∈ R n×n and a symmetric positive definite matrix W ∈ R n×n , maximizingthe sum of the Rayleighquotientx ? Dx andthe gener- alized Rayleigh quotient ...
  • 皇后问题局部极值启发式搜索方法 背景说明:皇后问题是算法领域的著名问题,本篇主要是对皇后问题背景下的局部搜索问题的研究。

    8-Queens Problem皇后问题局部极值启发式搜索方法


    Backgrounds

    背景说明:

    皇后问题是算法领域的著名问题,问题的背景是在8*8的国际象棋棋盘上摆满8个皇后使之互相不能攻击(由于皇后在国际象棋规则中横、纵、斜三个方向均可以移动,即摆放要使得皇后两两不在同一横、纵、斜线上)。放松这个问题的变长条件,我们可以将其扩大至N*N的棋盘上摆放N个皇后,使之各自在三个方向上不重复出现。

    大一时刚刚学习编程语言的我们,在初次接到这道题时,主要是运用了递归函数和回溯的思想去解决问题,通过深度优先搜索(depth-first search, DFS)策略,依次增加棋盘上可以摆放的位置。若遇到第k步棋子摆放的所有可能位置都发生冲突,则退回到第k-1步,修改k-1的状态,以达到期望的最终解。


    State space of 8-queens problem

    皇后问题的状态空间:

    状态空间,这里理解成所有摆放的可能性情况。

    皇后问题广义的解的状态空间,在没有过多限制条件的情况下是非常庞大的。

    在64个互不重复的位置上选择无差别的8个位置,总计应该有C(64, 8)=4,426,165,368 (44亿)大小的状态空间. 图示为较为集中的分布情况和可以完全分散开的分布情况。

      

            图 1.1 皇后分布紧凑          图 1.2 皇后分布松散 

    注意这里C(N*N, N)表示在N*N种方案中选择N个的方案数

    不同约束条件下的不同状态空间规模:

    实际上,全局解状态空间的大小也由我们的限制条件决定,下面提出两种解决方案。但总体上策略一致地是,我们用一个长度为8的一维数组表示每一列皇后的位置,这种表示方法首先就压缩掉了很大一部分的解的状态空间——因为相当于即使在横向上有皇后可能会重合,这种每列一个位置数的表示方式放弃了同一列上出现多个皇后的情况。图示为横纵均不重合分散分布和允许行重复的分布情况。  

      

           图 2.1 允许同行冲突 图 2.2 保证横纵均不冲突


    但上述两种表示方法出于求解出尽可能多的局部极值和尽快求解出最优解的考虑,依然是稳妥的。此时状态空间的大小至多已经被压缩至8^8=16,777,216 (1670万)种,相比于44亿已经极大地压缩了解的数量。

    (1).第一种表示方法,我们在允许行重合的情况下,每次操作只在该列上对棋子进行上下移动,并考虑每次移动对于总冲突数的影响。

    此种排列的状态空间为8^8=16,777,216.

    (2).第二种表示方法,我们对1~8的数字进行全排列,得到{a1, a2, …, ai, …, a8}用于表示第i 列上的皇后位置,这样可以保证行列均不重复。

    此种排列的状态空间为8!=40,320. 


    局部极值和全局最值的理解:

    全局最值,是对于每一个搜索问题需要找到的最终结果。对于不同的问题,其问题的离散/连续性可能不同,解的数量、状态分布也不相同。皇后问题本身是一个离散问题(解的状态空间数量是有限个数的,虽然解空间总规模可能非常大),且解的数量也是巨大的。至少对于8皇后,不存在唯一的状态。后续运行程序可知,在不考虑翻转、旋转等价的情况下,有92组互不相同的解。即在横纵均不重合的40,320规模的状态空间(或在更大范围的状态空间中)有92组全局最值。

    局部极值,是我们通过某种途径趋向全局最值时经过的、不能再用当前算法得出相对更进一步的结果(而必须由重启、重新打乱而继续开始一次新的搜索)时,“卡住”的状态。


    图 3 搜索路径与全局最优、局部最优的关系

    (state space为构成搜索路径的n维解空间的子集,objective function为目标函数)

    对于同一个问题,我们考虑f(x1, x2, …, xi, …, xn)=K的目标函数,希望通过启发式搜索的方式寻找全局K的最优值,即在n维解空间上,我们得到由目标函数f(.)决定的决策-目标曲面。在这个n+1维曲面S(X<x1, x2, …, xn>, K)上,不同的约束条件和算法就相当于通过不同的n维路径Ls<x1s, x2s, …, xns>去逼近全局最优K.

    其中Ls是某种启发式搜索算法的搜索路径,对于特定的搜索算法,Ls是n维解状态空间的有向路径。记Ls上的每一步解为Xj, 并设该条路径从某一初始状态X0出发,经过p步可以到达全局或局部最优值(其中任意Xj都是n维状态空间上的点),则该路径Ls可以表示为Ls={X0, X1, …, Xj, …, Xp}. 由于状态空间的不平整性,显然通过不同的算法定义到的局部最优解(局部极值)也是不同的。

    注意这里的定义包括后续目标函数的定义、估值函数的定义、最终的求解目的等等


    Evaluation function

    估值函数:

    对于启发式搜索,我们需要在巨大的解的状态空间中尽可能快地找出尽可能多的全局最优解,这就需要我们在做每一步选择(选择继续扩大解的已知范围,使之尽快铺展至n维或跳转至下一个解的暂留状态)时,能有效评判该解对于“发现最优解”这一目的的贡献意义的函数。或者我们也可以这样理解:搜索算法在在搜索路径的每一个解Xj停留时,并不知道该解是否处于到达全局最优的必经路径上,因此搜算算法会根据一个估值函数,对“该步操作可以到达全局最优或局部最优”这一事件做出概率分析。当搜索路径Ls走到Xj-1这一步时,会对所有可能的后续状态集合{Xj’|Xj’=Xj’1, Xj’2,…}进行估值,最后选出估值最大的一个Xj’作为Xj下一步。后续也据此方法继续走下去。

    当最终走到某一步,估值函数对所有可能的后继都没有最优估值或正向估值时,就说明搜索路径已经到达极值。此时检查目标函数K即可判断当前解是否是全局最优。

    需要说明的是,对于全局最优,我们也可以有不同的定义方法。在本问题中,我们选择用“冲突行、列或斜线上棋子数减一”表示冲突情况,即当某行、列或斜线上皇后棋子数目为0个或1个时,认为没有冲突;皇后棋子数目为2个及以上时,认为冲突数是皇后数目-1. 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    int  getChange( int  i,  int  j) {
         int  pi = position[i], pj = position[j];
         // i-pi, j-pj
         // left: i+pi==j+pj, right: i-pi==j-pj
         // left[i+pi-1] left[j+pj-1] right[N-(i-pi)] right[N-(j-pj)]
         if  (i+pi == j+pj) {
             // same left line
             return  ((left[i+pi-1]>=3)+1)
             + (right[N-(i-pi)]>=2)
             + (right[N-(j-pj)]>=2)
             - ((right[N-(i-pj)]>=1)+1)
             - (left[i+pj-1]>=1)
             - (left[j+pi-1]>=1);
         else  if  (i-pi == j-pj) {
             // same right line
             return  ((right[N-(i-pi)]>=3)+1)
             + (left[i+pi-1]>=2)
             + (left[j+pj-1]>=2)
             - ((left[i+pj-1]>=1)+1)
             - (right[N-(i-pj)]>=1)
             - (right[N-(j-pi)]>=1);
         else  {
             // (i, j) no collision
             return  (left[i+pi-1]>=2)
             + (left[j+pj-1]>=2)
             + (right[N-(i-pi)]>=2)
             + (right[N-(j-pj)]>=2)
             - (left[i+pj-1]>=1)
             - (left[j+pi-1]>=1)
             - (right[N-(i-pj)]>=1)
             - (right[N-(j-pi)]>=1);
         }
    }
        

    图 4 估值函数

    当发现两列(i, j) 可以被操作时,我们通过计算“交换后”与“交换前”两种状态(两个不同的解的停留点)之间目标函数即总冲突数的变化量,判断该步操作是否应被执行。其中我们考虑了交换前的(i, j)两列处于左斜向冲突(交换后存在右斜向冲突)、右斜向冲突(交换后存在左斜向冲突)、无冲突(交换后彼此依然无冲突)三种情况。 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    void  statistic() {
         clear();
         for  ( int  i=1; i<=N; ++i) {
             int  pi=position[i];
             left[i+pi-1]++;
             right[N-(i-pi)]++;
         }
    }
         
    void  reCount() {
         statistic();
         Collision = 0;
    leftCollision = 0, rightCollision = 0;
         for  ( int  i=1; i<=2*N-1; ++i) {
             if  (left[i]>=2) leftCollision += left[i]-1;
             if  (right[i]>=2) rightCollision += right[i]-1;
         }
         Collision = leftCollision + rightCollision;
    }

    图 5 目标函数定义,显然Collision=0时取得全局最优解


    Main Algorithm

    状态空间压缩:

    我们选择交换不重复的列来完成全局最优的搜索结果,这样可以在本身就比较小的状态空间(40,320种)范围内得出结果。并且对于全局最优来说,由于所有的全局最优解均包含在该状态空间中,全局最优的分布率也最高。

    重启策略:

    初始化时生成1到8的升序列并进行时间复杂度O(n)的遍历,每次随机出一个位置,将其与当前扫到的位置做一次交换。之后每次重启仅对现有列进行一次随机交换即可得到一个新的解空间上的X0作为新的一条搜索路径Ls的起点。

    搜索策略: 

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    public :
         void  solve( bool  ptl) {
             bool  moveOn =  true ;
             int  EFT = 0;
             rebootCount = 0;
             reCount();  // initialization process
             startCTime =  clock ();
             while  (!finish()) {  // Terminate
                 reboot();
                 rebootCount++;
                 reCount();
                 EFT=0; moveOn= true ;
                 while  (moveOn) {
                     moveOn =  false ;
                     for  ( int  i=1; i<=N; ++i) {
                         for  ( int  j=1; j<i; ++j) {
                             // -------------beginforeach i, j
                             EFT = 0;
                             // for any 'i' col & 'j' col
                             if  (collapsed(i, j)) {  // (i, j) pair collapsed
                                 EFT = getChange(i, j);
                             else  EFT = -1;
                             if  (EFT>0) {  // still trends to minimal
                                 moveOn =  true ;
                                 // do position swapping
                                 // maintain the array 'left' & 'right'
                                 execSwap(i, j);
                             }
                             reCount();
                         // -------------------------- endforeach i, j
                         }
                     }
                 }
                 if  (ptl) {
                     puff();
                     prinf();
                 }
             }
             endCTime =  clock ();
             costCTime = endCTime - startCTime;
         }

    图 6 主要搜索方法

    用两个主要变量控制搜索进程:

    int EFT用来衡量每次可交换列的对(i, j)交换后的估值函数值;

    bool moveOn判断当前循环是否有效,即是否有正向估值产生。

    每当一步搜索不能产生正向估值,即无论如何移动也不能使总的冲突数减小时,一个新的局部极值就可能被找到了。这时我们若判断出总冲突数降至0则可以进一步说明这个极值也是一个全局最优。


    Local search strategy

    贪心与模拟退火:

    爬山问题中,我们判断当前搜索方向是否还会继续趋向于更大值。若当前搜索方向上不再出现上升态势,搜索会停止并驻留在当前解上(即局部最优解),所以爬山策略是一种严格贪心算法,不保证可以找到全局最优。

    模拟退火算法正是通过增加一定几率的反向跳转趋势,使得搜索位置得以“逃脱”局部最优,获得移动到全局最优的可能;通过控制反向跳转几率的收敛,使得解有更大可能性最终停留在全局最优解上。

    局部搜索策略:

    皇后问题中,局部最优解的分布较为稠密,全局最优解的比例也相对较高。这使得模拟退火等一般意义上的概率收敛不一定能取得较好的效果。


    Simple Statistic

    简单统计结果:

    下面给一个运行的例子 

    50000次运行求出最优解,平均每个最优解的求得只需要2.6次重启,即经过2.6个局部最优解到达了全局最优;启动最多的也只有25次;平均耗时极短。  


    程序打印出的一些局部极值示意图微笑




    ps. 最近新开通了博客 会放一些学习过程中的心得体会(和一些作业或者报告)

    如果你觉得有用,欢迎分享,谢谢!

    展开全文
  • 局部最小值和全局最小值

    千次阅读 2018-06-29 15:19:10
    基于梯度得搜索是使用最为广泛得参数寻优方法。在此类方法中,我们从...若误差函数在当前点的梯度为零,则已达到局部极小,更新量将为零,这意味着参数的迭代更新将在此停止。显然,如果误差函数仅有一个局部极小,...
  • close figure(2) syms x y f=2.*x.^4+y.^4-2.*x.^2-2.*y.^2+3; fx=diff(f,x); fy=diff(f,y); fxx=diff(f,x,2); fyy=diff(fy,y); fxy=diff(fx,y); fx fy fxx fyy fxy [a,b]=solve(fx,fy) ...tmp=2...
  • 局部极值的问题的全局最优问题还没有根本最优解决,做的也是多找一些可能是全局最优解的局部最优结,再从中选择最打的那个,认为是最优的...... 单这样的问题就很复杂了,实际中有会有更多的问题..... 我有几个
  • 算法将布谷鸟全局搜索能力与Powell方法的局部寻优性能有机地结合,并根据适应度值逐步构建精英种群候选解池在迭代后期牵引Powell搜索的局部优化,在保证求解速度、尽可能找到全局极值点的同时提高算法的求解精度。...
  • 极值局部性概念 最值 :全局性概念 比较 : 4 ,极值定理 : 极值与导数 : 1 ,如果 y = f(x) 在 x0 出有极值,并且可导,那么 导数 = 0 2 ,反之不成立 3 ,结论 : 极值,一定导数等于 0 ,导数等于 0 ,...
  • 一、问题描述 无约束的多维极值问题一般描述如下公式...对于实际问题来说,这是不矛盾的,因为实际问题都存在一定的应用背景使用条件,局部最优点并不多,甚至有时局部最优点就是全局最优点,所以实际问题可以根据实
  • 局部特征融合为全局特征笔记

    千次阅读 2018-03-29 20:46:00
    图像特征分为全局特征和局部特征两种,其中全局特征代表了图像的整体表现特性,比如颜色直方图,而局部特征代表了图像的局部特性,往往能够从一幅图片中提取出若干个数量不等的局部特征, 这些局部特征组合起来代表...
  • 全局描述符表GDT和局部描述符表LDT

    千次阅读 2018-09-24 14:14:33
    段描述符  IA-32架构的处理器访问内存都是采用“段基址:段内偏移地址”的形式,即使到了保护模式也不例外。其次,实模式脆弱的安全性也是... 一个段描述符是8字节大小,它描述了一个内存段的地址范围各种属性。...
  • 算法将布谷鸟全局搜索能力与Powell方法的局部寻优性能有机地结合,并根据适应度值逐步构建精英种群候选解池在迭代后期牵引Powell搜索的局部优化,在保证求解速度、尽可能找到全局极值点的同时提高算法的求解精度。...
  • 最大稳定极值区域MSERs

    万次阅读 2013-10-28 17:36:23
    Harris兴趣点检测器是一种产生平移旋转不变结果的算法。...最大稳定极值区域(Maximally Stable Extremal Regions MSERs)是一种图像结构,不仅是是在平移旋转后,即便是经历相似仿射变换,它仍可被重复检测出来。
  • 空中签名序列长,为了解决传统的全局匹配方法造成的匹配慢、签名的局部信息丢失的问题,提出了对签名数据进行极值点分段再进行距离度量的方法。并针对传统DTW算法在极值点匹配中产生的不同极性极值点错匹配问题,...
  • 序 说起非线性最优化问题,我们并不陌生。初等数学里,二次函数的极值问题,就是典型的非线性最优化问题。比如,求函数 f(x)=x2−2x+4f...全局极小点、局部极小点 利用极值点一阶导数为零的特点,可以求出目标函数的
  • 无约束问题的极值条件

    千次阅读 2014-03-25 18:37:15
    有时候,我们希望根据一定的条件找到优化问题的极值点;另外一些时候,我们得到若干候选解,希望判断候选...1. 极值点的必要条件充分条件  一阶必要条件 设实值函数 在点 处可微,若是无约束优化问题 的局部极小点
  • 在数学分析中,函数的最大值最小值(最大值最小值)被统称为极值(极数),是给定范围内的函数的最大值最小值(本地 或相对极值)或函数的整个定义域(全局或绝对极值)。
  • 给出了判断局部极值和全局极值点的方法,解决了第二个问题。应用相对微分/差分法解连续离散非线性规划,在搜索过程中一旦满足了两个充分条件之一,就达到了极值点。根据搜索方向很容易确定极值点是极大点还是极...
  • 极值:动态规划、贪心、回溯。 这里涉及的内容很多,不是一篇博客能说清楚的,这里宏观上,对它们做个对比总结。 1.1 插入排序 直接插入排序:对待排序的记录逐个进行处理,每个新记录与同组那些已经排好序的...
  • 即不能检测出小于一定内在尺度的社区,并提出了基于极值优化模块密度来检测复杂网络社区结构的启发式算法,通过调整局部极值来优化全局的变量,使算法具有更好的持续搜索跳出局优解的能力.通过人工网络现实网络实验...
  • 全局图像特征与局部图像特征区别

    万次阅读 2018-03-10 09:41:41
    全局特征是指图像的整体属性,常见的全局特征包括颜色特征、纹理特征形状特征,比如强度直方图等。由于是像素级的低层可视特征,因此,全局特征具有良好的不变性、计算简单、表示直观等特点,但特征维数高、计算量...
  • 凸函数的极值和凸规划 文章目录 凸函数的极值和凸规划 凸函数的极值 凸规划 定理:由f、g、h凹凸性得是否为凸规划 证明:凸规划的最优解集必是凸集 线性规划的整体最小点与局部极小点 凸规划局部极小点性质 性质1:...
  • matlab求最值(极值

    万次阅读 多人点赞 2016-09-25 18:26:55
    根据我看过的材料来说 ,求最值(极值)无非有下面几种情况。 1、求简单函数的单一最值(极值) clear clc t=-100:0.001:100;%初值:增量:终值 symsx; y=x/(x*x+1); f=inline(y);%内联函数 max=max(f(t)) min=...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 6,539
精华内容 2,615
关键字:

局部极值和全局极值