zhanxw

Nov 072013
 

安装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的心得。

Oct 182013
 

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的测试,还可以参考这个链接

Sep 262013
 

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.

Sep 082013
 

矩阵最大特征值的求法

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的。

 

 

Sep 052013
 

Bioinformatics Tools (build for Windows)

1. Tabix (Download)

2. Bgzip (Download)

These tools are build using MingW, with necessary changes to source codes:

  • Use _WIN32 definition to make sure fopen() use “rb” to open files, not “r”. (fopen is equivalent to _bgzf_open)
  • Link with zlib-2.1.8
  • Tested with example data.
Aug 282013
 

对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),把清理资源的任务留给操作系统。

 

Aug 212013
 

给作者编号

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来实现更自然、自动的布局。

 

Jul 082013
 

让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
可以看出使用四个线程和一个线程耗时差不多一样。但如果细看矩阵计算时间(这里忽略),四个线程还是能提速不少的。
但是考虑到使用方便,一个线程在大多数情况下应该就够用了。

Jun 292013
 

在SnowLeopard上安装MacFusion

Install MacFusion on Snow Leopard

Snow Leopard (OS X 10.6)已经是比较老的的系统。在这个系统上装MacFusion需要些额外的准备:

1. 装Fuse for OS X 链接

2. 装MacFusion 链接

之后,要按照这个方法(链接)把sshfs 拷贝到MacFusion 的路径里。

May 302013
 

如何做Research

Research Principles Revealed

 

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

1. 同事和合作者

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

2. 注重框架里的假设

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

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

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