【论文笔记】标准正交基和投影在分布式凸优化中的应用

写在前面

原论文标题:Distributed convex optimization via continuous-time coordination algorithms with discrete-time communication.

之所以想要写这篇文章,是因为最近在看的一篇分布式优化论文(Kia 2015)[1],其中有这么一个问题。

Πn=In1n1n1nT\Pi_n=I_n-\frac{1}{n}1_n 1_n^T。为了简便表示,可以定义rRnr\in\mathbb R^nRRn×(n1)R\in\mathbb R^{n\times (n-1)},且满足以下条件:

r=1n1nrTR=0RTR=In1RRT=Πnr=\frac{1}{\sqrt{n}} 1_n,\,r^TR= 0,\,R^TR=I_{n-1},\,RR^T=\Pi_n。

那么如何得到期望的RR

显然这里,[r,R][r, R]nn维空间的一组标准正交基(orthonormal basis)。现在已知rr需要求出所有RR的列向量。这就需要线性代数中标准正交基和投影的知识。

正交基与投影

一个nn维空间中任何一组线性无关(linearly independent)的向量,都是这个nn维空间的一组基。当这组基的向量两两垂直(即eiTej=0e_i^Te_j=0),则称为正交基(orthogonal basis)。而标准正交基只是将正交基又添加了一个条件,模长为1(即ei=1\|e_i\|=1)。一个空间可以有无数组基向量,正交基和标准正交基也同样是有无数组的。

一维投影相当于把一个向量投影到另一个向量,可以进一步求取正交基。以下图片来自Jakob_Hu的博客[2]

已知向量uuvv,按照如上图所示方法可以求取投影p=uTvu2up=\frac{u^Tv}{\|u\|^2}u。求出pp后,进一步求出正交向量vpv-p非常容易。这就解决了在二维空间中求取一组正交基的问题,正交基为ppvpv-p

高维投影和Gram-Schmidt过程

为了解决第一节提到的问题,我们需要一个能够通过任意维度的一组基构造空间的正交基的算法。

三维空间

以三维向量为例,假设存在一组三维向量,需要求出这三个向量所在空间的正交基,其中两个已经处理得到相互正交的向量p1p_1p2p_2。再加上向量ww构成三维空间的一组基。

下一步我们需要做wwp1p_1p2p_2构成的二维子空间的投影pp,则wpw-p即我们想要的正交向量。

  1. 分别做wwp1p_1p2p_2上的投影,得到wTp1p12p1\frac{w^Tp_1}{\|p_1\|^2}p_1wTp2p22p2\frac{w^Tp_2}{\|p_2\|^2}p_2,则p=wTp1p12p1+wTp2p22p2p=\frac{w^Tp_1}{\|p_1\|^2}p_1+\frac{w^Tp_2}{\|p_2\|^2}p_2
  2. 计算wpw-p,得到wwTp1p12p1wTp2p22p2w-\frac{w^Tp_1}{\|p_1\|^2}p_1-\frac{w^Tp_2}{\|p_2\|^2}p_2即第三维的正交向量。

三维空间求正交基的整个过程可以看做是先求出相应的二维空间的正交基,进一步求取三维空间正交基。

四维及以上空间

与三维空间相似,先求取低维度空间的正交基,在其基础上进行高一维度正交基的求取。

给出任何一组nn维空间的基,正交基的过程都可以通过逐一维度的计算得到。任何一个维度的向量都减去它在低维度空间中已经正交的向量的投影,这一过程就是Gram-Schmidt过程

Gram-Schmidt过程构造标准正交基

回到第一节提到的问题,现在我们已知nn维空间的任意一个向量rr,如何利用Gram-Schmidt过程构造一组标准正交基。

% 维度为4
n = 4;
r = ones(n,1)/n;

% r'的零空间
N = @(A) eye(size(A,2))-pinv(A)*A;
R1 = N(r');
R1 = R1./vecnorm(R1);

% [r R1(:,1)]'的零空间
R2 = N([r R1(:,1)]');
R2 = R2./vecnorm(R2);

% [r R1(:,1) R2(:,2)]'的零空间
R3 = N([r R1(:,1) R2(:,2)]');
R3 = R3./vecnorm(R3);
R = [R1(:,1) R2(:,2) R3(:,3)];

由此构建第一节需要的RR

后记

这里最后提一下,按照前面这样定义的RR会有什么好的性质。

首先,自不必说,有,

r=1n1nrTR=0RTR=In1RRT=Πnr=\frac{1}{\sqrt{n}}\boldsymbol 1_n,\,r^TR=\boldsymbol 0,\,R^TR=I_{n-1},\,RR^T=\Pi_n。

然后,P=[r,R]P=[r,R]显然是标准正交基组成的矩阵。由于PP满秩,所以

PTP=PPT=InRTP=[0n1,In1]ΠnP=RRTP=R[0n1,In1]=[0n1,R]P^TP=PP^T=I_n,\,R^TP=[0_{n-1},I_{n-1}],\, \Pi_n P=RR^TP=R[0_{n-1},I_{n-1}]=[0_{n-1},R]。

下面给出该性质在第一节所看论文中的应用。

Section 4.1中的应用

定义图G\mathcal G是强连通平衡图,其拉普拉斯矩阵为LL。现有动力学如下

v˙=αβLxx˙=αf(x)βLxv\begin{aligned} \dot v&=\alpha\beta Lx,\\ \dot x&=-\alpha\nabla f(x)-\beta Lx-v, \end{aligned}

其中1nTv=01_n^Tv=0,故也有rTv=0r^Tv=0。令vˉ=αf(xˉ)\bar v=-\alpha\nabla f(\bar x)

定义坐标变换

u=vvˉy=xxˉu=Pwy=Pz\begin{aligned} u&=v-\bar v,\,&y &=x-\bar x,\\ u&=Pw,\,&y&=Pz。 \end{aligned}

由于拉普拉斯矩阵的性质,PTLP=diag(0,RTLR)P^TLP=\operatorname{diag}( 0,R^TLR)。同时,rTPw=rT(vvˉ)=rTvˉr^TPw=r^T(v-\bar v)=r^T\bar vRTPw=w2:nR^TPw=w_{2:n}

定义h=f(y+xˉ)f(xˉ)h=\nabla f(y+\bar x)-\nabla f(\bar x),有

w˙1=0w˙2:n=αβRTLRz2:nz˙1=αrThz˙2:n=αRThβRTLRz2:nw2:n\begin{aligned} \dot w_1&=0,\\ \dot w_{2:n}&=\alpha\beta R^TLRz_{2:n},\\ \dot z_1&=-\alpha r^Th,\\ \dot z_{2:n}&=-\alpha R^Th-\beta R^TLRz_{2:n}-w_{2:n}。 \end{aligned}

Section 5.2中的应用

z~2:n(t)=z2:n(tk)z2:n(t)=RTP(z(tk)z(t))=RT(x(tk)x(t))\tilde z_{2:n}(t)=z_{2:n}(t_k)-z_{2:n}(t)=R^TP(z(t_k)-z(t))=R^T(x(t_k)-x(t)),则有

z~2:nT(t)z~2:n(t)=(x(tk)x(t))TRRT(x(tk)x(t))=(x(tk)x(t))TΠ(x(tk)x(t))\tilde z_{2:n}^T(t)\tilde z_{2:n}(t)=(x(t_k)-x(t))^TRR^T(x(t_k)-x(t))=(x(t_k)-x(t))^T\Pi(x(t_k)-x(t))。

由于Π\Pi是幂等的,即Π=Π2\Pi=\Pi^2。因此,z~2:n(t)=Π(x(tk)x(t))\|\tilde z_{2:n}(t)\|=\|\Pi(x(t_k)-x(t))\|

p=[zT,w2:nT]Tp=[z^T,w_{2:n}^T]^T,同样的由于Π\Pi的幂等性,有

pTp=zTz+w2:nTw2:n=yTPPTy+wPTRRTPw=(xxˉ)T(xxˉ)+(vvˉ)TΠ(vvˉ)p^Tp=z^Tz+w_{2:n}^Tw_{2:n}=y^TPP^Ty+wP^TRR^T Pw=(x-\bar x)^T(x-\bar x)+(v-\bar v)^T\Pi(v-\bar v)

因此,p=xxˉ2+Π(vvˉ)2\|p\|=\sqrt{\|x-\bar x\|^2+\|\Pi(v-\bar v)\|^2}

题外话

此外,简单记录下第一节所看论文的另外一些收获。

两个定理

证明收敛性和有界性用到的两个定理[3]

Lemma 3.4. (Comparison lemma) Consider the scalar differential equation

u˙=f(t,u),u(t0)=u0\dot u=f(t,u),\,u(t_0)=u_0

where f(t,u)f(t,u) is continuous in tt and locally Lipschitz in uu, for all t0t\geq 0 and all uJRu\in J\subset \mathbb R. Let [t0,T)[t_0,T) (TT could be infinity) be the maximal interval of existence of the solution u(t)u(t) and suppose u(t)Ju(t)\in J for all t[t0,T)t\in[t_0,T). Let v(t)v(t) be a continuous function whose upper right-hand derivative D+v(t)D^+v(t) satisfies the differential inequality

D+v(t)f(t,v(t)),v(t0)u0,D^+v(t)\leq f(t,v(t)),\,v(t_0)\leq u_0,

with v(t)Jv(t)\in J for all t[t0,T)t\in[t_0,T). Then, v(t)u(t)v(t)\leq u(t) for all t[t0,T)t\in[t_0,T).

可以知道,当v˙(t)u˙(t)\dot v(t)\leq \dot u(t)v(t0)u(t0)v(t_0)\leq u(t_0)时,v(t)u(t)v(t)\leq u(t)

Theorem 4.10. Let x=0x=0 be an equilibrium point for x˙=f(t,x)\dot x=f(t,x) and DRnD\subset \mathbb R^n be a domain containing x=0x=0. Let V:[0,)×DRV:[0,\infty)\times D\to \mathbb R be a continuously differentiable function such that

k1xαV(t,x)k2xα,Vt+Vxf(t,x)k3xα,\begin{aligned} k_1\|x\|^\alpha\leq V(t,x)&\leq k_2\|x\|^\alpha,\\ \frac{\partial V}{\partial t}+\frac{\partial V}{\partial x}f(t,x)&\leq -k_3\|x\|^\alpha, \end{aligned}

t0\forall t\geq 0 and xD\forall x\in D, where k1k_1, k2k_2, k3k_3, and aa are positive constants. Then, x=0x=0 is exponentially stable. If the assumptions hold globally, then x=0x=0 is globally exponentially stable.

可以知道,当上述定理中的条件满足时,我们有V˙k3k2V\dot V\leq -\frac{k_3}{k_2}V,即不光指数收敛,且收敛速度不小于k3k2\frac{k_3}{k_2}

由Lemma 3.4可知,

V(t,x(t))V(t0,x(t0))e(k3/k2)(tt0)V(t,x(t))\leq V(t_0,x(t_0))e^{-(k_3/k_2)(t-t_0)}。

因此,

x(t)[V(t,x(t))k1]1/α[V(t0,x(t0))e(k3/k2)(tt0)k1]1/α[k2x(t0)αe(k3/k2)(tt0)k1]1/α=(k2k1)1/αx(t0)e(k3/k2a)(tt0)\begin{aligned} \|x(t)\|&\leq \left[\frac{V(t,x(t))}{k_1}\right]^{1/\alpha}\leq\left[\frac{V(t_0,x(t_0))e^{-(k_3/k_2)(t-t_0)}}{k_1}\right]^{1/\alpha}\\ &\leq \left[\frac{k_2\|x(t_0)\|^\alpha e^{-(k_3/k_2)(t-t_0)}}{k_1}\right]^{1/\alpha}=\left(\frac{k_2}{k_1}\right)^{1/\alpha}\|x(t_0)\|e^{-(k_3/k_2a)(t-t_0)}。 \end{aligned}

微分方程求解

使用Lemma 3.4证明收敛性,本质是把一个复杂的微分方程求解问题,转化成一个相对容易的微分方程来求解。

第一节所看论文中证明有界性时的一个问题为例。已知

x˙a(1+x)+b(1+x)2,x(0)=0\dot x\leq a(1+x)+b(1+x)^2,\,x(0)=0

现在要证明xx在时间τ\tau内有界ζ\zeta,并求出τ\tau的值。

应用Lemma 3.4,若存在y˙=a(1+y)+b(1+y)2\dot y=a(1+y)+b(1+y)^2y(0)=0y(0)=0,那么xyx\leq yt0\forall t\geq 0

分离变量,求解微分方程,参考步骤[4]如下

t+C=dt=dya(1+y)+b(1+y)2=dyb(1+y+a2b)2a24b=1bdzz2a24b2,(z=1+y+a2b)=1a(dzza2bdzz+a2b)=1alnza2b1alnz+a2b=1aln1+y1aln1+y+ab\begin{aligned} t+C=\int dt &=\int \frac{dy}{a(1+y)+b(1+y)^2}\\ &=\int \frac{dy}{b\left(1+y+\frac{a}{2b}\right)^2-\frac{a^2}{4b}}\\ &=\frac{1}{b}\int \frac{dz}{z^2-\frac{a^2}{4b^2}} ,\quad \left(z=1+y+\frac{a}{2b}\right)\\ &=\frac{1}{a}\int \left(\frac{dz}{z-\frac{a}{2b}}-\frac{dz}{z+\frac{a}{2b}}\right)\\ &=\frac{1}{a}\ln \left|z-\frac{a}{2b}\right|-\frac{1}{a}\ln \left|z+\frac{a}{2b}\right|\\ &=\frac{1}{a}\ln \left|1+y\right|-\frac{1}{a}\ln \left|1+y+\frac{a}{b}\right| \end{aligned}

假设y,a,b0y,a,b\geq 0。令y=t=0y=t=0,可得aC=ln(ba+b)aC=\ln (\frac{b}{a+b}),则有

eat+aC=ba+beat=by+bby+a+be^{at+aC}=\frac{b}{a+b}e^{at}=\frac{by+b}{by+a+b}

变换为yy关于tt的函数

y(t,0)=(a+b)(eat1)beat+a+by(t,0)=\frac{(a+b)(e^{at}-1)}{-be^{at}+a+b}

再令y(τ,0)=ζy(\tau,0)=\zeta,得到

τ=1aln(1+aζa+b(1+ζ))\tau=\frac{1}{a}\ln\left(1+\frac{a\zeta}{a+b(1+\zeta)}\right)


  1. Kia, S. S., Cortés, J., & Martínez, S. (2015). Distributed convex optimization via continuous-time coordination algorithms with discrete-time communication. In Automatica (Vol. 55, pp. 254–264). Elsevier Ltd. https://doi.org/10.1016/j.automatica.2015.03.001 ↩︎

  2. Jakob_Hu. 线性代数(14)——正交性、标准正交基和投影. https://blog.csdn.net/Jakob_Hu/article/details/90813435 ↩︎

  3. Khalil, H. K. (2002). Nonlinear systems (3rd ed.). Prentice Hall, ISBN: 0130673897. ↩︎

  4. 积分(反导数)计算器. https://zs.symbolab.com/solver/integral-calculator ↩︎