Author: zhanxw

安装Theano

zhanxw Leave a comment

安装Theano
Install Theano

Theano是一个Machine Learning的工具箱，因为和Deep Learning相关，现在非常流行。

这篇文章简单记录一下安装Theano是遇到的问题。

1. import theano出错

>>> import theano
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/__init__.py", line 76, in <module>
    import theano.tests
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tests/__init__.py", line 6, in <module>
    import unittest_tools
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tests/unittest_tools.py", line 8, in <module>
    import theano.tensor as T
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tensor/__init__.py", line 10, in <module>
    import blas
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tensor/blas.py", line 393, in <module>
    StrParam(default_blas_ldflags()))
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/tensor/blas.py", line 332, in default_blas_ldflags
    blas_info = numpy.distutils.__config__.blas_opt_info
AttributeError: 'module' object has no attribute '__config__'

加入下面几行可以解决问题：

import numpy
import numpy.distutils
import numpy.distutils.__config__

2.Theano无法编译

运行示例程序是遇到编译错误：

Traceback (most recent call last):
  File "logistic_sgd.py", line 373, in <module>
    sgd_optimization_mnist()
  File "logistic_sgd.py", line 271, in sgd_optimization_mnist
    y: test_set_y[index * batch_size: (index + 1) * batch_size]})
...
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/gof/cmodule.py", line 1641, in compile_str
    return dlimport(lib_filename)
  File "/net/fantasia/home/zhanxw/python/lib/python2.7/site-packages/Theano-0.6.0rc3-py2.7.egg/theano/gof/cmodule.py", line 263, in dlimport
    rval = __import__(module_name, {}, {}, [module_name])
ImportError: ('/net/fantasia/home/zhanxw/.theano/compiledir_Linux-3.4.61-x86_64-with-debian-wheezy-sid-x86_64-2.7.5/tmphSYif1/baf8a6fdee34b135ea9110785e356489.so: undefined symbol: ATL_dptsyrk', '[_dot22(<TensorType(float64, matrix)>, W)]')

这里的提示信息是找不到ATL_dptsyrk函数。在没有管理员权限的服务器上，Theano找不到合适的ATLAS函数库。
解决方法是加入以下几行到 $HOME/.theanorc

[blas]
ldflags = -lf77blas -latlas -lgfortran

对于其他错误来讲，还有一些更一般的解决方法。

1. 可以打开Theano的日志信息，比如在程序中加入这两行：

import logging                                                                  |
#logging.getLogger("theano.gof.cmodule").setLevel(logging.DEBUG)

Theano在运行是可以输出更多的信息。

2. 保留临时文件
可以在 $HOME/.theanorc 中加入以下两行

[global]
nocleanup = False

这样可以保留Theano编译过程中的临时文件，便于调试。

题外话

在有管理员权限的Ubuntu系统， 可以用OpenBLAS来代替ATLAS，因为OpenBLAS是多线程的，而且速度更快。具体可以参考Theano 的手册。

现在Theano可以正常使用，希望接下来几篇blog有空能写写学习Deep Learning的心得。

R语言如何用索引

zhanxw Leave a comment

R语言如何用索引
How to use index in R

R语言有三种索引方式： [[, [, $
本文基于R语言的定义，用例子里说明这三种索引方式的使用条件和它们之间的细微差别。

首先，R里面有下面的索引方式：

x[i]
x[i, j]
x[[i]]
x[[i, j]]
x$a
x$"a"

总结并翻译一下：

1. 如果x是一个vector (包括matrix和array)，用[i]或者[[i]]得到的是第i个元素。而且一般用[而不是[[，用[[的话会失去变量的名字（name和dimnames）。比如：

aVec <- c(1,2,3)
names(aVec) <- c("x1", "x2", "x3")
> aVec[1]
x1 
 1 
> aVec[[1]]
[1] 1

2. 如果x是一个list，那么[i]得到的是一个list，用[[i]]得到的是第i个元素。aList[1]得到的是一个list，aList[[1]]得到的一个vector。

aList <- list(1, 2)
> aList[1]
$x1
[1] 1

> aList[[1]]
[1] 1

3. 如果索引是一个vector（比如一次选取多个元素），只能用[，不能用[[。选出来的元素类型，索引一个或者多个是一样效果。

4. 用$加一个变量来索引list，这个变量会按字面选取，不会先求值再选取。比如： x<-1; aList[[x]] 和 aList[[1]]不一样。5. 用$来索引list之外的元素，得到的是NULL关于各种index的测试，还可以参考这个链接。

Ridge regression

zhanxw Leave a comment

Ridge regression 相关的资料

以前学Ridge regression，大概记住了如何求回归系数 \(\beta^{ridge}\)，但有些细节不了解。我找到下面两个slide感觉又学到了新东西。

 

第一个讲义 【Link or RidgeRegression – Ryan Tibshirani】

http://www.stat.cmu.edu/~ryantibs/datamining/lectures/16-modr1.pdf

Slide 6：Shrink之后的回归系数看起来什么样？

Slide 7：为什么要center和normalize predictors？

Slide 13：当真正的回归系数很大的时候，Ridge regression还有好的表现么？或者说当true coefficient很大的时候，应该向什么方向调整\(\lambda\)可以减少prediction error？

 

第二个讲义 【Link or RidgeRegression-Rudyregularization】

http://www-stat.stanford.edu/~owen/courses/305/Rudyregularization.pdf

Slide 17：Ridge regression的回归系数的bias具体是多少？

Slide 23-24：Ridge regression和Principal component analysis （PCA）的关系是什么？

Slide 37：选择Ridge regression的tuning parameter时一般用cross validation。如果用GCV（Generalized cross validation）的话，可以有close form.

矩阵最大特征值的求法

zhanxw Leave a comment

矩阵最大特征值的求法

Largest Eigenvalue of Matrix

 

矩阵的特征值分解是在科学计算常见的线性代数方法。在一些情况下，我们不需要全部特征值，只需要一部分特征值。或者因为矩阵的一些性质，我们只对一部分特征值感兴趣。比如对于矩阵 B = A’ * A， 如果A是m by n的实矩阵，m很大，n很小。根据线性代数理论，A, B矩阵的Rank最多是n，并且B矩阵是半正定矩阵。对矩阵B做特征值分解时，其特征值都是非负实数，并且最多有n个正数。因此，如果求解B矩阵的特征值，只需要知道其最大的n个特征值，因为其余特征值都是零。

求解一般矩阵的特征值主要有两种：（1）特征多项式（2）Power Method（包括 Arnoldi方法和QR分解)。 对于实对称正定矩阵，可以使用Choleskey分解。这些方法可以求出所有的特征值。而对于我们特别关心的实对称(Hermitain/Symmetric)正定(Positive Definite)矩阵，比如协方差矩阵，这里我们关心的是其最大（或者最小）的特征值。这里常用的方法的方法叫做Linczos方法，可以看做是Arnoldi方法的特例。Linczos方法计算A的特征值是（1）先随机生成一个向量v，（2）再利用一系列向量： v, A*v, A*A*v, …， A^(n-1) * v 来得到一个特殊矩阵T，这个矩阵的大小是m by m，这里m往往远大于n。（3）得到T之后，可以以O(m^2)的复杂度来的到特征值。

实现Lanczos方法时会使用Multiple restart，最有名的实现是ARPACK。ARPACK是用Fortran实现的，对Fortran不熟悉就不容易直接调用ARPACK的函数。不过很多其他软件都对ARPACK的功能进行了封装。在R里面使用ARPACK，可以用igraph包，R-blogger有一个这方面的例子（链接）。其他软件中一般也都有用ARPACK来实现求一部分特征值的功能：Matlab的eigs，Scipy的eigs和eighs（分别对应实对称矩阵和Hermitian矩阵， 链接），Julia的eigs（链接）。

后话：现在开源的科学计算软件很多，当我们想寻找合适高效的数值计算工具时，可以从源代码出发。比如我们知道Matlab有eigs函数，就可以用“julia eigs”来看julia语言里是怎么实现eigs的。

 

 

对RAII的思考

zhanxw Leave a comment

对RAII的思考

Some random thoughts about RAII

写程序时常常需要申请系统资源，比如打开文件，申请一块内存。申请到这些资源后，在程序退出或者资源使用完毕后，应当正确的释放。如果不能正确释放，会造成一系列问题。比如申请的内存没有释放造成内存泄露Memory Leak，申请的进程锁没有被解锁Unlock，造成进程间的死锁DeadLock。在C++语言里，解决这类资源管理问题的管用手法是RAII (Resource Acuiquistion Is Initialization)

这篇笔记是对RAII的一点思考。

 

1. 什么是RAII

简单来讲，把获取资源的代码放到类的构造函数里，把释放资源的代码放到析构函数里。比如用ofstream file(“output.txt”) 可以打开文件，当file变量不起作用是，文件会被自动关闭。

比如下面这张图（from：The RAII Programming Idiom），看看这里面有多少地方需要写释放资源的代码。如果使用RAII，这些其实地方都不用留代码。

RAII Example
RAII Example

 

 

2. RAII的优缺点

RAII的好处是利用C++语言优势安全、正确的管理资源。同时RAII是C++建议的资源获取方式，这种代码可以被广大C++用户理解。

不方便之处是，使用RAII有一些陷阱。比如不要用RAII一次获取多个资源。

 

3.为什么C++有RAII

C++语言保证了一个类构造之后，析构函数会被自动调用。这个使用方式与资源管理的方式相似。因此可以用类的生命周期来管理资源。

 

4. 为什么C/Java/Python没有RAII

C语言没有原生的构造和析构函数，获取的资源不能有任何自动机制来释放。

Java/Python有语言中的支持，即Dispose Pattern。举例来说就是 try…catch..finally语句。使用者只要把释放资源的语句写到finally，资源就会被释放。

 

4. RAII 和Exception的关系

RAII和Exception紧密相关，更确切的说，构造函数和异常这两个特性在某种程度上互相依赖。

对于构造函数来说（获取资源的语句在构造函数里），构造函数没有返回值，因此想知道资源是否成功获取是不能从函数返回值来判断的，唯一可以用的手法是在资源获取失败时抛出异常。也就是说构造函数需要使用异常。

另一方面，使用异常之后，需要用构造函数来管理资源。因为异常抛出以后，很可能处理异常的代码和异常发生的代码不在一个层次（异常在Call Stack上逐层向上）。为了实现异常安全（Exception Safe），应该使用构造函数（另一个选择是智能指针，但智能指针有智能指针的问题，详见C++FAQ的讨论）。

对已有的C++代码来说，实现或检查代码是不是异常安全不是一个的简单人物。这种情况下，异常这个特性往往会被禁用（比如Google C++ style guide）。如果异常被禁用了，我们就没法从构造函数本身获知资源是否成功获取，那是不是说我们没法使用RAII特性呢？

答案是否定的。我们可以在获取资源后，用其他的类函数来检查资源获取是否成功。比如ostream::is_open()就可以检查文件是否被正常打开。

 

5. 怎么绕开RAII

在C等不提供RAII支持的语言里，可以直接绕开RAII，即保证获取资源后，程序的每一个出口都有释放资源的语句。

这种方法有可能造成多处重复的资源释放代码，或者使用goto语句把所有程序跳转到一处资源释放代码。

 

6. 实践中怎么用RAII

实践中除了把资源获取的语句写到构造函数，把资源释放的语句写到析构函数，还应当注意：

1）获取多个资源时，可以写在多个类的构造函数里，使得每一个类的构造函数对应一个资源。这样在任何资源获取失败时，已经获得的资源会得到释放

2）有时候获取资源失败等于程序失败（Fatal Condition），这种情况下可以直接退出（exit），把清理资源的任务留给操作系统。

 

给作者编号

zhanxw Leave a comment

给作者编号

Author Number

http://zhanxw.com/author

 

最近写了一个网页，主要解决写文章给作者编号的问题。这个问题怎么来的呢？现在的文章需要提供作者，作者工作单位和作者的贡献，但是提供的方式是给每个作者单位按照作者出现的顺序编号。比如我的文章有79个作者，现在突然要加一个新的合作者，把他拍到第50位，那么从第51位到最后一位作者的编号可能都需要改变。这个过程很容易出错。我想写个软件给所有人用，顺便凑个热闹，学习一下jQuery + Django，最后的作品放到我的网页(link)。在这个过程中，还学习了Nginx, Emacs web-mode, 这篇笔记主要是把我的经验记录下来。

 

1. 学习Django

Django是Python语言实现的Web架构，它最开始是用于展示新闻（Newroom），国内的豆瓣也用Django。作为Python的“粉丝”，我学一下Django，希望以后可以把有用的工具放到网上。

想入门，最好最省时间的方法是Django tutorial，就在Django的官方网页。这个Tutorial一共有六个部分，介绍的常用Django的功能，包括： 模型（Model），模板（Template），视图（View），静态文件（Static），管理界面（Admin），网址分发（URL Dispatching）等等。学习之后，我感到Django的强项是利用Python语言做到来简化数据库使用。对于一个简单的App，基本开发的流程是： （1）设定project 的settings.py 和urls.py，注册新的app ；（2）设定urls.py来确定网址和视图（View）的关系；（3）用HTML开发模板（/template/app/XXX.html）同时实现各种需要的视图（views.py）

架设Django，一般简易使用uWSGI。我使用的是Nginx处理静态页面，同时用uWSGI处理Django相关界面。简易先在本机用uWSGI调好程序，再放到服务器上并配置Nginx和uWSGI。

Djano默认使用Unicode，而我们一般都用str()，这是ASCII编码。两种编码对于字符串处理来讲（包括join, index）没有什么区别，但在print的时候，Unicode的字符应该先转成UTF-8字符，比如： print u”aaa”.encode(‘UTF-8’)

 

2. 学习jQuery

用jQuery的主要目的是用Ajax来更新网页数据，这样用户用起来有行云流水的感觉。jQuery有.get 和 .post两个方法，使用起来和访问网页很相似。不过这里面有几个陷阱。第一个是GET方式适合小数据，对于Django来讲，默认的大小是4096个字符，大一点的数据应该用POST；第二个陷阱就是POST，为了安全（CSRF: cross-site request forgery），Django要求POST的数据必须有csrftoken，一般的HTML表格Form必须有csrf的标记才会被Django接收。

为了调试Ajax的输入输出，简易用Firefox下的FireBug或者最新版本的Firefox，它们都可以显示Ajax请求的参数和返回值。最新版本的Firefox还可以给Javascript设置断点，这就更加降低了调试的难度。

jQuery的功能很多，这个网页把jQuery 1.9以及之前版本的功能用不同颜色区分出来，很方便查询（link）。

在这个网页中，我还用了handsontable，这是一个在jQuery的基础上开发的实用的javascript spreadsheet，语法简单，使用方便。

3. 一点Nginx经验

对Nginx来讲，我们的配置是用location语句来把特定的网址传给uWSGI进程。注意Nginx设置里如果有多个location语句，Nginx并不一定选择最先匹配的网址，而是选择最长的网址。比如：

location ~* /author/static/.+.(jpg|jpeg|gif|css|png|js|ico|xml)$ {
root /var/django/zhanxw/; # STATIC_ROOT
expires 30d;
}

location ~* ^.+.(jpg|jpeg|gif|css|png|js|ico|xml)$ {
access_log        off;
expires           30d;
root /var/www;
}

如果第一个location写成： “location /author/static”，那么.jpg之类的文件会被第二个location处理，这就不是Django static文件的正确处理方式。

 

4. 一点Emacs web-mode经验

用Emacs写Django的模板HTML，最好用的不是django-html-mode，而是web-mode。只有这个mode可以识别 “{% static ‘polls/index.html’ %}”这样的记号，并正确缩进。

此外web-mode可以自动补全HTML tab，比如你在<p>后面打</，web-mode会帮你补全</p>。还可以用C-c C-e b 和C-c C-e e跳到一对tag的最前面（beginning）和最后面（end）。

 

5. 一般性经验

最开始设计应该以最少功能，最小实现为好，不要一下子把界面设计复杂。可以想几个用例（Use Case），保证最基本最重要的功能，其他功能应该越少越好。

网页的布局应该少用<br/>这种硬回车。在有Bootstrap的情况下，完全可以用<div>和<p>来用更少的HTML tag来实现更自然、自动的布局。

 

让R更快

zhanxw Leave a comment

让R更快
Speed up R

用R 3.0.1来做基准，用Rscript运行R-benchmark-25.R计时。

1) 基准时间：

197.16u 1.74s 199.71r 1703184kB 0 Rscript R-benchmark-25.R

之后下载OpenBLAS，编译，然后把编译好的libopenblas.so 替换R自带的libRblas.so.
2) 测试单线程的速度：

Single thread OpenBLAS

105.83u 108.53s 58.04r 3072256kB 0 Rscript R-benchmark-25.R

 

3) 下面试试四个线程：

Four threaded OpenBLAS

export OPENBLAS_NUM_THREADS=4

79.39u 23.67s 58.54r 2122688kB 0 Rscript R-benchmark-25.R
可以看出使用四个线程和一个线程耗时差不多一样。但如果细看矩阵计算时间（这里忽略），四个线程还是能提速不少的。
但是考虑到使用方便，一个线程在大多数情况下应该就够用了。

如何做Research （Research Principles Revealed）

zhanxw Leave a comment

如何做Research

Research Principles Revealed

 

看了Jennifer Widom research_principles 的幻灯，觉得有些方面值得学习。

1. 同事和合作者

交流思想，互相启发，优势互补。比如Jennifer 擅长Details，而合作者Stefano注重Intuition。

2. 注重框架里的假设

学习知识的时候有不同的层次。一般来讲是从抽象到具体，比如一个领域的知识分划到不同的章节，每个章节讲具体的知识细节。

这就需要仔细品味每个章节的特点和章节之间的异同。

对于做研究来讲，更重要的是知识的假设。特别是这些假设不成立的时候，如何发展或变通现有的理论，这样或许能做出新的研究。