袋熊的树洞

介绍

ConfigSpace是一个用于管理算法参数空间的Python包，主要用于算法参数选择任务。一些AutoML库，例如SMAC3、BOHB以及auto-sklearn，会用到该包。项目主页为：https://github.com/automl/ConfigSpace。

注明：本文章相关代码在Gist

初始化

使用ConfigSpace包时通常要创建一个参数空间实例

1
2
3
4

import ConfigSpace as CS
import ConfigSpace.hyperparameters as CSH

cs = CS.ConfigurationSpace()

这个参数空间集合实例cs包含所有参数的设置

整数参数和浮点参数

本节开始将介绍如何配置算法的参数空间，这里举例的算法为SVM分类算法，算法具体实现为sklearn.svm.SVC。由SVC类介绍可以知道两个参数：

1. C为惩罚参数，数据类型为浮点数，且$C \ge 0$
2. max_iter为最大迭代次数，数据类型为整数

假设要限定C的取值范围为$[0, 1]$以及max_iter的取值范围为$[10, 100]$，可以用UniformFloatHyperparameter和UniformIntegerHyperparameter设定参数范围

1
2

param_c = CSH.UniformFloatHyperparameter(name='C', lower=0, upper=1)
param_max_iter = CSH.UniformIntegerHyperparameter(name='max_iter', lower=10, upper=100)

设定完参数空间后，需要添加到参数空间集合实例cs中

1
2

cs.add_hyperparameter(param_c)
cs.add_hyperparameter(param_max_iter)

此时可以使用cs的sample_configuration方法进行采样获得一组随机的参数

1

cs.sample_configuration()

此时输出类似下面这种情况

1
2
3

Configuration:
  C, Value: 0.7114185317566737
  max_iter, Value: 84

Categorical参数和参数之间的联系

由sklearn.svm.SVC类介绍可知，算法核类型由参数kernel控制

• kernel限定算法的核类型，取值主要有'linear'，'poly'，'rbf'，'sigmoid'

此时可以用CategoricalHyperparameter来代表参数kernel

1
2
3

param_kernel = CSH.CategoricalHyperparameter(name='kernel', choices=['linear', 'poly', 'rbf', 'sigmoid'])

cs.add_hyperparameter(param_kernel)

每一种核还有相应的参数设置（设定SVC类对应的参数），即

• Linear核$K(x, y)=x^Ty$，无参数
• Poly核$K(x, y)=(\gamma x^Ty + r)^d$，其中参数$\gamma$对应gamma，参数$r$对应coef0，参数$d$对应degree
• RBF核$K(x, y)=\exp(-\gamma \Vert x - y\Vert^2)$，其中参数$\gamma$对应gamma
• Sigmoid核$K(x, y)=\tanh(\gamma x^T y + r)$，其中参数$\gamma$对应gamma，参数$r$对应coef0

首先创建参数degree、coef0以及gamma的参数空间

1
2
3
4
5

param_degree = CSH.UniformIntegerHyperparameter(name='degree', lower=2, upper=4)
param_coef0 = CSH.UniformFloatHyperparameter(name='coef0', lower=0, upper=1)
param_gamma = CSH.UniformFloatHyperparameter(name='gamma', lower=1e-5, upper=1e2)

cs.add_hyperparameters([param_degree, param_coef0, param_gamma])

有前面的描述可以知道不同的核对应不同的参数，也就是说核参数和核类型参数之间是由关联的

• 参数degree关联Poly核
• 参数coef0关联Poly核和Sigmoid核
• 参数gamma关联Poly核、RBF核和Sigmoid核

要想表示这种参数之间的关系，可以使用EqualsCondition以及OrConjunction，即

1
2
3
4
5
6
7
8

cond1 = CS.EqualsCondition(param_degree, param_kernel, 'poly')
cond2 = CS.OrConjunction(CS.EqualsCondition(param_coef0, param_kernel, 'poly'),
                         CS.EqualsCondition(param_coef0, param_kernel, 'sigmoid'))
cond3 = CS.OrConjunction(CS.EqualsCondition(param_gamma, param_kernel, 'rbf'),
                         CS.EqualsCondition(param_gamma, param_kernel, 'poly'),
                         CS.EqualsCondition(param_gamma, param_kernel, 'sigmoid'))

cs.add_conditions([cond1, cond2, cond3])

其中

1

CS.EqualsCondition(param_degree, param_kernel, 'poly')

意思为参数kernel值为'poly'时，设定参数degree的值。如果有多个条件，需要用OrConjunction来OR这些条件

1
2

cond2 = CS.OrConjunction(CS.EqualsCondition(param_coef0, param_kernel, 'poly'),
                         CS.EqualsCondition(param_coef0, param_kernel, 'sigmoid'))

意思为当参数kernel值为'poly'时，设定参数coef0值，或者当参数kernel值为'sigmoid'时，设定参数coef0值。

禁止参数取值组合出现

前面我们设定了sklearn.svm.SVC类某些参数的参数空间，假如SVC的核选择的是Linear核，即参数kernel取值为'Linear'，此时SVM变成了LinearSVM。如果SVC类的LinearSVM实现为sklearn.svm.LinearSVC，这时可以用LinearSVC类参数进一步控制算法的运行过程。
注：这里只是假设一种情况，即SVC类有LinearSVC类的全部参数，真实情况是SVC类并没有LinearSVC类的全部参数。

LinearSVC类部分参数如下

• penalty设置正则项类型，数据类型为字符串，取值为'l1'或者'l2'
• loss设置损失函数类型，数据类型为字符串，取值为'hinge'或者'squared_hinge'
• dual设置算法是否求解对偶问题，数据类型为布尔值，实际可以替换成字符串类型

首先根据这三个参数设置参数空间

1
2
3
4
5

param_penalty = CSH.CategoricalHyperparameter(name='penalty', choices=['l1', 'l2'], default_value='l2')
param_loss = CSH.CategoricalHyperparameter(name='loss', choices=['hinge', 'squared_hinge'], default_value='squared_hinge')
param_dual = CSH.CategoricalHyperparameter(name='dual', choices=['True','False'], default_value='False')

cs.add_hyperparameters([param_penalty, param_loss, param_dual])

当核类型为Linear核时，这三个参数才会被设置，因此要进行参数关联

1
2
3
4
5

cond1 = CS.EqualsCondition(param_penalty, param_kernel, 'linear')
cond2 = CS.EqualsCondition(param_loss, param_kernel, 'linear')
cond3 = CS.EqualsCondition(param_dual, param_kernel, 'linear')

cs.add_conditions([cond1, cond2, cond3])

这里限定一些参数组合不能出现

• 参数penalty取值'l1'，参数loss取值'hinge'
• 参数dual取值'False'，参数penalty取值'l2'，参数loss取值'hinge'
• 参数dual取值'False'，参数'penalty'取值'l1'

要禁止出现某些参数组合，可以使用ForbiddenEqualsClause，如果有多个组合，需要使用ForbiddenAndConjunction进行OR

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

penalty_loss = CS.ForbiddenAndConjunction(
        CS.ForbiddenEqualsClause(param_penalty, 'l2'),
        CS.ForbiddenEqualsClause(param_loss, 'hinge')
)
dual_penalty_loss = CS.ForbiddenAndConjunction(
        CS.ForbiddenEqualsClause(param_dual, 'False'),
        CS.ForbiddenEqualsClause(param_penalty, 'l2'),
        CS.ForbiddenEqualsClause(param_loss, 'hinge')
)

penalty_dual = CS.ForbiddenAndConjunction(
        CS.ForbiddenEqualsClause(param_dual, 'False'),
        CS.ForbiddenEqualsClause(param_penalty, 'l1')
)

cs.add_forbidden_clause(penalty_loss)
cs.add_forbidden_clause(dual_penalty_loss)
cs.add_forbidden_clause(penalty_dual)

安装桌面环境和VNC服务端

首先更新包列表

1

$ sudo apt-get update

安装桌面环境Xfce

1

$ sudo apt-get install xfce4 xfce4-goodies

安装VNC服务端

1

$ sudo apt-get install tightvncserver

设置VNC连接密码设置以及生成配置文件

首先执行vncserver命令来设置VNC连接密码以及生成VNC配置文件

1

$ vncserver

执行命令后会要求设置连接密码，显示以下内容

1
2
3
4

You will require a password to access your desktops.

Password:
Verify:

设置完密码后，命令会生成VNC配置文件并启动一个VNC实例

1
2
3
4
5

New 'X' desktop is your_hostname:1

Creating default startup script /home/your_username/.vnc/xstartup
Starting applications specified in /home/your_username/.vnc/xstartup
Log file is /home/your_username/.vnc/your_hostname:1.log

配置文件在下面目录里面

1

/home/your_username/.vnc/

第一次运行vncserver命令会自动启动VNC实例，分配到:1上，对应端口为5901 (端口5901=5900+1，如果是:2，则端口为5902，以此类推)。由于要配置VNC，所以先要关闭VNC实例

1

$ vncserver -kill :1

关闭成功后会显示以下信息

1

Killing Xtightvnc process ID 30095

配置VNC

要配置的文件为xstartup，该文件在$HOME/.vnc里面，即

1

/home/your_username/.vnc/

首先备份原始配置文件

1

$ mv ~/.vnc/xstartup ~/.vnc/xstartup.bak

然后创建新的配置文件

1

$ touch ~/.vnc/xstartup

编辑该文件，添加以下内容

1
2
3
4
5
6

#!/bin/sh

xrdb $HOME/.Xresources
unset SESSION_MANAGER
unset DBUS_SESSION_BUS_ADDRESS
startxfce4 &

网上很多教程会省略销毁那两个环境变量，即没有

1
2

unset SESSION_MANAGER
unset DBUS_SESSION_BUS_ADDRESS

我实际操作时发现如果没有销毁变量，连接上VNC时，画面是全灰的，销毁变量后，显示就正常了。为了保证VNC配置文件能够生效，赋予该文件执行权限

1

$ chmod +x ~/.vnc/xstartup

VNC连接

启动VNC实例，即执行命令

1

$ vncserver

执行成功后，显示

1
2
3
4

New 'X' desktop is your_hostname:1

Starting applications specified in /home/your_username/.vnc/xstartup
Log file is /home/your_username/.vnc/your_hostname:1.log

查看端口开启情况，可以看到5901端口已经开启了

1

$ ss -ltn

此时可以用VNC连接该电脑了，连接地址格式为ip:port

参考

[1] How to Install and Configure VNC on Ubuntu 18.04

[2] VNC server on Ubuntu 18.04 Bionic Beaver Linux

[3] 2018-06-07【Ubuntu 18.04 搭建VNC服务器】

介绍

有些机器学习算法，比如最小二乘支持向量机 (Least Squares Support Vector Machine, LSSVM)，需要计算核矩阵 (Kernel Matrix)。核矩阵的定义如下

$$
\Omega = K(X, Y) = \begin{bmatrix}
K(x_1, y_1) & K(x_1, y_2) & \cdots & K(x_1, y_n) \
K(x_2, y_1) & K(x_2, y_2) & \cdots & K(x_2, y_n) \
\vdots & \vdots & \ddots & \vdots \
K(x_m, y_1) & K(x_n, y_2) & \cdots & K(x_m, y_n)
\end{bmatrix}
$$

其中$K(\cdot, \cdot)$为核函数，常见的有线性核以及RBF核，矩阵$X \in R^{m \times d}$和$Y \in R^{n \times d}$为数据矩阵，每一行代表一个样本，定义为

$$
\begin{aligned}
X & = [x_1 , x_2 , \cdots , x_m]^T \
Y & = [y_1 , y_2 , \cdots , y_n]^T
\end{aligned}
$$

核矩阵的规模为$m \times n$，可以看出核矩阵的规模是由数据个数决定的，如果按照定义一个个生成核矩阵每一个元素，可想而知生成效率会极其低下，因此考虑用其他方式快速生成核矩阵。

快速生成核矩阵

本文只考虑核函数为RBF核的情况，其他情况有相类似的推导。RBF核又称高斯核，其定义为

$$
K(x, y) = e^{-\frac{\Vert x - y \Vert^2}{2\sigma^2}}
$$

其中$\sigma > 0$为高斯核带宽，$\Vert\cdot\Vert$为$L_2$范数。此时核矩阵$\Omega$为

$$
\Omega=K(X, Y) = \begin{bmatrix}
e^{-\frac{\Vert x_1 - y_1\Vert^2}{2\sigma^2}} & e^{-\frac{\Vert x_1 - y_2 \Vert^2}{2\sigma^2}} & \cdots & e^{-\frac{\Vert x_1 - y_n \Vert^2}{2\sigma^2}} \
e^{-\frac{\Vert x_2 - y_1 \Vert^2}{2\sigma^2}} & e^{-\frac{\Vert x_2 - y_2 \Vert^2}{2\sigma^2}} & \cdots & e^{-\frac{\Vert x_2 - y_n \Vert^2}{2\sigma^2}} \
\vdots & \vdots & \ddots & \vdots \
e^{-\frac{\Vert x_m - y_1 \Vert^2}{2\sigma^2}} & e^{-\frac{\Vert x_m - y_2 \Vert^2}{2\sigma^2}} & \cdots & e^{-\frac{\Vert x_m - y_n \Vert^2}{2\sigma^2}}
\end{bmatrix}
$$

将自然指数运算提取到矩阵外面，则核矩阵$\Omega$表示为$\Omega=e^{-\frac{D}{2\sigma^2}}$，这里矩阵$D$为

$$
D = \begin{bmatrix}
\Vert x_1 - y_1 \Vert^2 & \Vert x_1 - y_2 \Vert^2 & \cdots & \Vert x_1 - y_n \Vert^2 \
\Vert x_2 - y_1 \Vert^2 & \Vert x_2 - y_2 \Vert^2 & \cdots & \Vert x_2 - y_n \Vert^2 \
\vdots & \vdots & \ddots & \vdots \
\Vert x_m - y_1 \Vert^2 & \Vert x_m - y_2 \Vert^2 & \cdots & \Vert x_m - y_n \Vert^2
\end{bmatrix}
$$

将矩阵$D$中每一个元素展开，则可得

$$
\begin{aligned}
\Vert x_i - y_j \Vert^2 & = (x_i - y_j)^T(x_i - y_j) \
& = x_i^T x_i - 2 x_i^T y_j + y_j^T y_j \
& = \Vert x_i \Vert^2 - 2 x_i^T y_j + \Vert y_j \Vert^2
\end{aligned}
$$

从而矩阵$D$可以分解为三个矩阵之和

$$
D = A - 2B + C
$$

其中矩阵$A$定义为

$$
A = \begin{bmatrix}
\Vert x_1 \Vert^2 & \Vert x_1 \Vert^2 & \cdots & \Vert x_1 \Vert^2 \
\Vert x_2 \Vert^2 & \Vert x_2 \Vert^2 & \cdots & \Vert x_2 \Vert^2 \
\vdots & \vdots & \ddots & \vdots \
\Vert x_m \Vert^2 & \Vert x_m \Vert^2 & \cdots & \Vert x_m \Vert^2
\end{bmatrix}_{m \times n}
$$

这里定义一个向量$\alpha$，其第$i$个元素为数据矩阵$X$第$i$行所有元素平方和，即

$$
\alpha = \begin{bmatrix}
\Vert x_1 \Vert^2 & \Vert x_2 \Vert^2 & \cdots & \Vert x_m \Vert^2
\end{bmatrix}^T
$$

则矩阵$A$可以表示为

$$
\begin{aligned}
A & = \alpha 1^T \
& = \begin{bmatrix}
\Vert x_1 \Vert^2 \
\Vert x_2 \Vert^2 \
\vdots \
\Vert x_m \Vert^2
\end{bmatrix} \begin{bmatrix}
1 & 1 & \cdots & 1
\end{bmatrix}
\end{aligned}
$$

矩阵$B$定义为

$$
B = \begin{bmatrix}
x_1^T y_1 & x_1^T y_2 & \cdots & x_1^T y_n \
x_2^T y_1 & x_2^T y_2 & \cdots & x_2^T y_n \
\vdots & \vdots & \ddots & \vdots \
x_m^T y_1 & x_m^T y_2 & \cdots & x_m^T y_n
\end{bmatrix}_{m \times n}
$$

很容易可以知道矩阵$B$由数据矩阵$X$和$Y$相乘而得

$$
B = XY^T
$$

矩阵$C$定义为

$$
C = \begin{bmatrix}
\Vert y_1 \Vert^2 & \Vert y_2 \Vert^2 & \cdots & \Vert y_n \Vert^2 \
\Vert y_1 \Vert^2 & \Vert y_2 \Vert^2 & \cdots & \Vert y_n \Vert^2 \
\vdots & \vdots & \ddots & \vdots \
\Vert y_1 \Vert^2 & \Vert y_2 \Vert^2 & \cdots & \Vert y_n \Vert^2
\end{bmatrix}_{m \times n}
$$

类似地，这里定义一个向量$\beta$，其第$i$个元素为数据矩阵$Y$第$i$行所有元素平方和，即

$$
\beta = \begin{bmatrix}
\Vert y_1 \Vert^2 & \Vert y_2 \Vert & \cdots & \Vert y_n \Vert^2
\end{bmatrix}^T
$$

则矩阵$C$可以表示为

$$
\begin{aligned}
C & = 1 \beta^T \
& = \begin{bmatrix}
1 \
1 \
\vdots \
1
\end{bmatrix} \begin{bmatrix}
\Vert y_1 \Vert^2 & \Vert y_2 \Vert & \cdots & \Vert y_n \Vert^2
\end{bmatrix}
\end{aligned}
$$

综上可以得到

$$
D = \alpha 1^T - 2 XY^T + 1 \beta^T
$$

实现

原始生成方法

按照核矩阵定义，一个个生成矩阵每一个元素，实现代码如下

1
2
3
4
5
6
7
8
9
10
11

def method1(X1, X2, sigma=1.0):
    num_samples1, dims = X1.shape
    num_samples2, dims = X2.shape
    mat = np.empty((num_samples1, num_samples2))
    
    for i in range(num_samples1):
        for j in range(num_samples2):
            vec = X1[i, :] - X2[j, :]
            mat[i, j] = np.exp(-0.5 / sigma**2 * np.dot(vec, vec))
    
    return mat

快速生成方法1

根据前面讨论，将核矩阵分解为一系列矩阵向量乘积，实现代码如下

1
2
3
4
5
6
7
8
9
10
11

def method2(X1, X2, sigma=1.0):
    X1_sum = np.sum(X1**2, axis=1)
    X2_sum = np.sum(X2**2, axis=1)
    mat1 = np.dot(np.reshape(X1_sum, (X1_sum.shape[0], 1)),
                  np.ones((1, X2_sum.shape[0])))
    mat2 = np.dot(np.ones((X1_sum.shape[0], 1)),
                  np.reshape(X2_sum, (1, X2_sum.shape[0])))
    mat = mat1 + mat2 - 2 * np.dot(X1, X2.T)
    mat = np.exp(-0.5 / sigma**2 * mat)
    
    return mat

快速生成方法2

利用Numpy的广播机制，进一步提升核矩阵生成速度，实现代码如下

1
2
3
4
5

def method3(X1, X2, sigma=1.0):
    mat = np.sum(X1**2, 1).reshape(-1, 1) + np.sum(X2**2, 1) - 2 * np.dot(X1, X2.T)
    mat = np.exp(-0.5 / sigma**2 * mat)
    
    return mat

实验

实验代码在：https://gist.github.com/luowanqian/d2f655a95f77dfe088fd626f72d78966

这里做一个简单的实验，比对各种方法的效率。实验的数据矩阵$X$大小为$1000 \times 100$，数据矩阵$Y$大小为$5000 \times 100$，实验结果如下

由结果可以看出，快速生成方法速度提升是极大的。

问题

使用Tmux时经常会遇到的一个问题，那就是系统重启时，Session会被清除，每次打开电脑都要重启建一个Session，然后创建一堆Window以及Pane，这极大地降低了Tmux使用效率，因此需要想一个办法能够保存Tmux的Session。

解决方案

该问题的解决方案是安装一个叫Tmux Resurrect的插件。Tmux要安装插件，可以通过Tmux Plugin Manager这个插件进行安装，该插件相当于一个插件管理系统，可以快速地安装、更新以及删除插件。

Tmux Plugin Manager

安装Tmux Plugin Manager插件可以参考该插件的GitHub：tmux-plugins/tpm。安装很简单，首先clone插件到本地

1

$ git clone https://github.com/tmux-plugins/tpm ~/.tmux/plugins/tpm

然后修改.tmux.conf，在文件最底部添加以下内容

1
2
3
4
5
6

# List of plugins
set -g @plugin 'tmux-plugins/tpm'
set -g @plugin 'tmux-plugins/tmux-sensible'

# Initialize TMUX plugin manager (keep this line at the very bottom of tmux.conf)
run '~/.tmux/plugins/tpm/tpm'

重新加载配置文件

1
2

# type this in terminal if tmux is already running
$ tmux source ~/.tmux.conf

然后就可以在Tmux中使用快捷prefix + I(注意这里的I是大写)安装配置文件.tmux.conf中定义的插件了。

补充：

如果Tmux安装了gpakosz/.tmux，在文件.tmux.conf.local中配置好Tmux Plugin Manager后，在Tmux中使用prefix + I安装插件会没有效果，相关的讨论见“run ‘~/.tmux/plugins/tpm/tpm’” has no effect in .tmux.conf.local #61，具体的解决方案是用另一种写法写set -g @plugin，即将

1
2

set -g @plugin 'tmux-plugins/tpm'
set -g @plugin 'tmux-plugins/tmux-sensible'

改写为

1
2
3
4

set -g @tpm_plugins '          \
  tmux-plugins/tpm             \
  tmux-plugins/tmux-sensible
'

解决方案来源于：Help, tpm not working!

Tmux Resurrect

安装好Tmux Plugin Manager后，就可以安装Tmux Resurrect插件(tmux-plugins/tmux-resurrect)了。在Tmux配置文件.tmux.conf中添加

1

set -g @plugin 'tmux-plugins/tmux-resurrect'

然后使用快捷键prefix + I就可以安装插件了。

保存和恢复

安装完Tmux Plugin Manager，就可以使用下面快捷键保存和恢复Session了。

• prefix + Ctrl-s - save
• prefix + Ctrl-r - restore

平常重启电脑后，首先在终端打开Tmux，然后使用prefix + Ctrl-r就可以恢复保存的Session了，如果有多个Session可以使用prefix + s来选择Session。

参考

1. 保存和恢复 Tmux 会话

介绍

看本文之前建议看文章”Typed Memoryviews“，该文章主要介绍的是Memoryviews，使用Memoryviews可以直接访问Numpy数组的memory buffer，直接操作Numpy数组的数据，同时是连接C数组与Numpy数组之间的桥梁。相关代码在：CythonDemos/Demo1

BLAS

BLAS，全称Basic Linear Algebra Subprograms，即基础线性代数子程序库，里面有大量已经编写好的关于线性代数运算的程序。BLAS库的API分三大类，分别为Level 1、Level 2以及Level 3。

• Level 1: 函数处理单一向量的线性运算以及两个向量的二元运算
• Level 2: 函数处理矩阵与向量的运算，同时也包含线性方程组求解
• Level 3: 函数处理矩阵与矩阵之间的运算

API的全部介绍可以参考官方文档：blasqr.pdf

Cython调用BLAS函数

Cython要调用BLAS函数，主要是通过BLAS库的C接口，接口的具体定义可以参考现有BLAS库的 cblas.h，里面函数的命名和 blasqr.pdf 中函数命名是一一对应的，只不过函数前多了cblas_前缀。这里演示如何调用函数cblas_dgemv，在头文件 cblas.h 中函数声明为

1

void cblas_dgemv(OPENBLAS_CONST enum CBLAS_ORDER order, OPENBLAS_CONST enum CBLAS_TRANSPOSE trans, OPENBLAS_CONST blasint m, OPENBLAS_CONST blasint n, OPENBLAS_CONST double alpha, OPENBLAS_CONST double *a, OPENBLAS_CONST blasint lda, OPENBLAS_CONST double *x, OPENBLAS_CONST blasint incx, OPENBLAS_CONST double beta, double *y, OPENBLAS_CONST blasint incy);

去除一些宏以及替换定义的数据类型，简单表示就是

1

void cblas_dgemv(enum CBLAS_ORDER order, enum CBLAS_TRANSPOSE trans, int m, int n, double alpha, double *a, int lda, double *x, int incx, double beta, double *y, int incy)

这个函数主要是实现的功能计算矩阵与向量相乘，在 blasqr.pdf 中对应的是Level2中xGEMV函数

$$
y = \alpha Ax + \beta y
$$

在Cython调用BLAS函数首先要使用cdef extern引入头文件，然后声明函数

1
2
3
4
5
6
7
8
9
10
11
12
13

cdef extern from 'cblas.h':
    enum CBLAS_ORDER:
        CblasRowMajor=101
        CblasColMajor=102
    enum CBLAS_TRANSPOSE:
        CblasNoTrans=111
        CblasTrans=112
        CblasConjTrans=113
        CblasConjNoTrans=114
    void dgemv "cblas_dgemv"(CBLAS_ORDER Order, CBLAS_TRANSPOSE TransA,
                             int M, int N, double alpha, double *A, int lda,
                             double *X, int incX, double beta,
                             double *Y, int incY) nogil

关于enum中CblasRowMajor这些标识符的取值，是按照 cblas.h 中定义来取值。函数声明后就可以编写一个wrapper函数来调用BLAS函数，之所以要写一个wrapper函数，是因为声明的BLAS函数只能在C中进行调用，我们还要使得其能够在Python中进行调用，wrapper函数定义如下：

1
2
3
4
5
6
7
8

cpdef blas_dgemv(double[:, ::1] A, double[::1] x, double[::1] y,
                 int M, int N, double alpha, double beta):
    cdef double *A_ptr = &A[0, 0]
    cdef double *x_ptr = &x[0]
    cdef double *y_ptr = &y[0]

    dgemv(CblasRowMajor, CblasNoTrans, M, N,
          alpha, A_ptr, N, x_ptr, 1, beta, y_ptr, 1)

通过cpdef可以使得定义的函数可以被Python调用。这里说一下如何将MemoryViews转化成C的指针，对于向量x可以通过对第一个元素进行取址来进行转换&x[0]，对于矩阵A则是对其第一行第一列的元素进行取址操作&A[0, 0]。

由于BLAS库中将矩阵数据看成是一个一维数组，所以在调用BLAS函数时需要指定矩阵是Row Major还是Column Major，即参数CBLAS_ORDER Order，这里设定是Row Major，因此函数blas_dgemv中参数A限定是Row Major，即声明为double[:, ::1]。

整个pyx文件的内容为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

cdef extern from 'cblas.h':
    enum CBLAS_ORDER:
        CblasRowMajor=101
        CblasColMajor=102
    enum CBLAS_TRANSPOSE:
        CblasNoTrans=111
        CblasTrans=112
        CblasConjTrans=113
        CblasConjNoTrans=114
    void dgemv "cblas_dgemv"(CBLAS_ORDER Order, CBLAS_TRANSPOSE TransA,
                             int M, int N, double alpha, double *A, int lda,
                             double *X, int incX, double beta,
                             double *Y, int incY) nogil

cpdef blas_dgemv(double[:, ::1] A, double[::1] x, double[::1] y,
                 int M, int N, double alpha, double beta):
    cdef double *A_ptr = &A[0, 0]
    cdef double *x_ptr = &x[0]
    cdef double *y_ptr = &y[0]

    dgemv(CblasRowMajor, CblasNoTrans, M, N,
          alpha, A_ptr, N, x_ptr, 1, beta, y_ptr, 1)

剩下就是将pyx进行编译，然后就可以在Python中调用blas_dgemv函数了。

参考

1. Calling BLAS from Cython

1. Pandas备忘录系列文章

1. Pandas备忘录

2. 介绍

Note: 如果没有特别说明，pd指的是pandas，np指的是numpy。

1
2

import pandas as pd
import numpy as np

3. 查找

3.1. contains函数

主要是用Series.str.contains这个函数，如果参数regex设为True，则查找字符串解析为正则表达式，设为False则用常规查找方式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

>>> s1 = pd.Series(['Mouse', 'dog', 'house and parrot', '23', 'frog'])
>>> s1
0               Mouse
1                 dog
2    house and parrot
3                  23
4                frog
dtype: object
>>> s1.str.contains('og', regex=False)
0    False
1     True
2    False
3    False
4     True
dtype: bool
>>> s1[s1.str.contains('og', regex=False)]
1     dog
4    frog
dtype: object

4. DataFrame增加行

4.1. append函数

可以使用append函数给DataFrame增加行

1
2
3
4
5
6
7
8
9
10
11
12
13
14

>>> df = pd.DataFrame([[1, 2], [3, 4]], columns=list('AB'))
>>> df
   A  B
0  1  2
1  3  4
>>> df2 = pd.DataFrame([[5, 6]], columns=list('AB'))
>>> df2
   A  B
0  5  6
>>> df.append(df2, ignore_index=True)
   A  B
0  1  2
1  3  4
2  5  6

4.2. 使用特定索引名称增加行

使用DataFrame的loc函数可以指定索引名称增加行，前提是得知道列的排列顺序，因此在创建DataFrame时要用columns参数指定列的排列顺序

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

>>> employees = pd.DataFrame(
...     data={'Name': ['John Doe', 'William Spark'],
...           'Age': [23, 24]},
...     index=['Emp001', 'Emp002'],
...     columns=['Name', 'Age'])
>>> employees
                 Name  Age
Emp001       John Doe   23
Emp002  William Spark   24
>>> employees.loc['Emp003'] = ['Sunny', 45]
>>> employees
                 Name  Age
Emp001       John Doe   23
Emp002  William Spark   24
Emp003          Sunny   45

Weighted random sampling with replacement，中文翻译为带权重的有放回地采样，简单来说就是每一个样本都有一个权重，然后根据样本权重大小有放回地采样数据，样本的权重越大，样本被选择的概率也就越大。该问题换种方式描述就是，给定一个权重数组 $w$，大小为 $N$，然后根据权重大小，在下标数组 ${0, 1, \ldots, N-1}$ 中进行有放回地采样，下标 $i$ 被选中的概率为 $w[i]/\text{sum}(w)$。

这个采样方法一个应用是在Boosting算法中，Boosting算法每次训练完模型时，根据模型分类的情况改变数据样本的权重，然后用带权重的样本训练下一个模型，如果模型的训练不支持带权重的样本，则进行带权重的样本采样。

本文主要将博文 WEIGHTED RANDOM SAMPLING WITH REPLACEMENT WITH DYNAMIC WEIGHTS 内容进行解析 (文章拷贝：Copy)。算法实现代码在：adefazio/sampler

Rejection Sampling

这个思想很简单，给定一个权重数组 $w \in R^N$，权重最大值为 $w_\max$，然后在下标数组 ${ 0, 1, \ldots, N-1 }$ 中随机采样得到下标 $i$，接下来就是判断该下标是否被accept，首先根据分布 $X \sim U(0, 1)$ 获得一个随机数 $x$，然后判断 $w_\max * x \le w[i]$，如果符合该条件，则accept该下标 $i$，否则就reject该下标。

算法实现的代码如下：

1
2
3
4
5
6
7
8
9
10

w = [1,4,2,5] # Some data
w_max = max(w)
n = len(w)

while True:
	idx = random.randrange(n)
	u = w_max*random.random()
	if u <= w[idx]:
		break
print idx

算法思想可以用下面图表示：

橘色区域代表accept区域，灰色区域代表reject区域。这个算法有个缺点就是当数据不平衡时，性能会很低，因为如果小部分数据权重很大时，$w_\max$ 会变的很大，数据权重小的下标很难被accept，而且权重小的占大部分，这导致算法会经过很多次循环。

Level Rejection Sampling

这个是在Rejection Sampling之前对下标进行分层 (level)，根据权重 $w[i]$ 的大小将下标 $i$ 分到不同层中，第 $k$ 层权重的范围为 $[2^k, 2^{k+1}]$，由于每一层内的权重的差距不会偏大，对该层进行Rejection Sampling效率会较高。

该方法的思想如下图所示：

如图所示，总共分了5层，每一层存储的是数组下标。下标分层后，在进行采样时，首先随机选择某一层，然后进行Rejection Sampling获得最终想要的下标。

模拟实验

贴出一段实验代码：

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

import numpy as np
from sampler import Sampler
import matplotlib.pyplot as plt


if __name__ == '__main__':
    weights = np.array([1, 1, 3, 5, 2], dtype='d')
    num_entries = len(weights)
    normalized_weights = weights / np.sum(weights)

    weight_sampler = Sampler(num_entries, max_value=100, min_value=1)
    for i in range(num_entries):
        weight_sampler.add(i, weights[i])

    num_samples = 10000
    distro = np.zeros(num_entries)

    for i in range(num_samples):
        idx = weight_sampler.sample()
        distro[idx] += 1

    normalized_distro = distro / np.sum(distro)

    print('Sample distribution: {}'.format(distro))
    print('Weights distribution: {}'.format(weights))
    print('Normalized sample distribution: {}'.format(normalized_distro))
    print('Normalized weights distribution: {}'.format(normalized_weights))

    # plot frequency
    plt.bar(np.arange(num_entries),
            normalized_weights, width=0.35,
            label='Normalized Weights Distribution')
    plt.bar(np.arange(num_entries)+0.35,
            normalized_distro, width=0.35,
            label='Normalized Sample Distribution')
    plt.legend()
    plt.show()

输出为：

1
2
3
4

Sample distribution: [ 843.  824. 2513. 4182. 1638.]
Weights distribution: [1. 1. 3. 5. 2.]
Normalized sample distribution: [0.0843 0.0824 0.2513 0.4182 0.1638]
Normalized weights distribution: [0.08333333 0.08333333 0.25 0.41666667 0.16666667]

得到的结果如下图所示：

图中横坐标代表数组下标，蓝色柱代表各个下标的权重大小 (归一化)，橘色柱代表经过大量采样时，各个下标被采中的比例。从图中可以看出，采样结果基本符合权重分布规律，权重越大，被采中的概率也就越大。

在使用h5py这类库读取HDF5里面的数据时，数据的存储方式为row-major order，而一些机器学习算法的实现在使用数据时需要的存储方式为column-major order，此时需要对数据进行转换。如果数据量很小，转换是没有问题的，如果数据量巨大，进行转换时，由于需要多一倍内存来存放转换的数据，因此可能导致内存溢出。本篇博文主要是重新实现h5py库中的 read_direct 函数，使得该函数在读取HDF5文件时，可以使得读取到内存的数据是以column-major order方式存储。

相关实现代码在: HDF5 gist

思路

基本思路就是一行行地读取HDF5文件里面的数据，然后存储到内存的数组中，内存的数组是以column-major order方式存储数据。这里假设HDF5文件里面的数据存储是连续的，没有进行分块，复杂的情况等以后再实现，这里只是应对最简单的情况。

实现

实现主要是 Python + C 形式，Python部分主要提供一个函数接口用于调用，C部分主要是使用HDF5 C API对HDF5数据进行读取，两者之间的相互调用使用的是Cython。

C

这部分需要使用HDF5库提供的C API对HDF5文件数据进行读取。主要使用的API是 H5Sselect_hyperslab，这个函数可以设定读取存储空间的哪一部分数据或者将数据存在存储空间的哪一部分，这里的存储空间是一个抽象概念，既可以表示HDF5文件里面内容，也可以表示内存申请的存储空间，这里我们分别对HDF5文件内容和内存申请的空间分别进行设定。实现代码如下：

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

/*
 * 从H5文件中读取矩阵数据，以column order存储数据
 *
 * filename: H5文件名
 * dataset_name: 数据集名
 * buffer: 数据内存存储处，大小为 num_rows * num_cols
 * num_rows: 待读取的矩阵的行数
 * num_cols: 待读取的矩阵的列数
 * offset_x: 读取的数据在H5文件中的矩阵的x方向上的偏移 (第几行)
 * offset_y: 读取的数据在H5文件中的矩阵的y方向上的偏移 (第几列)
 *
 */
int read_direct_f(char *file_name, char *dataset_name, double *buffer,
                  int num_rows, int num_cols, int offset_x, int offset_y)
{
    int i;

    hid_t file_id, dataset_id;
    hid_t dataspace_id, memspace_id;
    herr_t status;

    hsize_t dims_mem[1];

    hsize_t offset[2];
    hsize_t count[2];
    hsize_t stride[2];
    hsize_t block[2];
    hsize_t offset_mem[1];
    hsize_t count_mem[1];
    hsize_t stride_mem[1];
    hsize_t block_mem[1];

    // open the hdf5 file and the dataset
    file_id = H5Fopen(file_name, H5F_ACC_RDONLY, H5P_DEFAULT);
    dataset_id = H5Dopen2(file_id, dataset_name, H5P_DEFAULT);

    // get file dataspace
    dataspace_id = H5Dget_space(dataset_id);

    // create memory space
    dims_mem[0] = num_rows;         // avoid numeric overflow
    dims_mem[0] = dims_mem[0] * num_cols;
    memspace_id = H5Screate_simple(1, dims_mem, NULL);

    // specify size and shape of subset to read
    offset[0] = offset_x;
    offset[1] = offset_y;
    count[0] = 1;
    count[1] = num_cols;
    stride[0] = 1;
    stride[1] = 1;
    block[0] = 1;
    block[1] = 1;
    offset_mem[0] = 0;
    count_mem[0] = num_cols;
    stride_mem[0] = num_rows;
    block_mem[0] = 1;

    for (i=0; i<num_rows; i++) {
        offset[0] = offset_x + i;
        offset_mem[0] = i;

        // select subset from dataspace
        status = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET,
                                     offset, stride, count, block);
        if (status < 0) {
            H5Eprint2(H5E_DEFAULT, NULL);
            return -1;
        }

        status = H5Sselect_hyperslab(memspace_id, H5S_SELECT_SET,
                                     offset_mem, stride_mem,
                                     count_mem, block_mem);
        if (status < 0) {
            H5Eprint2(H5E_DEFAULT, NULL);
            return -1;
        }

        // read data
        status = H5Dread(dataset_id, H5T_NATIVE_DOUBLE,
                         memspace_id, dataspace_id,
                         H5P_DEFAULT, buffer);
        if (status < 0) {
            H5Eprint2(H5E_DEFAULT, NULL);
            return -1;
        }
    }

    // close the resources
    H5Sclose(memspace_id);
    H5Sclose(dataspace_id);
    H5Dclose(dataset_id);
    H5Fclose(file_id);

    return 0;
}

基本实现思路就是每次就只读取HDF5文件中矩阵的一行，然后写入到内存空间中，每次读取都要用函数 H5Sselect_hyperslab 函数设定读取区域，同样，每次写入都要使用函数 H5Sselect_hyperslab 设定写入区域。这个实现思想参考了HDF5官方教程 Reading From or Writing To a Subset of a Dataset，教程中仅仅设定了文件的读取区域，没有设定写入区域。

Python

这部分使用Cython来实现Python中调用前面所写的C语言代码，这部分代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13

def read_direct_fortran(str file_name, str dataset_name,
                        double[::1, :] data, int offset_x, int offset_y):
    cdef int num_rows
    cdef int num_cols
    cdef bytes file_name_bytes = bytes(file_name, encoding='utf-8')
    cdef bytes dataset_name_bytes = bytes(dataset_name, encoding='utf-8')
    cdef int ret;

    num_rows = data.shape[0]
    num_cols = data.shape[1]

    ret = read_direct_f(file_name_bytes, dataset_name_bytes, &data[0, 0],
                        num_rows, num_cols, offset_x, offset_y)

问题描述

Given a 32-bit signed integer, reverse digits of an integer.

Example 1:

1
2

Input: 123
Output: 321

Example 2:

1
2

Input: -123
Output: -321

Example 3:

1
2

Input: 120
Output: 21

Note:

Assume we are dealing with an environment which could only store integers within the 32-bit signed integer range: $[−2^{31}, 2^{31} − 1]$. For the purpose of this problem, assume that your function returns 0 when the reversed integer overflows.

Related Topics: Math

原问题: 7. Reverse Integer

中文翻译版: 7. 反转整数

解决方案

解决思路就是：利用求模得到数字的每一位，求位的顺序要从低位到高位，然后根据获得的数位构造出逆数。这里有个困难就在求解过程中会发生数字溢出，一个解决方案就是用更多字节的整数类型来表示逆数，题目中使用的整数int是四个字节，此时可以用八字节的long long代替，然后判断数是否在32位整数表示的范围即可。代码如下：

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

#include <iostream>
using namespace std;

class Solution {
public:
    int reverse(int x) {
        int mod;
        long long num;
        bool postive;

        num = 0;
        postive = x >=0 ? true : false;
        while (x) {
            mod = x % 10;
            x = x / 10;
            num = num * 10 + mod;
            if ((postive && num > 0x7FFFFFFF)
                || (!postive && num < (signed int)0x80000000))
                return 0;
        }

        return num;
    }
};

int main()
{
    int x = 120;

    Solution solu;
    cout << "Revert: " << solu.reverse(x) << endl;

    return 0;
}

介绍

斯特林公式 (Stirling’s approximation) 是一条用来取n的阶乘的近似值的数学公式。常见的形式为

$$
n! \approx \sqrt{2 \pi n} (\frac{n}{e})^n
$$

对n的阶乘取自然对数，可以得到

$$
\ln{(n!)} \approx \frac{1}{2} \ln{(2 \pi)} + (n + \frac{1}{2}) \ln{(n)} - n
$$

一个会经常遇到的问题，那就是计算 $n! / a_n$，单独计算 $n!$ 需要很大的计算量，如果将问题转换为

$$
\frac{n!}{a_n} = e^{\ln{(n!)} - \ln{(a_n)}}
$$

$\ln{(n!)}$ 可以用前面介绍的近似公式进行求解，此时计算量就会小很多。

例子

这里举一个例子来感受下这个近似公式是如何缩减大量计算量。一个经典的概率问题－Birthday Problem

Birthday Prolem

Calculate the probability p that at least two people in a group of k people will have the same birthday, that is, will have been born on the same day of the same month but not necessarily in the same year.

可以很快地想到

至少两个人的生日相同的概率 = 1 - 所有人的生日不相同的概率

则计算公式为

$$
p = 1 - \frac{P_{365, k}}{365^k}
$$

其中 $P_{365, k}$ 为排列公式，等于 $365! / (365 - k)!$，则概率 $p$ 为

$$
p = 1 - \frac{365!}{(365 - k)!365^k}
$$

如果按照这个式子一个个计算里面的每一项时，需要很大的计算量，尤其是 $365!$，一般计算器估计算不出来。此时，我们将式子用指数，对数进行改写，可写为

$$
p = 1 - e^{\ln{(365)!} - \ln{(365-k)!} - k\ln{(365)}}
$$

其中 $\ln{(365)!}$ 和 $\ln{(365-k)!}$ 可以用Stirling公式进行近似计算，可以得到

$$
\begin{align}
& \ln{(365)!} - \ln{(365-k)!} - k \ln{(365)} \
\approx ; & \left[\frac{1}{2} \ln{(2 \pi)} + (365 + \frac{1}{2}) \ln{(365)} - 365 \right] \
& - \left[ \frac{1}{2} \ln{(2 \pi)} + (365 - k + \frac{1}{2}) \ln{(365-k)} - 365 + k \right] \
& - k \ln{(365)} \
\approx ; & (365 + \frac{1}{2}) \ln{(365)} - (365 - k + \frac{1}{2}) \ln{(365-k)} - k - k \ln {(365)}
\end{align}
$$

可见，阶乘和指数部分变成一堆对数的乘积以及求和，这样计算量就小了很多。这里贴一下用Python实现的概率计算程序

1
2
3
4
5
6
7
8
9
10
11
12

import math
import sys


if __name__ == '__main__':
    k = int(sys.argv[1])
    sum = ((365 + 0.5) * math.log(365) -
           (365 - k + 0.5) * math.log(365 - k) -
           k - k * math.log(365))
    prob = 1.0 - math.exp(sum)
    print('Probability (k = {}) is {}'.format(k, prob))

用程序计算了一些值，得到下面列表

从表格中可以看到，$k = 70$ 时，概率很接近 $1$ 了，原以为至少需要100多人才可能达到0.999这个概率值，可是现在看来是不需要达到100人就可以了。这让我想起了在班里，有人的生日和我是同一天，整个班的人数大概是50人，根据表格，可以得到概率为0.970，由此可见班里同学生日和我同一天的概率还是挺高的。