一些非光滑凸优化算法

在很多凸优化问题中会遇到目标函数非光滑的情况,例如稀疏优化中的\ell_1范数的正则项,将约束条件转化为可行集合的指示函数的无约束优化问题,目标函数为多个子目标的最大值的优化问题等.

我们将会考虑三种不同的非光滑优化算法:次梯度方法,近似点梯度方法和加速的近似点梯度方法.它们的误差和迭代次数k的关系大致上为\mathcal O(1/\sqrt{k}), \mathcal O(1/k)\mathcal O(1/k^2).

次梯度方法可以看成是通常针对光滑函数的梯度方法的一种扩展,它通过推广梯度的概念,给出了一个和梯度方法类似的迭代格式.对于一些目标函数能够分解为一个光滑项和非光滑项的情况,近似点梯度方法通过引入近似点算子来处理非光滑项,同时对于光滑项仍然使用梯度方法,它的迭代分别两步,首先用梯度方法更新光滑项,之后再将非光滑项近似点算子作用在更新后的结果上.如果光滑项是一个强凸的函数,那么这个方法能够达到线性的收敛率.加速的近似点梯度方法是由Nesterov提出的一个方法,它能够解决和近似点梯度相同的问题,但是能够有更快的收敛速度.

我们在第2节中会给出一部分后续需要的概念和命题.在第3节和第4节中我们会重点考察这三种优化算法的收敛性和收敛速度的证明,同时在第5节还会介绍针对一类目标函数可分但是约束条件不可分的优化问题的交替方向乘子法.最后在第6节和第7节中,我们在LASSO和离散最优传输问题中针对不同的方法进行测试.

(more…)

Read More

球坐标下Laplace算子的公式

前几天给人找了几道微积分题目,里面有一个是计算球坐标下面的Laplace算子的表达式。这当然可以暴力用复合函数的求导来算,但是可能比较复杂。事实上可以有一种积分的办法来算出这个结果。

我们知道,在标准坐标下,Laplace算子的形式为

 \Delta u = \frac{\partial^2u}{\partial x^2} + \frac{\partial^2u}{\partial y^2} + \frac{\partial^2u}{\partial z^2}

现在计算在球坐标

\begin{cases}
x = r\sin\theta\cos\varphi \\
y = r\sin\theta\sin\varphi \\
z = r\cos\theta
\end{cases}

下Laplace算子的表达式。首先计算其Jacobi矩阵

 \frac{\partial(x, y, z)}{\partial(r, \theta, \varphi)} = \begin{pmatrix}
\sin\theta\cos\varphi & \sin\theta\sin\varphi & \cos\theta \\
r\cos\theta\cos\varphi & r\cos\theta\sin\varphi & -r\sin\theta \\
-r\sin\theta\sin\varphi & r\sin\theta\cos\varphi & 0
\end{pmatrix}^T

对应逆矩阵为

 \frac{\partial(r, \theta, \varphi)}{\partial(x, y, z)} = \begin{pmatrix}
\sin \theta\cos\varphi & \sin\theta\sin\varphi & \cos\theta \\
\frac{1}{r}\cos\theta\cos\varphi & \frac{1}{r}\cos\theta\sin\varphi & -\frac{1}{r}\sin\theta \\
-\frac{\sin\varphi}{r\sin\theta} & \frac{\cos\varphi}{r\sin\theta} & 0
\end{pmatrix}

注意到这是一个正交矩阵。

我们利用Laplace算子的积分表达式,考虑

 \begin{aligned}\frac{ \partial{u(r, \theta, \varphi)}}{\partial(x, y, z)} \cdot \frac{\partial{v(r, \theta, \varphi)}}{\partial(x, y, z)} &= \frac{\partial{u(r, \theta, \varphi)}}{\partial(r, \theta, \varphi)}\frac{\partial(r, \theta, \varphi)}{\partial(x, y, z)}^T\frac{\partial(r, \theta, \varphi)}{\partial(x, y, z)}\frac{\partial{v(r, \theta, \varphi)}}{\partial(r, \theta, \varphi)} \\ &= u_rv_r + \frac{1}{r^2}u_\theta v_\theta + \frac{1}{r^2\sin^2\theta}u_\varphi v_\varphi \end{aligned}

对于任意的$v \in C_c^\infty$,利用分部积分

 \begin{aligned}
\int \Delta u v dxdydz &= -\int \nabla u \cdot \nabla v dxdydz \\
&= -\int \left( u_r v_r + \frac{1}{r^2} u_\theta v_\theta + \frac{1}{r^2\sin^2\theta}u_\varphi v_\varphi \right) r^2\sin\theta drd\theta d\varphi \\
&= \int \left( \left (u_r r^2\sin\theta \right )_r + ( u_\theta \sin\theta )_\theta + u_{\varphi\varphi} \right) v drd\theta d\varphi \\
\end{aligned}

再根据

 \int \Delta u v dxdydz = \int \Delta u v r^2\sin\theta drd\theta d\varphi

这样就可以得到

 \begin{aligned}
\Delta u &= \frac{1}{r^2\sin\theta}\left( \left (u_r r^2\sin\theta \right )_r + ( u_\theta \sin\theta )_\theta + u_{\varphi\varphi} \right) \\
&= \frac{1}{r} \frac{\partial^2}{\partial r^2}(ru) + \frac{1}{r^2\sin \theta} \frac{\partial}{\partial \theta}\left( \sin \theta \frac{\partial u}{\partial \theta} \right) + \frac{1}{r^2\sin^2 \theta} \frac{\partial^2 u}{\partial \varphi^2}
\end{aligned}

Read More

Bernstein多项式和连续函数的多项式近似

开学前随手翻了一下数值分析的教材,发现里面提到了一个定理

设$I = [0, 1]$,$f \in C(I)$,那么存在一串多项式 $\{ p_n \}$,满足 $\lim_{n \to \infty} \sup_{x \in I} |p_n(x) - f(x)| = 0$。

定理中的 $p_n$ 可以构造性的找出来,下面的Bernstein多项式就是一种可能

 p_n(x) = \sum_{i=0}^n {n \choose i} x^i(1-x)^{n-i} f(\frac{i}{n})

这有一个概率论版本的证明:

首先考虑一串独立的服从Bernoulli分布的随机变量$X_1, X_2, \cdots$,满足$\mathbb P[X_i = 1] = x$。令$S_n = \frac{1}{n} \sum_{i=1}^n X_n$,那么有

 \mathbb E[f(S_n)] = \sum_{i=0}^n {n \choose i} p^i(1-p)^{n-i} f(\frac{i}{n}) = p_n(p)

那么对于$\delta > 0$有如下估计

 \begin{aligned} |p_n(x) - f(x)| & = |\mathbb E[f(S_n) - f(x)]| \leq \mathbb E[|f(S_n) - f(x)|] \\ &\leq \mathbb E[|f(S_n) - f(x)| \chi_{\{|S_n - x| < \delta\}}] + 2\sup_{y\in I} |f(y)| \mathbb P[|S_n - x| \geq \delta] \end{aligned}

由于$f \in C(I)$,其在$I$上一致连续并且有界,对于$\varepsilon > 0$,有$\delta > 0$使得当$|a - b| < \delta$时有$|f(a) - f(b)| < \varepsilon$。另外,根据Chebyshev不等式并且注意到$\mathbb E[S_n] = x$

 \mathbb P[|S_n - x| \geq \delta] \leq \frac{\text{Var}[S_n]}{\delta^2} = \frac{x(1-x)}{n\delta^2} \leq \frac{1}{4n\delta^2}

最终我们有如下估计

 \sup_{x \in I} |p_n(x) - f(x)| \leq \varepsilon + \frac{\sup_{y\in I}|f(y)|}{2n\delta^2}

这样最后就可以得到相应的结论。同时,这个方法还可以构造出一系列类似的多项式。

Read More

Carmér‘s Theorem

我们知道,两个独立的正态随机变量加起来还是一个正态的随机变量。但是反过来呢?是不是两个独立的随机变量加起来是个正态,它们就都是正态呢?Carmér定理说的就是这样一件事情:假设X_1, \cdots, X_n是一组独立的随机变量,并且X_1+\cdots+X_n是服从正态分布的,那么X_1,\cdots,X_n都是服从正态分布的。

这是这学期概率论课上老师给的bonus问题,这个定理原始论文是“Über eine Eigenschaft der normalen Verteilungsfunktion”,一篇用德文写的。然后我看着它的公式大概猜了一下这篇文章在干啥…… 证明见下面正文,应该是对的吧。感觉这个证明把这学期实分析、复分析学的东西都用上了一些。

主要的思路是,先证明特征函数\mathbb E[e^{itX}]可以定义在整个复平面上,并且是一个解析函数。同时,证明这个函数的增长阶是2并且没有零点,之后用Hadamard因子分解定理说明它的形式。

还有一个叫做Raikov定理,和这个定理非常相似,把定理中的“正态分布”换成“Poisson分布”就是Raikov定理的内容。证明和这个定理也很像,不过需要先说明两个随机变量都是离散的。

(more…)

Read More

Pólya’s Urn

Pólya’s Urn是这样一个模型:假设现在有一个盒子,里面有r个红球和g个绿球,现在不断地随机从盒子取出一个球,之后再放回去,外加放入额外的c个和取出的球颜色相同的球。现在问,盒子中绿色球的比例在最后是否会趋向与某个分布?

实际上最后这个比例会是一个参数是g/cr/c的Beta分布。

首先我们设S_n = r + g + nc是第n次操作后盒子中所有球的数目。再设G_n, R_n分别是第n次操作后盒子中绿色球、红色球的数目。以及令X_n = G_n / S_n, \mathcal F_n = \sigma(G_k, S_k, 0 \leq k \leq n)。计算条件期望

 \begin{aligned}
\mathbb E[G_{n+1} | \mathcal F_n] &= \frac{G_n}{S_n}(G_n + c) + \frac{R_n}{S_n}G_n \\
&= (c + S_n) X_n = S_{n+1} X_n
\end{aligned}

可以看出\mathbb E[X_{n+1}|\mathcal F_n] = X_n\{X_n\}是一个鞅。再由于0\leq X_n \leq 1,我们知道X_n几乎处处收敛到一个随机变量X。现在考虑其分布。

首先可以计算出在第n步及以前,有m \leq n次取出绿色球的概率,并且写为紧凑的形式如下

 \begin{aligned}
&\mathbb P\left[ X_n = \frac{mc + g}{nc + g + r} \right] \\&= \frac{n!}{m!(n-m)!}\cdot \frac{\displaystyle\prod_{k=0}^{m-1}(k + g/c) \prod_{k=0}^{n - m - 1}(k + r/c)}{\displaystyle\prod_{k=0}^{n - 1}(\frac{g+r}{c} + k)}\\&= \frac{\Gamma(n + 1)}{\Gamma(m + 1)\Gamma(n - m + 1)}\cdot \frac{\Gamma(m + g/c)\Gamma(n - m + r/c)}{\Gamma(n + (g + r)/c)\text{B}(r/c, g/c)}
\end{aligned}

其中\Gamma\text{B}分别是Gamma函数和Beta函数。

接下来利用Stirling公式\Gamma(x) = \sqrt{2\pi/x} (x/e)^x (1 + O(1/x)),可以得到

 \frac{\Gamma(x + a)}{\Gamma(x + b)} = x^{a - b} \left( 1 + O\left( \frac{1}{x} \right) \right)\text{ as } x \to \infty

那么对于固定的\alpha, \beta, m满足0 < \alpha < \beta < 1以及\alpha n < m < \beta n,可以有

 \begin{aligned}
&\mathbb P\left[ X_n = \frac{mc + g}{nc + g + r} \right] \\&= \frac{1}{\text{B}(r/c, g/c)} n^{1 - (g + r)/c} (n - m)^{r/c - 1} m^{g/c - 1} \left( 1 + O\left( \frac{1}{n} \right) \right)\\
&= \frac{1}{n \text{B}(r/c, g/c)} \left (1 - \frac{m}{n} \right )^{r/c - 1} \left( \frac{m}{n} \right)^{g/c - 1} \left( 1 + O\left( \frac{1}{n} \right) \right)
\end{aligned}

对于所有满足\alpha n < m < \beta nm求和

 \lim_{n\to \infty} \mathbb P[\alpha \leq X_n \leq \beta] = \frac{1}{\text{B}(r/c, g/c)}\int_{\alpha}^\beta (1 - x)^{r/c - 1}x^{g/c - 1} dx

最后可以得到X的分布是\text{Beta}(g/c, r/c)

Read More

一道有意思的题目

前一段时间在做题的时候看到一道比较有意思的题目,大意是这样的

在单位圆周C上任意给出n个点w_1, \cdots, w_n,要求证明在这个圆周C上存在一个点w使得它到前面n个点距离的乘积为1。也就是

 \prod_{i=1}^n | w - w_i | = 1

这道题目是从 Stein 的 Complex Analysis 上来的,一眼看过去似乎和复分析没有任何关系… 但是却可以利用极大模原理来证明w的存在性。

大概是首先将w_1, \cdots, w_n放在复平面上,看成一个复数。然后定义函数

 f(z) = \prod_{i=1}^n (z - w_i)

这个函数可以定义在整个\mathbb C上并且是一个整函数。考虑其在z = 0处的模|f(0)| = \prod_{i=1}^n |w_i| = 1,那么按照极大模原理,f(z)在开单位圆\mathbb D内部不能达到极大模(除非它是常数),因此边界C = \partial \mathbb D上一定存在一个点z_0使得|f(z_0)| \geq 1。之后再利用|f(z)|的连续性,以及|f(w_1)|=1和介值定理就可以说明w的存在性。

证明还是挺有意思的,我想了一会儿似乎没想到直接通过几何的办法该怎么证…

Read More

多维正态分布中的边际分布、条件分布及Bayes公式

这个是概率论的一个期末小论文,大概删掉了开头结尾一些废话放在这里。

这篇文章主要是想要证明几个关于多维正态分布的很有用的定理,一是已知一个多维正态分布的联合分布,求边际分布以及条件分布的公式。二是已知一个条件分布和先验分布求后验分布的Bayes公式。这几个公式在进行统计的时候会比较有用。

预警:下面会有大量的公式,如果有公式恐惧症的话看每一节开头的结论就好了。

(more…)

Read More

三维空间中的旋转:旋转矩阵、欧拉角

考虑这样一个问题:如何计算三维空间中一个点绕着某一条向量旋转一个特定角度之后的坐标?旋转矩阵、欧拉角和四元数都是用来解决这个问题的方法。关于四元数的内容可以看平面等距变换、三维空间的旋转和四元数

接下来我们来讨论一下旋转矩阵和欧拉角这两个方法,并且我们选取右手坐标系作为我们的坐标系。

旋转矩阵

首先,对于一个三维空间的点 P(x, y, z),要将其绕 z 轴旋转 \theta 角度是可以很简单地用旋转矩阵来表示的

 \begin{equation}
R_z(\theta) =
\begin{bmatrix}
\cos \theta & -\sin \theta & 0 \\
\sin \theta & \cos \theta & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
\end{equation}

类似地,绕另外两个坐标轴旋转的矩阵可以表示如下(注意 R_y 有些不同)

 \begin{equation}
R_x(\theta) =
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos \theta & -\sin \theta \\
0 & \sin \theta & \cos \theta \\
\end{bmatrix}
\end{equation}

\begin{equation}
R_y(\theta) =
\begin{bmatrix}
\cos \theta & 0 & \sin \theta \\
0 & 1 & 0 \\
-\sin \theta & 0 & \cos \theta \\
\end{bmatrix}
\end{equation}

对于这三个特殊的旋转轴我们已经有了解决方案了,那么对于任意的轴 p 要怎么办呢?我们可以将这个旋转分解:

  • 将整个坐标轴旋转,使得旋转轴 pz 轴重合
  • 再将点 Pz 轴旋转 \theta
  • 再将整个坐标轴旋转回原位

3d-angle-rotation-matrix

如图,我们可以用两个角 \phi\psi 表示一个旋转轴的位置(这里认为旋转轴是个单位向量

 \begin{equation}
\left \{ \begin{aligned}
x &= \sin \phi \cos \psi \\
y &= \sin \phi \sin \psi \\
z &= \cos \phi
\end{aligned} \right.
\end{equation}

这样整个旋转就可以表示如下(最先进行的旋转它对应的旋转矩阵在最右侧)

R_z(\psi)R_y(\phi)R_z(\theta)R_y(-\phi)R_z(-\psi)

或者利用 R(-\alpha) = R^{-1}(\alpha) = R^T(\alpha),可以改写为

R_z(\psi)R_y(\phi)R_z(\theta)R_y^T(\phi)R_z^T(\psi)

经过化简,就可以得到最终的旋转矩阵

\begin{equation}
\begin{bmatrix}
\cos \theta + x^2(1 - \cos \theta) & xy(1 - \cos \theta) - z\sin \theta & xz(1 - \cos \theta) + y\sin \theta \\
yx(1 - \cos \theta) + z\sin \theta & \cos \theta + y^2(1 - \cos \theta) & yz(1 - \cos \theta) - x\sin \theta \\
zx(1 - \cos \theta) - y\sin \theta & zy(1 - \cos \theta) + x\sin \theta & \cos \theta + z^2(1 - \cos \theta)
\end{bmatrix}
\end{equation}

我们将旋转矩阵左乘需要旋转的向量就可以得到旋转后的结果了!

旋转矩阵一个很方便的地方是它可以沿着任意轴任意角度的旋转,但是,旋转矩阵缺点是它需要有 9 个元素来表示一个旋转,而且矩阵的乘法也比较慢。

欧拉角

euler-angle

欧拉角是用三个旋转角度 \alpha, \beta, \gamma 来标示旋转的。如图,图中蓝色坐标系是起始的坐标系,红色的坐标系是最后旋转完成的坐标系。整个旋转分为三个步骤:

  • 将坐标系绕 z 轴旋转 \alpha
  • 将旋转后坐标系绕自己本身的 x 轴(也就是图中的 N 轴)旋转 \beta
  • 将旋转后坐标系绕自己本身的 z 轴旋转 \gamma

由于绕不同的轴旋转所得到的欧拉角是不同的,所以欧拉角在使用的时候必须要先指明旋转的顺序,这里使用的是“zxz”的顺序。

欧拉角表示的旋转转换成旋转矩阵就是

 \begin{equation} R_z(\alpha)R_x(\beta)R_z(\gamma) \end{equation}

需要注意,这里的后面两次旋转并不是在原本固定的坐标系下的旋转。在旋转矩阵中,每一次旋转的叠加都是在左边乘上对应的旋转矩阵,然而,在这里乘法的顺序是相反的。

这可以这样来理解,假设最原始的固定的坐标系是 C_0

  • 假设有一个和 C_0 重合坐标系 C_3,先将 C_3C_0z 轴旋转 \gamma
  • 假设有一个和 C_0 重合坐标系 C_2,将 C_2 和前一步旋转后的 C_3 一起绕 C_0x 轴旋转 \beta
  • 假设有一个和 C_0 重合坐标系 C_1,将 C_1 和前一步旋转后的 C_2, C_3 一起绕 C_0z 轴旋转 \alpha

然后我们仅看这三个坐标系的关系:C_0 绕自己的 z 轴旋转 \alpha 角就可以和 C_1 重合;C_1 绕自己的 x 轴旋转 \beta 角可以和 C_2 重合,这是因为 C_1C_2 在最后一步是一起旋转的,它们的相对位置不会改变;同样可以知道 C_2 绕自己的 z 轴旋转 \gamma 角就可以和 C_3 重合。

在这里 C_1, C_2 相当于是前面欧拉角旋转的前两步得到的坐标系。

这样的话从 C_0C_3 的变换就相当于之前欧拉角的旋转变换,因此按照这个过程,旋转矩阵就是按照上面的顺序相乘了。

Read More

反演魔术:反演原理及二项式反演

首先我来说说什么是反演(inversion),对于一个数列\{f_n\}来说,如果我们知道另外一个数列\{g_n\},满足如下条件

 g_n = \sum_{i=0}^n a_{ni}f_i

反演过程就是利用 g_0, g_1, \cdots, g_n 来表示出 f_n

 f_n = \sum_{i = 0}^n b_{ni}g_i

本质上来说,反演其实是一个解线性方程组的过程,但是你观察后发现,这个矩阵实际上是一个三角矩阵,那必然存在着快捷的方法。对于使得这两个反演公式成立的这些系数应该满足什么条件呢?我们现在就来探讨一番!

我们首先讨论一下在什么情况下能够比较容易建立反演公式,之后介绍二项式反演以及它的两个应用,其中一个是错位排列问题

(more…)

Read More

特殊多项式在整点上的线性插值方法

首先我来解释一下题目是什么东西,假如给你一个 k 次多项式 P(x),并且已知了 P(0), P(1), \cdots, P(k),要求计算出 P(n), n \in \mathbb N 的值,由于已经知道了 k + 1 个点值,那么要插值出系数是可以在 \mathcal O(k^2) 完成,如果使用 FFT 的话是可以在 \mathcal O(k \log^3 k) 的时间内完成的[1],但是对于这个特别的问题,我们有一个 \mathcal O(k) 的算法!这好像是杜教先提出来的,然后还有这篇文章

(more…)

Read More