共轭梯度法初步

共轭梯度法及其浅显分析

更多内容请关注我的个人博客
浩源小站

背景

无处不在的线性方程组,昭示着一个优秀的算法能带来的巨大效益。对于简单的低阶方程组,只需要用普通的高斯消元法就能得到相当不错的结果,但是对于大型的方程组,常常与之相伴的就是稀疏性,如果还坚持使用高斯消元法,那么将带来的是空间上的巨大浪费。为此,迭代法求解方程组便应运而生。

共轭梯度法就是其中的佼佼者。

由来

二次型的最小值

我们熟知的线性方程组,常常可以写成Ax=b的形式,实际上当A是实对称矩阵时,就是二次型f(x)=\frac{1}{2} x^TAx-b^Tx+cx的导数为0时的表达式。

f(x)=\frac{1}{2}x^TAx-b^Tx+c\\ f'(x)=\frac{1}{2}A^Tx+\frac{1}{2}Ax-b

A​为实对称矩阵的时候,f'(x)=0​就等价于Ax=b​.

而我们求解线性方程组的问题,也可以转化成求解minf(x)

举个例子,当

A=\left[ \begin{matrix} 3&2\\ 2&6 \end{matrix} \right], b=\left[ \begin{matrix} 2\\8 \end{matrix} \right],c=0


时,f(x)​的图如下

image

等高线图如下

image

看得出,此时的函数具有唯一的最小值。

由代数学的知识,我们知道,当矩阵A​为正定、半正定、负定、以及不定的时候,方程组Ax=b​具有不同的解的情况。对应f‘(x)=0​则是不同的最小值情况。

下面这幅图说明了矩阵A的性质对于f(x)的影响,依次是正定、负定、半正定和不定的情况。

A的性质对于f的影响

可以看到,当A不定的时候,f(x)具有鞍点,这个时候无法通过导数为0来求得最值。

在接下来的共轭梯度方法中,我们都首先假定,A具有良好的性质,也就是对称和正定。

最速下降法(梯度下降法)

如何得到f(x)的最小值?

有一个形象的方法,从任意一点x_{(0)}出发,每次沿着当前位置下滑最快的方向走,也就是该点处的梯度方向。

从而得到一系列的解序列:x_{(1)},x_{(2)},\dots直到两次下降的差小于给定的误差限为止。

由上面的结论,f'(x_{(i)})=Ax_{(i)}-b,下面记误差(error)为e_{(i)}=x_{(i)}-x,残差(residual)为r_{(i)}=b-Ax_{(i)}

于是有

\begin{array}{**lr**} r_{(i)}=-f'(x_{(i)})\\ x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)} \end{array}

其中\alpha_{(i)}为第i步沿着方向r_{(i)}的步长

如何选取步长\alpha_{(i)}呢?朴素的想法是选取的步长尽量要让f(x_{i+1})的值最小,这样才能更快的到达f(x)的全局最小值。

从简单的微积分知识中,我们把f(x_{(i+1)})看作是\alpha_{(i)}的函数,要使得步长最合适,也就相当于

\frac{d}{d\alpha_{(i)}}f(x_{(i+1)})=0

根据链式法则

\begin{align*} \frac{d}{d\alpha_{(i)}}f(x_{(i+1)})&=f'(x_{(i+1)})^T\frac{d}{d\alpha_{(i)}}x_{(i+1)}\\ &=-r_{(i+1)}^Tr_{(i)} \end{align*}

也就是说,r_{(i)}r_{(i+1)}​正交

\begin{align*} r_{(i+1)}^Tr_{(i)}&=0\\ (b-Ax_{(i+1)})^Tr_{(i)}&=0\\ (b-A(x_{(i)}+\alpha_{(i)}r_{(i)}))^Tr_{(i)}&=0\\ (b-Ax_{(i)})^Tr_{(i)}-\alpha_{(i)}(Ar_{(i)})^Tr_{(i)}&=0\\ (b-Ax_{(i)})^Tr_{(i)}&=\alpha_{(i)}(Ar_{(i)})^Tr_{(i)}\\ \alpha_{(i)}&=\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}}. \end{align*}

也就是说,每一步的步长都可以根据当前的残差r_{(i)}来确定

总结一下,最速下降法就是:

\begin{align*} r_{(i)}&=b-Ax_{(i)}\\ \alpha_{(i)}&=\frac{r_{(i)}^Tr_{(i)}}{r_{(i)}^TAr_{(i)}}\\ x_{(i+1)}&=x_{(i)}+\alpha_{(i)}r_{(i)}. \end{align*}

这样迭代下去,直到达到事先给定的误差限\epsilon为止

当然,有的时候,出于减小复杂度的考虑,我们会换用这种手段来求r_{(i)}​

r_{(i+1)}=r_{(i)}-\alpha_{(i)}Ar_{(i)}

这个公式是在x_{(i+1)}=x_{(i)}+\alpha_{(i)}r_{(i)}的两边左乘-A再加上b得到的。

这样每次迭代过程我们只需要算一次矩阵乘法,可以减少运算。

雅可比迭代法

在开始共轭梯度法前,我们还需要再来看看它的同行们都是怎么干活的。

对于方程Ax=b,我们把A分成两个矩阵DE,其中D是对角阵,而E是剩下的部分,也就是说:

\begin{align*} Ax&=b\\ (D+E)x&=b\\ Dx&=b-Ex\\ x&=D^{-1}(b-Ex). \end{align*}

这里由于D是对角阵,所以能保证它的可逆性。

这里为了便于分析,把上面的式子写作

x=Bx+z

其中B=-D^{-1}E,z=D^{-1}b

下面考察该方法迭代过程的收敛性,依旧把误差(error)记作e_{(i)}

\begin{align*} x_{(i+1)}&=Bx_{(i)}+z\\ &=B(x+e_{(i)})+z\\ &=Bx+z+Be_{(i)}\\ &=x+Be_{(i)}\\ \end{align*}

故而

e_{(i+1)}=Be_{(i)}

这里判断迭代过程的敛散性主要看矩阵B的性质,如果B是实对称矩阵,那么e_{(i)}就可以用B的特征向量线性表示,迭代过程就可以写成:(下面以二阶矩阵为例)

\begin{align*} e_{(i+1)}&=Be_{(i)}\\ &=B^ie_{(0)}\\ &=B^i(x_1+x_2)\\ &=\lambda_1^ix_1+\lambda_2^ix_2 \end{align*}

其中\lambda_1,\lambda_2分别是矩阵B的特征向量x_1,x_2的特征值。

由此看得出,如果B的最大的特征值小于1,这个迭代过程就收敛,不然迭代过程会发散。

而一个矩阵的最大特征值也就是它的谱半径(spectral radius)\rho(B)

对于B不是实对称矩阵的情况,也有类似的结论,只是说明过于复杂,在此不提。

总之,雅可比迭代法收敛的充分条件是迭代矩阵B的谱半径\rho(B)<1​

下面还是以最初的方程举例子,来分析雅可比迭代的具体过程,这有助于改进迭代过程,得到更好的算法

迭代过程现在是

x_{(i+1)}=\left[ \begin{matrix} 0&-\frac{2}{3}\\ -\frac{1}{3}&0 \end{matrix} \right]x_{(i)}+ \left[ \begin{matrix} \frac{2}{3}\\ -\frac{4}{3} \end{matrix} \right]

此时B的特征向量以及与其对应的特征值分别为

$$
x_1=\left[
\begin{matrix}
\sqrt2\1
\end{matrix}
\right],\lambda_1=-\frac{\sqrt2}{3}\
x_2=\left[
\begin{matrix}
-\sqrt2\1
\end{matrix}
\right],\lambda_1=\frac{\sqrt2}{3}\

$$

下图是迭代过程中,二次型函数f(x)的等高线图,注意两个箭头分别表示特征向量

image

而下图则是迭代过程

image

考虑每一次迭代与特征向量的关系

image

我们会发现,每一次的迭代都朝着两个特征向量的线性组合方向前进。

最速下降法的分析

前面我们得到了最速下降法的迭代过程,并且得到了一个十分有意思的结论

$$

r_{(i+1)}^Tr_{(i)}=0

$$

这意味着,最速下降法每一次的迭代过程,下降的方向都与上一次的方向正交。

用具体例子来表述是这样的:

image

可以看到,(b)是最差的情况,需要蜿蜒前进很多次才接近解,而其他比较好的情况则两次迭代就收敛了.

但即便迭代的次数不同,每一次迭代的方向都与上一次正交,这就是最速下降法的特点。

在(b)中还有一个情况值得注意,由于这是一个二维的图形,故而如果每一次都和上一次正交,那么相邻一次迭代的方向是平行的!

换句话来说,在最速下降法中,有一个很糟糕的现象,为了收敛到解附近,同样的迭代方向可能走了不止一次。

这个现象不只是最速下降法中会出现,雅可比迭代法中也同样的出现了,表现就是每一次的方向都是特征向量的线性组合,而且大多数情况下,前一次迭代过的特征向量方向上的分量,在下一次的迭代中还继续存在。

共轭的想法

如果我能选取一系列线性无关的方向向量\{d_{(0)},d_{(1)},d_{(2)},\dots,d_{(n-1)}\},沿着每个方向只走一次,最后就能到达解x处,是不是能很好的解决前面雅可比迭代和最速下降法的问题?

最直接的想法来源于直角坐标系,如果每一个方向都是正交的,既自然,有具有很好的性质,是不是能得到很好的结果呢?

如果这样选取了方向,那么应该有误差e_{(i+1)}与方向d_{(i)}正交

\begin{align*} d_{(i)}^Te_{(i+1)}&=0\\ d_{(i)}^T(e_{(i)}+\alpha_{(i)}d_{(i)})&=0\\ \end{align*}

所以

\alpha_{(i)}=-\frac{d_{(i)}^Te_{(i)}}{d_{(i)}^Td_{(i)}}

但是这里步长大小依赖于当前的误差e_{(i)},而不知道正确解x是无法求出误差来的,故而这种想法虽然很自然,但是无法得到有用的结果。

考虑到矩阵具有变换向量的作用,不妨记x^TAy=0为向量x,y关于矩阵A共轭。

接下来选取的一系列方向向量都关于矩阵A两两共轭,\{ d_{(0)},d_{(1)},\dots,d_{(n-1)}\}

我们从上一个想法得到灵感,要求这里的e_{(i+1)}也与d_{(i)}共轭,这是很自然的,只需要想象是对原本的空间做了一个变换,原本不正交的向量在这里正交,那么这里的要求就是合乎情理的了。

下面也一样,对步长大小进行简单的分析

\begin{align*} \frac{d}{d\alpha_{(i)}}f(x_{(i+1)})&=f'(x_{(i+1)})^T\frac{d}{d\alpha_{(i)}}x_{(i+1)}\\ f'(x_{(i+1)})^Td_{(i)}&=0\\ d_{(i)}^Tr_{(i+1)}&=0\\ d_{(i+1)}^TA(e_{(i)}+\alpha_{(i)}d_{(i)})&=0\\ \end{align*}

所以得

\alpha_{(i)}=\frac{d_{i}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}}

这里得到的步长不仅可以算,而且可以证明,这样迭代,对于n阶的方程组,至多n步就可以收敛到正确解。

简单的证明

首先,把误差e_{(0)}表示为方向向量的线性组合e_{(0)}=\sum_{j=0}^{n-1}\delta_jd_{(j)}

两侧左乘d_{(k)}^TA 有:

\begin{align*} d_{(k)}^TAe_{0}&=\sum_{j=0}^{n-1}\delta_{j}d_{(k)}^TAd_{(j)}\\ &=\delta_{k}d_{(k)}^TAd_{(k)} \end{align*}

上式利用了之前提到的共轭性质,由此得

\delta_{k}=\frac{d_{(k)}^TAe_{(k)}}{d_{(k)}^TAd_{(k)}}

与上面的步长比较,有

\alpha_{(k)}=-\delta_k

于是对于第i步迭代的误差e_{(i)}

\begin{align*} e_{(i)}&=e_{(0)}+\sum_{j=0}^{i-1}\alpha_{(i)}d_{(j)}\\ &=\sum_{j=0}^{n-1}\delta_{(j)}d_{(j)}-\sum_{j=0}^{i-1}\delta_{(j)}d_{(j)}\\ &=\sum_{j=i}^{n-1}\delta_{(j)}d_{(j)} \end{align*}

由此可知,当i=n的时候,e_{(n)}=0

Q.E.D

并且我们能从上面的证明中得到这样一个事实,第i迭代时,有关前i-1次迭代的方向上都已经没有误差了,也就是这样的方法,每次都完全消除某个方向上的误差,这样至多n步得到精确解就是很自然的了。

共轭梯度法

上面的由来中,已经充分做好了准备来得到共轭梯度法,现在只剩一个问题,那一组两两共轭的方向向量该如何选取?

选取了一组基后,只需要进行Gram-Schmidt正交化就可以得到一组正交基,那么同样的也可以用这样的方法来得到一组共轭的方向向量。

问题是,如果随便选取一组线性无关的方向基,就进行Gram-Schmidt共轭化是否足够好呢?

我们总希望这组方向基不仅容易得到,而且具有很好的性质,减少计算量。

考虑到之前得到的结论,每一次迭代之间的残差都是相互正交的,我们不妨就选残差

\{ r_{(0)},r_{(1)},r_{(2)},\dots,r_{(n-1)} \}

作为共轭化之前的基。由于使用共轭化的方向向量来迭代至多只有n步,且每步都将该方向上的误差消灭掉,故而这一组基不仅线性无关,而且由于它们是残差,还具有正交的良好性质。

推导

首先,对于第i​+1次迭代后的残差r_{(i+1)}​

\begin{align*} r_{(i+1)}&=-Ae_{(i+1)}\\ &=-A(e_{(i)}+\alpha_{(i)}d_{(i)})\\ &=r_{(i)}-\alpha_{(i)}Ad_{(i)}. \end{align*}

接着我们来推导一下Gram-Schmidt共轭化的公式

i次的方向向量为d_{(i)},我们希望d_{(i)}r_{(i)}+\{ d_{(0)},d_{(1)},d_{(2)},\dots,d_{(i-1)}\}中得到

,也就是

d_{(i)}=r_{(i)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)},(i>k)

在上式的两侧同左乘上d_{(i)}^TA

\begin{align*} d_{(i)}^TAd_{(j)}&=r_{(i)}^TAd_{(j)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)}^TAd_{(j)}\\ 0&=r_{(i)}^TAd_{(j)}+\beta_{ij}d_{(j)}^TAd_{(j)},\quad i>j\\ \beta_{ij}&=-\frac{r_{(i)}^TAd_{(j)}}{d_{(j)}^TAd_{(j)}} \end{align*}

下面来结合残差r_{(j+1)}=r_{(j)}-\alpha_{(j)}d_{(j)},继续做一些推演

在两侧同左乘r_{(i)}^T

$$
\begin{align*}
r_{(i)}Tr_{(j+1)}&=r_{(i)}Tr_{(j)}-\alpha_{(j)}r_{(i)}^TAd_{(j)}\
\alpha_{(j)}r_{(i)}TAd_{(j)}&=r_{(i)}Tr_{(j)}-r_{(i)}^Tr_{(i)}\

\end{align*}
$$

故而

r_{(i)}^TAd_{(j)}= \left\{ \begin{align*} \frac{1}{\alpha_{(i)}}r_{(i)}^Tr_{(i)},\quad i&=j\\ -\frac{1}{\alpha_{(i-1)}}r_{(i)}^Tr_{(i)},\quad i&=j+1\\ 0,\quad otherwise \end{align*} \right.

结合\beta_{ij}的值,有

$$
\beta_{ij}=
\left{
\begin{align}
\frac{1}{\alpha_{i-1}}\frac{r_{(i)}Tr_{(i)}}{d_{(i-1)}TAd_{(i-1)}}&,\quad i=j+1\
0&,\quad i>j
\end{align
}

\right.
$$

考虑到完全可以把\beta_{ij}写成\beta_i,因为j只能取固定的值。

神奇的事情发生了,原本需要i\beta​才能确定的方向向量,在共轭化的条件下,只需要当前的数据和前一步的数据就可以得到,而不必存储之前所有走过的路径信息。

结合前面得到的\alpha_(i)=\frac{d_{(i)}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}}

\beta_{i}=\frac{r_{(i)}^Tr_{(i)}}{d_{(i-1)}^Tr_{(i-1)}}

实际上,上面的式子还能够写成更加漂亮的样子

先由d_{(i)}=r_{(i)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)},(i>k)

做其与r_{(j)}的内积有

\begin{align*} d_{(i)}^Tr_{(j)}&=r_{(i)}^Tr_{(j)}+\sum_{k=0}^{i-1}\beta_{ik}d_{(k)}^Tr_{(j)}\\ 0&=r_{(i)}^Tr_{(j)},\quad i<j \end{align*}

j=i

d_{(i)}^Tr_{(i)}=r_{(i)}^Tr_{(i)}

带入\beta_i中有

\beta_i=\frac{r_{(i)}^Tr_{(i)}}{r_{(i-1)}^Tr_{(i-1)}}

这样,Gram-Schmidt共轭化就完美的实现了,不仅实现了每一个方向只迭代一次,而且需要存储的数据只有上一步的残差。

下面是共轭梯度法涉及到的所有公式

\begin{align*} d_{(0)}&=r_{(0)}=b-Ax_{(0)}\\ \alpha_{(i)}&=\frac{r_{i}^Tr_{(i)}}{d_{(i)}^TAd_{(i)}}\\ x_{(i+1)}&=x_{(i)}+\alpha_{(i)}d_{(i)}\\ r_{(i+1)}&=r_{(i)}-\alpha_{(i)}Ad_{(i)}\\ \beta_{(i+1)}&=\frac{r_{(i+1)}^Tr_{(i+1)}}{r_{(i)}^Tr_{(i)}}\\ d_{(i+1)}&=r_{(i+1)}+\beta_{(i+1)}d_{(i)} \end{align*}

最后给出共轭梯度法的Matlab实现

function x=Conjugate_gradient(A,b,x0)
x=x0;
d=b-A*x0;
r=d;
for k=0:size(A)-1
    if r==0
        break;
    end
    alpha=(r'*r)/(d'*A*d);
    x=x+alpha*d;
    temp=r'*r;
    r=r-alpha*A*d;
    belta=(r'*r)/(temp);
    d=r+belta*d;
end
end

分析

正交性是共轭梯度法成功的关键,注意到每一次迭代中,新的残差r_{(i+1)}与前面所有的r_{(i)}正交,如果一个r_{(i)}=0那么Ax_{(i)}=b方程已经解完了。

否则在n步迭代后,r_{(n)}和n个两两正交的向量\{ r_{(0)},r_{(1)},\dots,r_{(n-1)}\}所张成的空间正交,而r_{(0)},r_{(1)},\dots,r_{(n-1)}这n个向量是R^n空间中的所有正交向量,由此r_{(n)}必须是零向量,所以Ax_{(n)}=b

比较

共轭梯度法从某种程度上要简单于高斯消元法,不必考虑行和列的相消,而且代码实现也十分简洁。下面来比较一下共轭梯度法和高斯消元法在复杂度上的优势。

共轭梯度法的每一次迭代都要做一次矩阵向量乘法和一些向量内积的计算,复杂度为n^2,当做完n次迭代后,复杂度变成了n^3,而高斯消元法只是\frac{1}{3}n^3左右,从这一点上,当矩阵不是稀疏矩阵的时候,共轭梯度法没有优势,而当矩阵变得稀疏的时候,n大的惊人,高斯消元法如果要得到解,就必须要做完所有的运算才能得到解,这样需要消耗大量的资源。而共轭梯度法不一样,它每一次迭代都能在某个分量上得到解,并且可以通过残差来度量解的精确情况,我们并不需要将算法进行到底,只需要解达到精度要求就可以退出。

同样的,与其他迭代法相比,共轭梯度法又能在确定的步数内收敛,这就是共轭梯度法的优势。

但是,当矩阵变成病态矩阵的时候,由于每一步的误差的累计,很有可能导致方向向量出现偏差,从而出现很糟糕的结果,这也是共轭梯度法的一个缺陷,我们可以通过预条件处理来减小病态矩阵带来的误差。

预条件处理

预条件处理的思想是降低方程系数矩阵的条件数,方法是左乘一个矩阵。

M^{-1}Ax=M^{-1}b

其中M是可逆的n阶矩阵,称为预条件子

常用的预条件子有:

  • 雅可比预条件子:M=D,DA的对角矩阵
  • 高斯-塞尔德预条件子:M=(D+\omega L)D^{-1}(D+\omega U),其中A=L+D+U,分别是下三角、对角、上三角矩阵,\omega是介于0和2之间的常数
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 201,924评论 5 474
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,781评论 2 378
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 148,813评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,264评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,273评论 5 363
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,383评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,800评论 3 393
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,482评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,673评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,497评论 2 318
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,545评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,240评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,802评论 3 304
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,866评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,101评论 1 258
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,673评论 2 348
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,245评论 2 341

推荐阅读更多精彩内容