0和1

最近开始学习神经网络和深度学习，使用的是网上教程：http://neuralnetworksanddeeplearning.com/，这是学习心得第一讲，介绍经典的MNIST手写数字图片数据集。 MNIST（Modified National Institute of Standards and Technology database）数据集改编自美国国家标准与技术研究所收集的更大的NIST数据集，该数据集来自250个不同人手写的数字图片，一半是人口普查局的工作人员，一半是高中生。该数据集包括60000张训练集图片和10000张测试集图片，训练集和测试集都提供了正确答案。每张图片都是28×28=784大小的灰度图片，也就是一个28×28的矩阵，里面每个值是一个像素点，值在[0,1]之间，0表示白色，1表示黑色，(0,1)之间表示不同的灰度。下面是该数据集中的一些手写数字图片，可以有一个感性的认识。 MNIST数据集可以在Yann LeCun的网站上下载到：http://yann.lecun.com/exdb/mnist/，但是他提供的MNIST数据集格式比较复杂，需要自己写代码进行解析。目前很多深度学习框架都自带了MNIST数据集，比较流行的是转换为pkl格式的版本：http://deeplearning.net/data/mnist/mnist.pkl.gz，该版本把原始的60000张训练集进一步划分成了50000张小训练集和10000张验证集，下面以这个版本为例进行介绍。 pkl是python内置的一种格式，可以将python的各种数据结构序列化存储到磁盘中，需要时又可以读取并反序列化到内存中。mnist.pkl.gz做了两次操作，先pkl序列化，再gz压缩存储，所以要读取该文件，需要先解压再反序列化，在python3中，读取mnist.pkl.gz的方式如下： 1 2 3 4 5 import pickle import gzip f = gzip.open(‘../data/mnist.pkl.gz’, ‘rb’) training_data, validation_data, test_data = pickle.load(f, encoding=’bytes’) f.close() 这样就得到了训练集、验证集和测试集。将数据集序列化到文件中的方法也很简单，需要注意的是pickle在序列化和反序列化时有不同的协议，可以用protocol参数进行设置。 1 2 3 4 dataset=[training_data, validation_data, test_data] f=gzip.open(‘../data/mnist3.pkl.gz’,’wb’) pickle.dump(dataset,f,protocol=3) f.close() 我们从mnist.pkl.gz读取到的training_data, validation_data, test_data这三个数据的结构是一样的，每个都是一个二维的tuple。以training_data为例，training_data[0]是训练样本，是一个50000×784的矩阵，表示有50000个训练样本，每个训练样本是一个784的一维数组，784就是把一张28×28的图片展开reshape成的一维数组；training_data[1]是训练样本对应的类标号，大小为50000的一维数组，每个值为0~9中的某个数，表示对应样本的数字标号。 ...

VSCode是微软开源的一个很强大的IDE，可以支持几乎所有编程语言，而且是跨平台的，Linux用户终于可以用上宇宙最强IDE了。我最近在使用VSCode编写调试Python项目，其调试功能很强大，和VS上调试C++的感觉是一样的，强烈推荐。 VSCode还可以连接Github，进行版本控制。下面以我最近学习的深度学习项目为例，介绍下怎样在Ubuntu下使用VSCode连接Github。以我fork的repo为例：https://github.com/01joy/neural-networks-and-deep-learning。 连接Github有两种方式，一种是HTTPS，另一种是SSH，在每个repo页面的右边，有一个Clone or download按钮，可以获取到这两种连接方式的地址。HTTPS方式和网址类似，以HTTPS开头；SSH方式以git@githu.com开头。使用HTTPS连接比较简单，但是每次push的时候需要输入用户名和密码，比较麻烦，如果想记住密码，需要把用户名和密码以明文的形式保存到一个文件中，个人感觉不方便且不安全。下面以SSH连接为例进行介绍。 首先设置Github提交时的用户名和密码，一般设置成全局的：https://help.github.com/articles/setting-your-username-in-git/、https://help.github.com/articles/setting-your-commit-email-address-in-git/ 生成一对新的SSH公钥和私钥，并添加到ssh-agent中。注意生成的时候需要输入passphrase，这个passphrase不是Github的密码，自己随便取一个记住就好。https://help.github.com/articles/generating-a-new-ssh-key-and-adding-it-to-the-ssh-agent/ 把SSH的公钥添加到Github账号中：https://help.github.com/articles/adding-a-new-ssh-key-to-your-github-account/ 测试SSH连接是否成功：https://help.github.com/articles/testing-your-ssh-connection/ （可选）修改SSH密码，即第1步设置的passphrase：https://help.github.com/articles/working-with-ssh-key-passphrases/ 到这里，本级就能通过SSH连接Github了。 如果没有安装VSCode，可以直接通过Ubuntu的终端连接Github，步骤如下： 在本地创建一个和远程repo名称一样的空文件夹 终端cd到该文件夹内 git init # 在该文件夹内初始化 git remote add origin git@github.com:01joy/notes-on-writing.git # 使用repo的SSH地址 git pull origin master # 把远程代码拉到本地 修改代码 git add . # 在根目录执行，添加所有修改 git commit -m ‘comments’ # commit第7步添加的修改 git push origin master # 把第8步发布到远程 如果安装了VSCode，其实和直接用终端是一样的，在菜单栏的Terminal下新建一个终端，在这个终端内执行上述代码，如果在第4步出现”Enter password to unlock the private key”时，输入创建SSH时第2步的密码即可，只需一次，下次就不用再输入密码了。点击File的Open Folder打开本地repo文件夹。点击VSCode左边栏的Explorer可以在编辑器下修改代码。切换到左边栏的Source Control可以进行Git相关操作，修改的文件右边会出现一个M，点击这个M会出现diff视图；Source Control左边的右上角有三个点，点击这个按钮会出现很多Git操作，包括commit、push等，其实相当于调用上述代码，效果是一样的。 VSCode快捷键Ctrl+Shift+P会出现命令窗口，在里面输入commit、push等会出现相关操作的，能起到一定的加速效果，当然也可以自定义快捷键。 Have Fun!

上一篇博客主要介绍了逻辑回归的理论知识，这篇博客咱们用Python机器学习包sklearn中的LogisticRegression做一个分类的实例。 数据还是学生样本，只有两个特征，分别是两门课的分数score1和score2，类标号y表示这个学生是否能被录取。先上分类效果图： 完整的Python代码如下： 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 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 # -*- coding: utf-8 -*- """ Created on Wed Nov 08 17:49:41 2017 @author: zhenlin """ import numpy as np import pandas as pd from sklearn.cross_validation import train_test_split from sklearn.linear_model import LogisticRegression import matplotlib.pyplot as plt from sklearn.metrics import precision_recall_curve from sklearn.metrics import classification_report # 1. 构造数据 sample_number = 200 # 第一个高斯分布参数 mean1 = [0, 4] # 两个维度上的均值 cov1 = [[5, 3], [3, 10]] # 两个维度的协方差矩阵，必须满足对称半正定 # 第二个高斯分布参数 mean2 = [7, 5] cov2 = [[7, 2], [2, 15]] # 从两个二元高斯分布中随机采样数据点 class1_x1, class1_x2 = np.random.multivariate_normal(mean1, cov1, sample_number).T # .T表示转置 class2_x1, class2_x2 = np.random.multivariate_normal(mean2, cov2, sample_number).T # 两个高斯分布对应两个类标号 data = [[class1_x1[i],class1_x2[i],0] for i in range(sample_number)]+[[class2_x1[i],class2_x2[i],1] for i in range(sample_number)] # 填充到pandas中 data = pd.DataFrame(data,columns=['score1','score2','result']) score_data = data[['score1','score2']] result_data = data['result'] # 2. 训练模型 average_precision = 0 # 平均准确度 iters = 10 # 交叉验证次数 for i in xrange(iters): # 数据划分，80%用于训练，20%用于预测 x_train, x_test, y_train, y_test = train_test_split(score_data, result_data, test_size = 0.2) # 构造默认逻辑回归模型 model = LogisticRegression() # 训练 model.fit(x_train, y_train) # 预测 predict_y = model.predict(x_test) # 计算测试集上的准确度 average_precision += np.mean(predict_y == y_test) average_precision /= iters # 3. 绘制分类面 - 法1 x1_min, x1_max = score_data['score1'].min() - .5, score_data['score1'].max() + .5 def generate_face(prob): y = -np.log(1.0 / prob - 1.0) n = 500 x1 = np.linspace(x1_min, x1_max, n) # w1x1+w2x2+b=y x2 = (-model.coef_[0][0] / float(model.coef_[0][1])) * x1 + (y - model.intercept_) / float(model.coef_[0][1]) return x1, x2 pos_data = data[data['result'] == 1] neg_data = data[data['result'] == 0] plt.scatter(x = pos_data['score1'], y = pos_data['score2'], color = 'black', marker = 'o') plt.scatter(x = neg_data['score1'], y = neg_data['score2'], color = 'red', marker = '*') face_04_x1, face_04_x2 = generate_face(0.4) face_05_x1, face_05_x2 = generate_face(0.5) face_06_x1, face_06_x2 = generate_face(0.6) plt.plot(face_04_x1, face_04_x2) plt.plot(face_05_x1, face_05_x2) plt.plot(face_06_x1, face_06_x2) plt.xlim(score_data['score1'].min(), score_data['score1'].max()) plt.ylim(score_data['score2'].min(), score_data['score2'].max()) plt.xlabel('score1') plt.ylabel('score2') plt.legend(['prob_threshold = 0.4', 'prob_threshold = 0.5', 'prob_threshold = 0.6'], loc='center left', bbox_to_anchor=(1, 0.865)) plt.show() # 4. 绘制分类面 - 法2 pos_data = data[data['result'] == 1] neg_data = data[data['result'] == 0] h = 0.02 s1_min, s1_max = score_data['score1'].min() - .5, score_data['score1'].max() + .5 s2_min, s2_max = score_data['score2'].min() - .5, score_data['score2'].max() + .5 # 生成s1在[s1_min, s1_max]，且s2在[s2_min, s2_max]的网格数据点 # meshgrid含义参见：http://blog.sciencenet.cn/blog-791749-675394.html s1, s2 = np.meshgrid(np.arange(s1_min, s1_max, h), np.arange(s2_min, s2_max, h)) # 把两个坐标的值按列拼在一起构成二维数据点 Z = model.predict(np.c_[s1.ravel(), s2.ravel()]) # 绘制边界和散点 Z = Z.reshape(s1.shape) # 坐标点是(s1[i], s2[i])，对应颜色是Z[i]，颜色主题使用plt.cm.Paired plt.pcolormesh(s1, s2, Z, cmap = plt.cm.Paired) plt.scatter(x = pos_data['score1'], y = pos_data['score2'], color = 'black', marker = 'o') plt.scatter(x = neg_data['score1'], y = neg_data['score2'], color = 'red', marker = '*') plt.xlim(s1.min(), s1.max()) plt.ylim(s2.min(), s2.max()) plt.xlabel('score1') plt.ylabel('score2') plt.show() # 5. 评估模型 # 对于测试数据，模型输出1的概率 answer = model.predict_proba(x_test)[:,1] # 计算不同概率阈值下的P和R precision, recall, thresholds = precision_recall_curve(y_test, answer) # prob > 0.5的报告为1 report = answer > 0.5 print(classification_report(y_test, report, target_names = ['neg', 'pos'])) print('average precision: %f'%average_precision) # 6. 绘制PRC曲线 # step阶跃图，在点(recall[i],precision[i])进行跳变 plt.step(recall, precision, color='b', alpha=0.2, where='post') # 对PRC下方填充颜色 plt.fill_between(recall, precision, step='post', alpha=0.2, color='b') plt.xlabel('Recall') plt.ylabel('Precision') plt.ylim([0.0, 1.05]) plt.xlim([0.0, 1.0]) plt.title('2-class Precision-Recall curve') plt.show() 下面将逐模块介绍代码细节，大神可以略过。 ...

最近实验室在组织学习NG的机器学习视频，我也跟进了一下。讲到逻辑回归那一课，一直想不明白，逻辑回归到底是怎样分类的？逻辑回归的分类面在哪里？逻辑回归有像SVM的max margin那样清晰的推导过程吗？为什么需要Sigmoid函数？今天就让我们来一窥逻辑回归的始末。 假设有一堆学生样本，为了简单起见，只有两个特征，分别是两门课的分数score1和score2，类标号y表示这个学生是否能被录取。可视化如下图，黑点表示y=1即被录取，红点表示y=0即未被录取，可以看到黑点处在score1和score2都比较高的区域。我们的任务就是给定这些训练样本，需要确定一个分类面来划分这两类数据。 分类面的实质就是\(y=\mathbf{w^T x}+b\)，其中\(\mathbf{w}\)和\(\mathbf{x}\)都是向量，对应到本例中，展开为\(y=w_1x_1+w_2x_2+b\)。所以，寻找分类面的过程就是寻找倾斜度\(\mathbf{w}\)和截距\(b\)。 因为分类的结果是离散的，只有0或者1，可以用感知机来分类。即 但是怎样找这里的\(\mathbf{w}\)和\(b\)使得分类结果最好呢，我们需要定义一个优化的目标函数或者说损失函数。这里只能定义为分类错误的个数，只能一点点挪动超平面来找分类错误最少的超平面了，即只能用暴力枚举的方法来找\(\mathbf{w}\)和\(b\)。 另一种方法是令我们的分类面算出来的值和真实标号越接近越好，即最小化误差\((\mathbf{w^T}x^{(i)}+b-y^{(i)})^2\)，然后通过梯度下降求解\(\mathbf{w}\)和\(b\)。 但是这会有一个问题，对于上图数据，可以求到一个比较好的分类面，如下左图。但是如果数据中出现了如下右图右边的一个离群点或者噪声，为了最小化这个点的误差，即\((\mathbf{w^T}x^{(i)}+b-1)^2\)，导致分类面往右偏了，这一偏直接导致很多y=1的样本分类错误。所以这种最小化误差的方法对离群点很敏感。 假设分类超平面还是\(f(\mathbf{x})=\mathbf{w^T x}+b\)，我们希望的分类效果是这样的：\(f(\mathbf{x})=0\)是分类面；\(f(\mathbf{x})>0\)分类为1，且不管\(f(\mathbf{x})\)多大，都分为1；\(f(\mathbf{x})<0\)分类为0，且不管\(f(\mathbf{x})\)多小，都分为0。 因为类标号是离散的{0,1}，所以想到把\(f(\mathbf{x})\)映射到[0,1]之间，即\(g(f(\mathbf{x}))\)。为了满足上述条件，\(g(f(\mathbf{x}))\)需要满足： \(g(0)=0.5\)，即在分类面上无法判断类标号是0还是1 当\(f(\mathbf{x})>0\)时，\(g(f(\mathbf{x}))>0.5\) 当\(f(\mathbf{x})\rightarrow+\infty\)，\(g(f(\mathbf{x}))\rightarrow 1\)，且\(g'(f(\mathbf{x}))\rightarrow 0\)，即对于上右图右边的离群点，分类为1，且导数趋近于0，表示对其不敏感 当\(f(\mathbf{x})<0\)时，\(g(f(\mathbf{x}))<0.5\) 当\(f(\mathbf{x})\rightarrow-\infty\)，\(g(f(\mathbf{x}))\rightarrow 0\)，且\(g'(f(\mathbf{x}))\rightarrow 0\)，和第3点类似 满足上述性质的函数之一就是Sigmoid函数，其定义域为\([-\infty,+\infty]\)，值域为[0,1]，正好把原始的函数结果\(f(\mathbf{x})\)映射到了[0,1]，而概率的取值范围正好是[0,1]，所以Sigmoid函数正好可以作为分类为1的概率。 所以逻辑回归最终的形式就是： $$g(\mathbf{x})=\frac{1}{1+e^{-(\mathbf{w^T x}+b)}}$$分类面依然还是\(f(\mathbf{x})=\mathbf{w^T x}+b=0\)，因为\(f(\mathbf{x})=0\)时，\(g(\mathbf{x})=0.5\)，正好满足上述条件1。 Sigmoid函数的另一个优点是，它把原始函数值映射到了概率空间，这样就可以用极大似然的方法求解参数\(\mathbf{w}\)和\(b\)。下面用参数\(\mathbf{\theta}\)代表参数\(\mathbf{w}\)和\(b\)，用\(h_{\mathbf{\theta}}(\mathbf{x})\)代表\(g(f(\mathbf{x}))\)。则有： $$P(y=1|\mathbf{x};\mathbf{\theta})=h_{\mathbf{\theta}}(\mathbf{x})$$$$P(y=0|\mathbf{x};\mathbf{\theta})=1-h_{\mathbf{\theta}}(\mathbf{x})$$合并成一个式子就是： $$P(y|\mathbf{x};\mathbf{\theta})=h_{\mathbf{\theta}}(\mathbf{x})^y(1-h_{\mathbf{\theta}}(\mathbf{x}))^{1-y}$$由于所有样本独立同分布（I.I.D.），似然函数就是 $$L(\mathbf{\theta})=P(\mathbf{y}|X;\mathbf{\theta})=\prod\limits_{i}P(y^{(i)}|\mathbf{x}^{(i)};\mathbf{\theta})=\prod\limits_{i}h_{\mathbf{\theta}}(\mathbf{x}^{(i)})^{y^{(i)}}(1-h_{\mathbf{\theta}}(\mathbf{x}^{(i)}))^{1-y^{(i)}}$$最大化似然的含义就是，在给定样本\(X\)的情况下，我们想找一个参数\(\mathbf{\theta}\)，使得观测到类标号\(\mathbf{y}\)的概率最大。 最大化似然等价于最大化log似然，log展开之后就是： $$l(\mathbf{\theta})=logL(\mathbf{\theta})=\sum\limits_{i}y^{(i)}logh_{\mathbf{\theta}}(\mathbf{x}^{(i)})+(1-y^{(i)})log(1-h_{\mathbf{\theta}}(\mathbf{x}^{(i)}))$$而很巧的是，最大化log似然，其实等效于最小化log对数损失。对于单个样本，损失为： $$cost(h_{\theta}(\mathbf{x}),y) = \begin{cases}-log(h_{\theta}(\mathbf{x})) & \text {if y=1} \\ -log(1-h_{\theta}(\mathbf{x})) & \text{if y=0} \end{cases}$$即如果正确类标号是1，但算出来的\(h_{\theta}(\mathbf{x})\)很接近0的话，则损失\(-log(h_{\theta}(\mathbf{x}))\)就会很大。类标号为0的情况类似。把这两种情况合成一个式子就是： $$cost(h_{\theta}(\mathbf{x}),y) = -ylog(h_{\theta}(\mathbf{x})) – (1-y)log(1-h_{\theta}(\mathbf{x}))$$所有样本的损失之和就是： $$J(\mathbf{\theta})=cost(h_{\theta}(X),\mathbf{y}) = \sum\limits_{i}-y^{(i)}logh_{\mathbf{\theta}}(\mathbf{x}^{(i)})-(1-y^{(i)})log(1-h_{\mathbf{\theta}}(\mathbf{x}^{(i)}))$$所以最大化对数似然\(\max l(\mathbf{\theta})\)和最小化对数损失\(\min J(\mathbf{\theta})\)是等效的！最优化求解的方法用梯度下降和牛顿法都可以。 和Sigmoid很像的函数还有很多，比如y=arctan(x)也可以，不过Sigmoid有一个很好的特点，即它的导数可以由自身算出来，\(f'(x)=f(x)(1-f(x))\)。 传统的逻辑回归只能处理二分类问题，那么怎样将其扩展到解决多分类问题呢？有两种方法，一种方法是建立k个二元分类器，比如类标号有A,B,C，则建立3个二元分类器，分别是1）A和非A；2）B和非B；3）C和非C。训练每个2元分类器时，都对所有的数据重新标号为0或1，这样共需要训练k个二元分类器。 还有一种方法是直接将二元逻辑回归推广到多元回归，即Softmax回归，有关Softmax回归的内容，请参考此博客，非常详细。简单来说，二元逻辑回归的假设函数是：多元Softmax回归的假设函数在形式上有所不同：其中是模型的参数。请注意这一项对概率分布进行归一化，使得所有概率之和为1。 Softmax回归的损失函数如下，其实就是logistic回归损失函数的推广： 二元逻辑回归是多元Softmax回归在k=2时的特例，令k=2并利用Softmax回归参数冗余的特点，可以得到一般形式的二元逻辑回归假设函数，具体可以看原博客。 面对一个多元分类问题，是选择Softmax回归还是k个二元分类器呢，这取决于你的类别之间是否互斥，如果互斥，可以用Softmax回归，否则，请用k个二元分类器。 这就是逻辑回归的理论知识，下一篇博客将介绍逻辑回归在Python中的应用。 参考： http://www.cnblogs.com/sparkwen/p/3441197.html https://tech.meituan.com/intro_to_logistic_regression.html http://blog.csdn.net/bitcarmanlee/article/details/51165444 http://ufldl.stanford.edu/wiki/index.php/Softmax%E5%9B%9E%E5%BD%92

上一篇博客简单介绍了Linux平台下不同性能分析工具的特点，最后是Google出品的gperftools胜出。今天我们将举一个用gperftools对链接库进行性能分析的例子。 假设我们创建了一个TestGperftools项目，其中包括一个动态链接库模块complex和主程序main，build文件夹用于存放动态链接库和可执行程序，文件结构如下所示。 1 2 3 4 5 6 7 8 9 10 11 12 13 czl@ubuntu:~/czl/TestGperftools$ tree . ├── build ├── complex │ ├── complex.cpp │ ├── complex.h │ └── Makefile ├── main │ ├── main.cpp │ └── Makefile └── Makefile 3 directories, 6 files 其中complex借用了这篇博客的第三个公式计算圆周率π，主要用来模拟一个很耗时的操作。该模块会生成一个libcomplex.so动态链接库供外部调用。其complex.cpp函数如下： ...

这半年一直在研究pLink 2的加速策略，了解了相关的性能分析工具，现记录如下。 对软件进行加速的技术路线一般为：首先对代码进行性能分析（ performance analysis，也称为 profiling），然后针对性能瓶颈提出优化方案，最后在大数据集上评测加速效果。所以进行性能优化的前提就是准确地测量代码中各个模块的时间消耗。听起来很简单，不就是测量每行代码的运行时间吗，直接用time_t t=clock();不就好了，但是事情并没有那么简单。如果只进行粗粒度的性能分析，比如测量几个大的模块的运行时间，clock()还比较准确，但是如果测量的是运算量比较小的函数调用，而且有大量的小函数调用，clock()就不太准确了。 比如下面的一段代码，我最开始的性能分析方法是在fun1()~fun3()前后添加time_t t=clock()，然后作差求和的。但是3个fun()加起来的时间居然不等于整个while循环的时间，有将近50%的时间不翼而飞了！ 1 2 3 4 5 6 7 8 9 10 11 while (true) { if (fun1()) { for (int i = 0; i < k; ++k) { if (flag1) { fun2(); } } } else { fun3(); } } 一种可能是while循环以及内部的for循环本身占用了较多的时间，但是好像不太可能呀。还有一种可能是clock()测量有误，time_t只能精确到秒，clock_t只能精确到毫秒，如果一次fun*()的时间太短，导致一次测量几乎为0，那么多次的while和for循环调用累加的时间也几乎为0，导致实际测量到的fun*()时间远小于真实时间。所以自己用代码进行性能分析可能会有较大的误差，最好借助已有的性能分析工具。 ...

今天阅读《C++ Primer, 5e》的第二章，介绍C++的基本内置类型，觉得有一些平时工作容易出错的知识点，现摘录如下： 1 unsigned char c = -1; // 假设char占8比特，c的值为255 当我们赋给无符号类型一个超出它表示范围的值时，结果是初始值对无符号类型表示数值总数取模后的余数。例如，8比特大小的unsigned char可以表示0至255区间内的值，如果我们赋了一个区间以外的值，则实际的结果是该值对256取模后所得的余数。因此，把-1赋给8比特大小的unsigned char所得的结果是255。 1 signed char c2 = 256; // 假设char占8比特，c2的值是未定义的 当我们赋给带符号类型一个超出它表示范围的值时，结果是未定义的（undefined）。此时，程序可能继续工作、可能崩溃，也可能生成垃圾数据。 1 2 3 4 unsigned u = 10; int i = -42; std::cout << i + i << std::endl; // 输出-84 std::cout << u + i << std::endl; // 如果int占32位，输出4294967264 在第一个输出表达式里，两个（负）整数相加并得到了期望的结果。在第二个输出表达式里，相加前首先把整数-42转换成无符号数。把负数转换成无符号数类似于直接给无符号数赋一个负数，结果等于这个负数加上无符号数的模。unsigned (int)的取值范围是0~\(2^{32}-1\)，所以总数有\(2^{32}\)个数，-42%\(2^{32}\)=-42+\(2^{32}\)，u+i=10+(-42+\(2^{32}\))=4294967264。 1 2 3 unsigned u1 = 42, u2 = 10; std::cout << u1 – u2 << std::endl; // 正确：输出32 std::cout << u2 – u1 << std::endl; // 正确：不过，结果是取模后的值 当从无符号数中减去一个值时，不管这个值是不是无符号数，我们都必须确保结果不能是一个负值。 ...

随机矩阵是这样一类方阵，其元素为非负实数，且行和或列和为1。如果行和为1，则称为行随机矩阵；如果列和为1，则称为列随机矩阵；如果行和和列和都为1，则称为双随机矩阵。 前面我们介绍的谷歌矩阵和HMM中的转移矩阵都属于随机矩阵，所以随机矩阵也称为概率矩阵、转移矩阵、或马尔可夫矩阵。 随机矩阵有一个性质，就是其所有特征值的绝对值小于等于1，且其最大特征值为1。下面通过两种方法证明这个结论。 首先，随机矩阵A肯定有特征值1，即 $$\begin{equation}A\vec 1=1\times\vec 1\end{equation}$$其中的单位向量\(\vec 1=(\frac{1}{n},…,\frac{1}{n})^T\)，因为A的行和为1，所以上述等式成立。即1是A的特征值。 反证法 假设存在大于1的特征值\(\lambda\)，则有\(A\vec x=\lambda\vec x\)。令\(x_k\)是\(\vec x\)中最大的元素。又因为A的元素非负，且行和为1，所以\(\lambda\vec x\)中的每个元素都是\(\vec x\)中元素的凸组合，所以\(\lambda\vec x\)中的每个元素都小于等于\(x_k\)。 $$\begin{equation}a_{i1}x_1+a_{i2}x_2+…+a_{in}x_n=\lambda x_i\leq x_k\end{equation}$$但是如果\(\lambda>1\)，则\(\lambda x_k>x_k\)，和(2)式矛盾，所以\(\lambda\leq 1\)。又因为(1)式，所以A的最大特征值为1。 常规证法 设对称随机矩阵A的特征值\(\lambda\)对应的特征向量为\(x\)（为了简便，以下省略向量符号），则有\(Ax=\lambda x\)，即\(x^TAx=\lambda x^Tx\)，欲证明\(|\lambda|\leq 1\)，只需证明 $$\begin{equation}\lambda=\frac{< x, Ax >}{< x, x >}\leq 1\end{equation}$$根据定义有： $$\begin{equation}< x, Ax >=\sum_{i=1}^na_{ii}x_i^2+2\sum_{i < j, i\sim j}a_{ij}x_ix_j\end{equation}$$对于\(i < j, i\sim j\)，有： $$\begin{equation}a_{ij}(x_i-x_j)^2=a_{ij}x_i^2-2a_{ij}x_ix_j+a_{ij}x_j^2\end{equation}$$两边求和并移项得到： $$ \begin{equation} \begin{array} \displaystyle{2\sum_{i < j}}a_{ij}x_ix_j & = & \displaystyle{\sum_{i < j}a_{ij}x_i^2+\sum_{i < j}a_{ij}x_j^2-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_{i < j}a_{ij}x_i^2+\sum_{i < j}a_{ji}x_j^2-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_{i < j}a_{ij}x_i^2+\sum_{i > j}a_{ij}x_i^2-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_i(\sum_{j\neq i}a_{ij}x_i^2)-\sum_{i < j}a_{ij}(x_i-x_j)^2}\\ & = & \displaystyle{\sum_i(x_i^2(1-a_{ii}))-\sum_{i < j}a_{ij}(x_i-x_j)^2} \end{array} \end{equation} $$第2、3个等号都是因为A是对称矩阵，所以可以把\(a_{ij}\)替换为\(a_{ji}\)，然后互换\(i,j\)下标。最后一个等号是因为A的行和为1。 ...

马尔可夫聚类算法（The Markov Cluster Algorithm, MCL）是一种快速可扩展的基于图的聚类算法。它的基本思想为：在一个稀疏图G中，如果某个区域A是稠密的（是一个聚类），则在A中随机游走k步，还在A内的概率很大，也就是说，A内的k步路径（k-length path）很多。所以我们可以在图中随机游走k步，如果某个区域连通的概率很大，则该区域是一个聚类。随机游走的下一步只和当前所处节点有关，也就是说这是一个马尔可夫的随机游走过程。 我们用一个例子来演示马尔可夫聚类算法的过程。 上图是一个很小的网络，我们用肉眼大概能看出有三个聚类，分别是左边的{1,6,7,10}，中间的{2,3,5}和右边的{4,8,9,11,12}。我们用MCL看看结果如何。 为了随机游走，我们常用邻接矩阵来表示图，如果i,j有边，则N[i][j]=1，否则N[i][j]=0。又随机游走可能有自回路，所以加上单位矩阵I，得到矩阵N+I。 MCL有两个关键的步骤，分别是Expansion和Inflation。 Expansion就是不断对矩阵进行幂次运算，相当于随机游走。假设随机游走了2步，则得到如下图的关联矩阵\((N+I)^2\)，第1行第10列为4，说明1到10的2-length path有4条：1→6→10，1→7→10，1→1→10，1→10→10。随机游走k步之后，\((N+I)^k[i][j]\)越大，说明\(i\)和\(j\)之间的连通性越强。 $$\begin{equation}Expand(M)=M^k\end{equation}$$ Inflation是为了增强更强的连接，减弱更弱的连接，只有这样才能得到边界比较明确的聚类。MCL的做法是对元素做幂次运算，然后按列归一化，公式为： $$\begin{equation}(\Gamma_rM)_{pq}=\frac{(M_{pq})^r}{\sum_{i=1}^k(M_{iq})^r}\end{equation}$$参数经验值是\(k=r=2\)。不断做Expansion和Inflation操作，直到算法收敛，得到若干个聚类。中间过程请点此查看，下图为最终结果。 从图中可以看出，和1有边的只剩下6,7,10了，所以得到聚类{1,6,7,10}，同理能得到聚类{2,3,5}和{4,8,9,11,12} ，和我们肉眼得到的结果是一致的。 MCL算法的原理很简单，得到的聚类效果也不错。下面总结一下MCL的算法过程： 给定无向图G，Expansion和Inflation的参数\(k\)和\(r\) 生成G的邻接矩阵\(N\) 添加自回路，得到矩阵\(N+I\) 循环对\(N+I\)做Expansion和Inflation操作，即计算公式(1)和(2)，直到收敛 根据最终得到的矩阵，进行划分聚类 此算法是我在上《生物信息学中的算法设计》课上是学到的，当时觉得这个算法真是神奇，如此简单，但又如此有效，实在高明。查阅文献得知，此为Stijn van Dongen的博士论文，本博客的图片均来自其博士论文，想深入了解图聚类算法，请下载他的论文。

上一回介绍了HMM的解码问题，今天我们介绍HMM的学习问题和识别问题，先来看学习问题。 正如上一回结束时所说，HMM的学习问题是：仅已知观测序列\(\vec y\)，要估计出模型参数组\(\vec\lambda=(\mu,A,B)\)，其中\(\mu\)为初始概率分布向量，\(A\)为转移概率矩阵，\(B\)为发射概率矩阵。 算法设计 求解HMM的参数学习问题，就是求解如下的最优化问题： $$\begin{equation} P(\vec Y = \vec y|\hat \lambda)=\max\limits_{\vec \lambda} P(\vec Y = \vec y|\vec \lambda)\end{equation}$$也就是找一个参数\(\vec \lambda\)，使得模型在该参数下最有可能产生当前的观测\(\vec y\)。如果使用极大似然法求解，对于似然函数\(P(\vec Y=\vec y|\vec \lambda)=\sum\limits_{i_1,…,i_T}\mu_{i_1}b_{i_1y_1}a_{i_1i_2}…a_{i_{T-1}i_T}b_{i_Ty_T}\)而言，这个最大值问题的计算量过大，在实际中是不可能被采用的。为此，人们构造了一个递推算法，使其能相当合理地给出模型参数\(\vec \lambda\)的粗略估计。其核心思想是：并不要求备选\(\vec\lambda\)使得\(P(\vec Y=\vec y|\vec \lambda)\)达到最大或局部极大，而只要求使\(P(\vec Y=\vec y|\vec \lambda)\)相当大，从而使计算变为实际可能。 EM算法 为此，我们定义一个描述模型“趋势”的量\(Q(\vec\lambda^*|\vec\lambda)\)代替似然函数\(P(\vec Y=\vec y|\vec\lambda)\)，其定义为： $$\begin{equation} Q(\vec\lambda^*|\vec\lambda)=\sum\limits_{\vec x}P(\vec x,\vec y|\vec\lambda)\ln P(\vec x,\vec y|\vec\lambda^*)\end{equation}$$利用在\(0 < x < 1\)时，不等式\(\ln x\leq x-1\)成立，可以证明： $$\begin{equation} Q(\vec\lambda^*|\vec\lambda)-Q(\vec\lambda|\vec\lambda)\leq P(\vec Y=\vec y|\vec\lambda^*)-P(\vec Y=\vec y|\vec\lambda)\end{equation}$$由此可见，对于固定的\(\vec\lambda\)，只要\(Q(\vec\lambda^*|\vec\lambda)>Q(\vec\lambda|\vec\lambda)\)，就有\(P(\vec Y=\vec y|\vec\lambda^*)>P(\vec Y=\vec y|\vec\lambda)\)。于是想把模型\(\vec\lambda_m\)修改为更好的模型\(\vec\lambda_{m+1}\)，只需找\(\vec\lambda_{m+1}\)使得： $$\begin{equation}Q(\vec\lambda_{m+1}|\vec\lambda_m)=\sup_{\vec\lambda}Q(\vec\lambda|\vec\lambda_m)\end{equation}$$即只要把\(Q(\vec\lambda|\vec\lambda_m)\)关于\(\vec\lambda\)的最大值处取成\(\vec\lambda_{m+1}\)，就有\(P(\vec Y=\vec y|\vec\lambda_{m+1})>P(\vec Y=\vec y|\vec\lambda_m)\)。 这样得到的模型序列\(\{\vec\lambda_m\}\)能保证\(P(\vec Y=\vec y|\vec\lambda_m)\)关于\(m\)是严格递增的，虽然在这里还不能在理论上证明\(P(\vec Y=\vec y|\vec\lambda_m)\)收敛到\(\max_{\vec\lambda}P(\vec Y=\vec y|\vec\lambda)\)，但是当\(m\)充分大时，\(\vec\lambda_m\)也还能提供在实际中较为满意的粗略近似。 ...