Jan 21, 2019 - 伪随机数与任意分布的生成

背景

计算机软件里有大量需要生成随机数的场景,更准确一点,需要随机性。比如蒙特卡洛方法这样的模拟方法,Bootstrap中经验分布随机采样等等…如何依据想要的分布生成随机数,以及验证随机性是计算机时代大量工作的基础。

线性同余发生器(Linear Congruential Generator)

由递归式定义: X_n+1 = (a*X_n+c) mod m a称为乘数,c为步数, m为模数。

当给定x_0,也就是代码中经常需要设置的seed,随后的”随机数”就是确定的,因此称为伪随机。

在任意给定的<a,c,m>三元组中,由于X_i必然小于m,而X_i可以唯一确定X_i+1, 因此任意上述公式生成的递推序列,必然是有限循环的,并且周期小于等于m。

下面给出线性同余发生器,可取得最大理论周期m的三个条件:

  • gcd(c, m) = 1
  • 因式分解m所得的质数,是因式分解(a-1)所得的质数的子集。(Each prime factor of m is also a factor of(a-1) )
  • 如果m整除4, (a-1)整除4

具体证明略(因为我也不懂为什么)

分布

均匀分布

由线性同余发生器产生的数字,可以得到一个随机的值域为[0, m-1]的序列。当把产生的随机数/(m-1),就映射到了[0,1]的区间上

任意分布

以上,可以产生均匀分布,均匀分布的值可以理解成任意0-1的概率,于是可以求反函数,将目标分布的累积分布函数解出应有的函数参数值,将随机性映射到目标分布随机变量的取值。

Box-Muller到正态分布

转换到正态分布并不能使用累积分布函数取反函数的方法,因为正态分布的CDF没有闭式解。BM使用极坐标下的变换,使用两个独立的均匀分布随机变量得到一个二维联合正态分布,且两者独立,因此就可以得到正态分布。

Jan 2, 2019 - 【指南向】ubuntu18.04装机tensorflow-gpu

背景

找计算资源跑CNN demo. rmbp没有gpu,跑一个mini batch都要好几秒(单次trainning CPU狂转估计得跑十来分钟)。阿里云/腾讯云的GPU实例贵得吓人,有一个极客云很便宜,但是感觉很不正规。调查下来,Google Cloud注册送钱,生成一个最低配的实例,这些钱可以用很久。于是就他了。

我选的安装栈

  • ubuntu18
  • nvidia driver
  • cuda toolkits(主要是cuda9.0, cudnn)
  • miniconda
  • tensorflow-gpu

废话

选ubuntu是因为据说跑DL用得最多,出问题好找,选择了最新的ubuntu。 nvidia driver与cuda套件据说比较经常在nvidia官网找版本手动安装,但是我嫌麻烦,用的是ubuntu-drivers工具。装miniconda是为了方便装tensorflow-gpu系列和cuda系列。没有选择anaconda是因为整个套件太大了(12G)…磁盘占用是要钱的。总之,就是在拒绝任何手工编译的原则下,完成安装。tensorflow官网有指导ubuntu16安装gpu套件,因为发行版不匹配我没有尝试…

来了

首先在google cloud平台上生成一个VM实例,我选的是单核CPU+2G内存,(因为计划跑python…多核应该用处不大),GPU不是每个机房都有,选台湾,1个GPU+12G内存,Nvidia Tesla K80,基本就是最便宜的套路了。系统选ubuntu+20G SSD(磁盘选太小了其实)。(后来发现VM Instance关机的时候是不收钱的…其实可以创建一个配置更高的实例,不用的时候关机就好了)

这里有个麻烦的事情是,google cloud在IAM管理-配额里面可以看到GPU数量限制是0,需要在他的界面上请求加数量,(我加的是1),据说加多了会要求充值,不让白嫖…这个审核是人工的,可能得要一点时间,我差不多等了一个小时不到。

管理这个VM,从网页命令行也可以,我使用了gcloud在我本机的命令行远程操作,需要看一点谷歌的说明配置命令。

安装显卡驱动

参考https://linuxconfig.org/how-to-install-the-nvidia-drivers-on-ubuntu-18-04-bionic-beaver-linux安装N卡驱动,我这里的版本是nvidia-390,小版本390.77,然后重启。

主要的命令是

$ ubuntu-drivers devices
$ sudo ubuntu-drivers autoinstall

完事重启以后,

$ nvidia-smi

查看是否成功

运行nvidia-docker检查一下驱动

这里我选择先跑跑看docker验证一下驱动是不是装好了…免得再然后报错不知道到底是驱动的问题还是cuda套件的问题

顺着docker官方指南,安装docker-ce

https://docs.docker.com/install/linux/docker-ce/ubuntu/#set-up-the-repository

安装完docker-ce,按nvidia-docker readme安装支持

https://github.com/NVIDIA/nvidia-docker

跑完readme的安装与验证不算,我还跑了一个tensorflow官网上的docker实例

https://www.tensorflow.org/install/docker

$ docker run -it --rm tensorflow/tensorflow \
   python -c "import tensorflow as tf; tf.enable_eager_execution(); print(tf.reduce_sum(tf.random_normal([1000, 1000])))"

完美运行…心满意足去装cuda

安装tensorflow-gpu运行环境

主要需要安装的是cuda和cuDNN,cuDNN就是专门为深度神经网络加速的一个组件。安装了conda之后,这两个东西都可以傻瓜安装,无需编译。

安装conda

https://conda.io/miniconda.html

执行官方的安装脚本,顺着他安装。

完事之后,创建一个conda env

$ conda create -n test python=3.6 
$ conda activate test

然后参考https://mc.ai/install-tensorflow-gpu-to-use-nvidia-gpu-on-ubuntu-18-04-do-ai/ 最后部分,使用conda安装所有的python开发套件

$ conda install \
 tensorflow-gpu==1.10.0 \
 cudatoolkit==9.0 \
 cudnn=7.1.2 \
 h5py

按照文章中的简单python脚本,测试tensorflow在gpu上的运行

import tensorflow as tf
print("tf version = ", tf.__version__)
with tf.device('/gpu:0'):
    a = tf.constant([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], shape=[2, 3], name='a')
    b = tf.constant([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], shape=[3, 2], name='b')
    c = tf.matmul(a, b)

with tf.Session() as sess:
    print (sess.run(c))

再来一个tensorflow官方给的gpu demo(这个其实我没试)

https://www.tensorflow.org/guide/using_gpu

# Creates a graph.
with tf.device('/device:GPU:0'):
  a = tf.constant([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], shape=[2, 3], name='a')
  b = tf.constant([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], shape=[3, 2], name='b')
  c = tf.matmul(a, b)
# Creates a session with log_device_placement set to True.
sess = tf.Session(config=tf.ConfigProto(log_device_placement=True))
# Runs the op.
print(sess.run(c))

跑一个正式的demo

运行https://github.com/astorfi/TensorFlow-World/tree/master/codes/3-neural_networks/convolutional-neural-network 中CNN在MINIST数据集上的例子(也就是我在本机CPU版本上十几分钟都跑不完的demo)

注意,在train.sh中修改--log_device_placement置为True,观察gpu是否被使用了。

结果十几秒就跑完了一次trainning…相比之下,我RBMP上的3.1GHz的CPU也算是苹果笔记本里能打的了…跑浮点GPU厉害太多了…

总结

下次磁盘可以开大点…CPU,内存,GPU都是可以事后加的,磁盘分区不能直接扩容,挺麻烦的。

以及

白嫖,爽~

Sep 21, 2018 - cpython字典源码(一)

本文源代码基于 commit: 73820a60cc3c990abb351540ca27bf7689bce8ac

哈希表:开放定址法

As we all know,在需要快速搜索的场景下,哈希表提供O(1)的检索复杂度。哈希就可能会产生冲突。对冲突的处理基本上分为两大类:链表法(又称拉链法),开放定址检测法。 链表法,就是说在每个哈希槽(slot)上其实都是一个链表表头,所有落入这个槽的元素都接入这个链表。

开放定址检测法,本质就是多试几次,万一可以呢,发现本key的哈希槽被占用的时候,使用一个算法算一算下一个槽的位置,可以就占掉,不可以就再试下去。

个人觉得,这两种方法在实现和使用的时候,各有不同。链表法实现较为直观,但可能造成空间浪费。想象申请了n块内存,结果哈希值比较稀疏,值都落到了几个大链表里,就不划算。另一方面,链表法什么时机重hash也是一个问题(个人感觉),毕竟单个槽过长,搜索会退化成O(n)。开放定址检测不存在内存浪费,申请多少用多少,在扩容之前只会用这些申请的内容。重hash的时机也比较好掌握;但是在设计检测算法的时候也要注意,首要是不能产生死循环(在有足够空位的情况下,要能最终找到至少一个空位),其次二次检索的效率也要注意。比如线性检测就不太好,在哈希槽较满的情况下容易产生近O(n)的效率。

cpython中解决字典的哈希冲突使用的是开放定址法。

开放定址法的删除

链表法中要删除一个元素很容易,把他从哈希槽的链表里拿掉就好了。但开放定址法不能简单地删掉。举个例子,用线性探测(虽然不好)时,搜索其实构成了一个搜索顺序。比如搜索key为13, hash值为3(比如hash函数就是n%10)的元素能插入哪个位置,3,4,5都已经有元素了,那么搜索顺序会有

    3(have) --> 4(have) --> 5(have) --> 6(null)

ok, 找到了我的家(yeah),那么key 13就落户编号6的哈希槽。但是这个时候,用户删除了位置4的哈希槽元素,然后再去搜索key 13元素位置的时候会变成这样:

    3(have) --> 4(null) --> 5(have) --> 6(have 13)

到位置4的时候发现没有了。可能有人会问,那算法继续搜索下去就好了嘛,总会找到key 13的,但这是在知道这个key存在的情况下才有意义;如果算法这么实现的话,所有搜索不存在key的情形都会把搜索算法变成O(n),因此开放定址法的搜索止于找到该key或发现一个空位。如上,这个搜索空间就断裂了。(这个问题在cpython源码注释中也有说明)

所以开放定址法不能直接删除元素,而是应该标记删除,以确保搜索空间的完整连续性,构造出如下的搜索空间:

    3(have) --> 4(dummy) --> 5(have) --> 6(have 13)

其中,dummy就是cpython中使用的标记删除,算法跑到4处会继续向下搜索直到6的位置

cpython的搜索例程

上代码

static Py_ssize_t _Py_HOT_FUNCTION
lookdict(PyDictObject *mp, PyObject *key,
         Py_hash_t hash, PyObject **value_addr)

cpython中关于在dict中搜索的代码位于dictobject.c,函数签名如上。参数中传入一个操作的字典对象mp,查找的key,该key的hash值(传入是因为之前算过,就不要重复算了),最后一个参数是查找出来的值的引用,需要在函数内赋值(或者没找到置空)。返回值就是该key的位置,有三种可能:

  • DKIX_EMPTY
  • DKIX_DUMMY
  • 某个大于等于0的数字。

DKIX_EMPTY就是说没找到,DKIX_DUMMY就是说找到了这个key的位置,但是已经被删除了,最后就是找到的哈希槽编号的数值了。

size_t i, mask, perturb;
    PyDictKeysObject *dk;
    PyDictKeyEntry *ep0;

top:
    dk = mp->ma_keys;
    ep0 = DK_ENTRIES(dk);
    mask = DK_MASK(dk);
    perturb = hash;
    i = (size_t)hash & mask;

初始了一坨变量,其中hash & mask这个操作比较特殊。其中hash就是哈希值,mask等于dk_size-1。这个操作让i落在哈希槽范围内。这个地方当时粗看的时候我觉得很诡异,为啥与操作能到达跟取余一样的效果呢。这里不是说与操作的结果值跟取余是一样的,而是说能达到同样的收敛值域的作用。然后我做了一个不加验证的假设: cpython中字典的容量都是2^n,然后取2^n-1,让n位二进制位都是1,这样与操作,x & (2^n-1)的结果必然小于等于2^n-1。此外这段代码设置了一个tag:top,目的是在有必要的时候goto回来重新开始。

注意,这里cpython的实现是for(;;)死循环,实现上要求必要有返回。下面分几段代码分析一个例程。

    for (;;) {
        Py_ssize_t ix = dk_get_index(dk, i);
        if (ix == DKIX_EMPTY) {
            *value_addr = NULL;
            return ix;
        }
        if (ix >= 0) {
            ...
        }
        perturb >>= PERTURB_SHIFT;
        i = (i*5 + perturb + 1) & mask;
    }

首先在每个循环例程最开始,都会用当前的i来找内存该位置的状态,如果发现这个位置是空,就直接返回了。如果发现有东西,就会进一步做一些操作来确认这块内容是不是我想要的。如果又不是DKIX_EMPTY,又不是大于等于0的某个数,按照上述的约定,就一定是DKIX_DUMMY,也就是找到了一个已删除的内容,继续做探测。

   for (;;) {
        Py_ssize_t ix = dk_get_index(dk, i);
        if (ix == DKIX_EMPTY) {
            ...
        }
        if (ix >= 0) {
            PyDictKeyEntry *ep = &ep0[ix];
            assert(ep->me_key != NULL);
            if (ep->me_key == key) {
                *value_addr = ep->me_value;
                return ix;
            }
            if (ep->me_hash == hash) {
                PyObject *startkey = ep->me_key;
                Py_INCREF(startkey);
                int cmp = PyObject_RichCompareBool(startkey, key, Py_EQ);
                Py_DECREF(startkey);
                if (cmp < 0) {
                    *value_addr = NULL;
                    return DKIX_ERROR;
                }
                if (dk == mp->ma_keys && ep->me_key == startkey) {
                    if (cmp > 0) {
                        *value_addr = ep->me_value;
                        return ix;
                    }
                }
                else {
                    /* The dict was mutated, restart */
                    goto top;
                }
            }
        }
       ...
    }

以上代码就是找到了一个非空的哈希槽之后该做什么。 首先比较了搜索的key与该哈希槽存储的key是不是根本同一个指针。因为cpython会复用很多小对象。 其次会比较这个槽存储元素的hash值。因为经过哈希碰撞之后,每个槽存储的hash值并不一定是一手取与得到的。如果不一样,for循环就会进入下一个例程。has值一样,然后要比较该元素的key是否deepEqual当前的key。因为都说了哈希碰撞嘛,这个key跟key也可能是不一样但是hash值一样的。(连SHA1都可以碰撞的说…)。经过hash值一样,key相等比较之后,这个元素位”基本”上就是我们要找的元素了。

但是,这里cpython做了”看似”比较多余的事情,他比较了在代码中使用的字典还是不是要搜索的字典,被搜索到的key还是不是之前比较的key。这里用局部变量跟外部变量引用相比较,是要确定外部变量有没有被修改。这个代码为什么可以生效,目前我还没有明确的证据,因为就算这个if (dk == mp->ma_keys && ep->me_key == startkey)能成立,进入if之后就能保证这个if依然成立吗;返回值依然有可能是脏的。还需要进一步阅读代码。

检测算法

上一节的代码里,其实我真正关心的只有下面这段做检测的代码。从函数关系来说,这个探测是一种平方探测。(突然想起了时间序列里的AR(1)…)。perturb初值就是hash值,每次循环右移PERTURB_SHIFT = 5位参与运算

        perturb >>= PERTURB_SHIFT;
        i = (i*5 + perturb + 1) & mask;

怎么证明这个函数关系不会陷入死循环呢?恩,我也不知道…严格证明懒,留一个随缘未解之谜吧…

我写了一个go的小代码来看这个算法会”路过”那些哈希槽。

package main

import (
	"fmt"
)

func main() {
	hash := 9
	perturb := hash
	mask := 7
	i := hash & mask
	for j:=0;j<20;j++{
		fmt.Println(i)
		perturb >>= 5
		i = cal(i, perturb, mask)
	}
}

func cal(i int, perturb int, mask int) int {
	return (i*5 + perturb + 1) & mask
}

输出的搜索空间如下:

1 -->6 --> 7 --> 4 --> 5 --> 2 --> 3 --> 0

可以看到,在这个hash:9 ,mask:7输入下,可以一次性探测到所有的8个位置。后面的输出就是重复这个序列了。(是不是很神奇呢…)

Ending

好了,本文到此为止,字典的CRUD就只介绍R的核心部分。(因为写文章很花时间的…)其他内容,有缘有心情再见。