Category: 数值计算

从回转寿司到数值代数

最近一直在写毕业论文, 也想不起来更新博客. 发现很久没写了, 上来补一篇.
过年的时候和好友聚餐, 去吃回转寿司(配图从网上随便找的, 非去过的店铺, 仅供参考). 回转寿司里有一个很长的转盘, 厨师把寿司放到传送带上, 然后传送带不断旋转, 这样顾客能够随意挑选自己想要的寿司.

我的好朋友吃饱了以后就开始犯强迫症了. 他看转盘上的寿司摆放不均匀, 每次看到一个寿司到左边的寿司和到右边的寿司距离不相等, 他就把那个寿司挪一下位置, 让它刚好在左右两个寿司的正中间. 挪着挪着他就问我: 你说我这么一直挪着, 如果没其他人放寿司或者拿走寿司, 最终这些寿司能摆放均匀么?
我说我感觉上可以, 他也说感觉上可以. 我们怎么证明呢? 当时给了他一个解答, 后来又想到一个更有意思的解答. 仅仅用到一些数值代数知识.
记我们有$n$个寿司, 顺时针或者逆时针顺序编号, 第$i$个寿司到第$i+1$个寿司的距离是$d_{i}$. 因为是个转盘, 所以第$n+1$个寿司实际上就是第$1$个寿司, 为了描述方便, 后面凡出现$n+1$的时候都认为是$1$, 这样就不用特殊说明了. 另外$d_{1}+d_{2}+\cdots+d_{n}=L$, 也即转盘长度. 我好友的操作实际上就是对向量$d$左乘一个矩阵. 举个例子来说明一下, 为了写着方便, 以$n=3$为例, 他只会有挪动三个寿司的操作, 也就是三种不同操作, 分别相当于对$d$左乘的矩阵是这三种:
$
\left(
\begin{array}{ccc}
0.5 & 0.5 &0 \\
0.5& 0.5 &0 \\
0 & 0 & 1
\end{array}
\right)
$,
$
\left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & 0.5 & 0.5 \\
0 & 0.5& 0.5
\end{array}
\right)
$,
$
\left(
\begin{array}{ccc}
0.5 & 0 & 0.5 \\
0 & 1 &0 \\
0.5 & 0 & 0.5
\end{array}
\right)
$.
这种矩阵有什么特点? 稍微变型一下. 我们定义$n$个列向量$u_1, u_2, \ldots, u_n$: 其中$u_i$是这个样子的, 他的第$i$个分量是$1/\sqrt{2}$, 第$i+1$个分量是$-1/\sqrt{2}$, 其余分量全是$0$. 于是可以验证前面提到的矩阵实际上可以写成这样: $I-u_1u_1^T, I-u_2u_2^T, \ldots, I-u_nu_n^T$. 这个非常像Householder变换. 因为这几个矩阵都差不多, 我们拿出其中一个$I-u_iu_i^T$来分析一下.
很显然$I-u_iu_i^T$的秩是$n-1$, 它有唯一一个$0$特征值, $u_i$就是它的特征向量. 而其余特征值都是$1$.
我们需要想办法说明, 每一次挪动寿司都会导致一个量发生不可逆的变化, 这样才能证明, 不会操作了半天又退回到某个状态. 这里我们用向量的模长.
对于任意一个向量$d$, 如果将$I-u_iu_i^T$左乘$d$, 那么$d$的模长会变小, 除非$d$和$u_i$垂直. 具体是这样的:
$\| (I-u_iu_i^T)d \|^2 = \| d – u_iu_i^Td \|^2 = d^Td – 2(u_i^Td)^2 + (u_i^Td)^2 = d^Td – (u_i^Td)^2$
所以, 除非$d$和所有的$u_1,u_2,\ldots,u_n$都垂直, 否则在这$n$个矩阵的轮流作用下, $d$肯定是会一直变化的, 因为他模长变小了. 收敛到什么状态呢? 那就是$d$和所有的$u_1,u_2,\ldots,u_n$都垂直, 所以我们来看他们一列列堆起来构成的矩阵 $M=(u_1, u_2, \ldots, u_n)$, 如果收敛了, 只能是$M^Td=0$. 很显然$M$这个矩阵的秩至少是$n-1$, 因为我们观察他的前$n-1$列, 第$i$列的第$i$行元素非零, 其他列的那一行却是$0$, 所以线性无关. 因此$M^Td=0$的解空间至多是一维的. 容易找到一组解是$(1,1,\ldots,1)^T$, 所以所有解都是它的常数倍. 因此不停操作, 不停左乘矩阵, 不停迭代, 最终收敛到所有分量都相等, 也就是相邻寿司之间的距离都相等.
那从$d$出发开始迭代, 收敛到什么状态呢. 可以发现$d_1 + d_2 + \ldots + d_n=L$这个在操作过程中是保持不变的. 实际上看, 寿司距离加起来肯定永远是转盘长度. 从计算上看, 如果定义$e$是各个维度全是$1$的向量, 可以验证$e^T u_i = 0$ , 所以 $e^T (I-u_i u_i^T) d = e^T d – e^T u_i u_i^T d = e^Td = L $ , 这个就说明在每次变化后, $d_1 + d_2 + \ldots + d_n$不变. 所以最终收敛的时候, $d_i = L / n$.
也就是说, 经过我好友的不懈努力, 最终是可以让所有寿司均匀摆放的!

优化算法在应用问题中的常见技巧

这是今天去Hulu交流时用的PDF,比较简短。总结了在做应用问题时的一点点感觉。
点击下面链接下载。
Techniques for Optimization Methods in Applications

VPS、全局优化、Python、并行

说说最近的事,找不到一个合适的标题,就以若干关键字为题好了。

最近购买了一个VPS,新年特价,30刀一年,2G内存,续费仍然是30刀一年。对,你没看错,没有少写一个0
据说这家超售很严重,但是技术水平不错,不大看得出来。
有了VPS就可以继续肆无忌惮的挂网站、爬网页什么的了。实测CPU很给力,跑程序很快。但是毕竟是便宜货,刚入手两星期左右,今天差不多挂了10个小时才修复,看在价格上也就不说什么了。
科院选课助手 ishangke.net 目前已经迁移过去了。由于配置的提升,一些耗费内存、CPU的操作有明显变快。
预计一个月内可能会上线一个和音乐有关的小站,请期待。兴奋的是,以前从未发布过类似的东西。不过悲观的认为,不会有几个用户吧。
最近在想一个相关的问题,我如何通过网上能获得到的数据,进行音乐推荐呢?目前还没好主意。

这段时间在外访问,生活上比国内要无聊很多。这边的老师在做一些全局优化的问题,针对的是实际工程上的问题。我完全不懂背景,单从数学上来看,就是说有一个非常复杂的函数,复杂到求一次函数值需要解一大堆东西。函数本身也写不出表达式,只是给自变量能求出函数值而已。目的是用最短的时间尽可能找到函数最小值。由于求一次函数值花费时间很多,比如要将近一分钟。所以花费的时间就约等于函数值计算的次数。也就是要尽量少的找一些点求值,却希望找到最小值点。基本思想就是用一个函数去插值,但是为了求解到全局最优,又不能一直去找一堆距离很近的点,也要去尝试一些“未探索的区域”。所以对这两个目标进行一种权衡。

这边的程序都有Matlab和Python两个版本。老师说是Matlab要付费所以后来的程序都选择编写Python版本,于是果然看到了国外开始用Python进行科学计算的实例了。但是我用的太少,稍微补了一下Python和NumPy的基本知识,略有收获。第一步的任务只是用现有的程序去测试一个实际问题,感觉难度不大,具体意义大不大,由于不了解问题背景,我也没有发言权。我觉得有必要做点什么在国内宣传一下使用Python做科学计算,但是还没找到力所能及且自己觉得有意思的事情。而且我自己也才入门水平。

我感觉他这里的并行应该是希望能到集群上跑,但是他这里没有集群环境,用的是一些配置很不错的多核服务器。然后程序采用的是单机多进程,不是多线程。所以逻辑上来说移植到集群上并行跑问题不是太大,只要找到相应的工具应该就能办到。

随机样本选择——快速求解机器学习中的优化问题

前阵子去参加了数学规划会议,报告很多,人也很多。或者说报告和人过多了……
有少数感兴趣的报告,这里谈一下全场最后一个报告。报告人是Jorge Nocedal,就是著名的LBFGS的作者。
他关注的问题是一类机器学习中非常常见的优化模型:
$J(w) = \frac{1}{N} \sum_{i=1}^N l(f(w; x_i); y_i)$
他谈到了机器学习里面常用的随机梯度下降法(stochastic gradient descent, SGD)。在实际大数据的情况下,他比经典优化方法效果要好。因为经典优化方法算一次下降方向就已经需要很大计算量了。
SGD需要的一次迭代计算量非常小,但是代价是他非常可能不是下降方向,而且需要一个有经验的人选择很好的步长(学习速率)。并且,在理论上,SGD的性质少,所以在优化这个方向里没什么人提到。即使是在计算上,Nocedal也提到SGD有个很大的缺点是不容易并行。
以上这些现象我在学习推荐系统SVD模型的时候也有一些感触。
所以他采用了另外一个思路,思想简单,效果却很好。而且能很好的利用上经典优化算法。经典算法都需要计算梯度。他的思想就是:只选取一部分样本来计算梯度。
这样有一个近似的梯度计算方法,用
$J_S(w)=\frac{1}{|S|} \sum_{i \in S} l(f(w; x_i); y_i)$
代替原来目标函数,求它的梯度。
其中S是样本集,样本集的大小相对于N来说很小。但是他也不能太小。只要有了梯度,很多算法都能用上了,比如说他文章中提到梯度方法,Newton-CG方法。S中的元素可以在每次迭代的时候随机选取。
但是需要面临的一个问题就是如何确定S的大小。他给出了一个公式用来确定S的大小是否合适。
$\| \nabla J_S(w) – \nabla J(w) \|_2 \le \theta \| \nabla J(w) \|_2 $
这个式子的含义就是样本集上的近似梯度和真实梯度足够接近。当然这个又回到了老问题,就是真实梯度计算量太大。于是他又做了一些近似变换,使得只需要计算在样本集S上的方差即可近似估计上式是否满足。
$ \frac{ \| Var^2_S( l(w; i) ) \|_1 }{ \| S \| } \le \theta ^2 \| \nabla J(w) \|_2 ^2 $
这样就可以每次迭代验证一下S是否足够大,如果不够大就扩大一下。
这个技巧部分解决了我对大规模优化计算的一些困惑,他可以利用上很多经典优化算法,并且又能让每次迭代的运算量不会太大。这个技巧也使得并行计算容易用上去了。
后面考虑去使用一下,算算某个实际例子测试一下效果。
Nocedal这篇论文的链接:http://users.eecs.northwestern.edu/~nocedal/PDFfiles/dss.pdf

TopCoder Asia-Pacific 2012 解题思路

上月偶然看到网上有提到TopCoder Asia-Pacific 2012 (TopCoder亚太2012),比赛形式为马拉松,题目和传感器定位问题有些类似。因为我今年比较忙,没报名KDD CUP,也没报名TCO,GCJ,这次感觉应该有时间把这个题目做完,就尝试了一下。
题目在此 http://community.topcoder.com/longcontest/?module=ViewProblemStatement&rd=15137&pm=11942
题目大致是说有$n (100 \le n \le 2000)$ 个10维空间上的点,每一维坐标从[0, 1000]的整数中随机选取。对于每一个点,你知道到它最近和最远的点的距离,也知道其他点到这个点距离的远近排序。每个点的每个维度都有一定概率已知,需要算出不知道的维度。

Continue reading »

大规模问题快速算法的困惑

最近所里请了南京大学何炳生老师来做系列报告,他主要研究一阶算法,所谓一阶算法,只利用一阶最优性条件,通常借助一阶信息(梯度,Jacobi矩阵),不用二阶信息(Hession)。对于大规模问题来讲,这样做会节约很多计算量。当然他主要做理论,演示了如何利用他的方法简单和巧妙的证明一些算法的全局收敛性,感觉很有收获。

Continue reading »

为什么很多网络流问题总有整数解——全单模矩阵的性质

以前在Online Judge上遇到过很多问题需要转化为网络流模型,原问题一般是要求有整数解的。但是转化为网络流问题后虽然边上的最大流上限都是整数,其实并没有要求最大流也是整数,可求出来的最大流就是整数。我以前一直迷茫这个现象,并且我是很懒的人,一直没有严谨的想过这个问题或者查过这个问题。我们课题组平时很少讨论整数规划的问题,所以这个其实比较基本的问题我一直都不知道。其实你可以从直观感觉如果有非整数的流,那总应该能再怎么样多流一点变成整数,不过这显然是直观感觉,什么也说明不了。昨天看书看到有讲这些,终于理解了,在这总结一下。

Continue reading »

The Art of Lemon队的KDD CUP 2011 Track 2解决方案大致思路

随着KDD CUP 2011的结束,需要开始总结我们的解决方案了。我们在最终测试集Test2中排名第二,和在排行榜中测试集Test1上的排名是一致的。我先发一篇Blog大致总结一下我们的方案,一来自己回顾和理清整个过程便于后面详细的写Solution Paper,二来与大家分享我们队的成果。

Continue reading »

科学计算的新武器——Python

Python其实已经诞生了很多年,毫无疑问他是一个高级程序语言,写起来非常简洁,一行Python语句大约相当于六行C语句。Python起初并没有流行起来,但这些年非常火,究其原因可能是因为Python的发明者留大胡子了,呵呵。

Continue reading »

我理解的SVD

本科学数学的时候就发现,有一些定义、定理确实是需要仔细琢磨的。课本从来不会或者很少会告诉我们为什么需要某个定义,为什么如此定义,为什么要提出这样的定理,这个定理有什么意义。以前偶尔还会自己去思考一下其中的道理,但是慢慢的也变得懒了。最近看到有人转载matrix67的博文,觉得确实有必要重新的审视一下自己学过的东西。

Continue reading »