机器学习算法之PCA算法

前言

在机器学习中降维是我们经常需要用到的算法,在降维的众多方法中PCA无疑是最经典的机器学习算法之一,最近准备撸一个人脸识别算法,也会频繁用到PCA,本文就带着大家一起来学习PCA算法。

前置内容

要学会PCA算法,首先需要了解矩阵分解算法。而矩阵分解算法又分为特征值分解和SVD(奇异值)分解,这两个算法的目的都是提取出一个矩阵最重要的特征。

特征值分解

特征值,特征向量

如果一个向量v是矩阵A的特征向量,则一定可以表示为下面的形式:
Av=\lambda v
其中\lambda是特征向量v对应的特征值,一个矩阵的一组特征向量是一组正交向量。这里大概可以理解为矩阵A与向量v相乘实际上就是对向量v进行了一次线性变换(旋转或者拉伸),而该变换的效果为常数\lambda乘以向量v。当我们求特征值和特征向量的时候就是为了求矩阵A能使哪些向量(特征向量)只发生伸缩变换,而变换的成都可以用特征值\lambda表示。

特征值与特征向量的几何意义

上面提到一个矩阵乘以一个向量就相当于一个线性变化可能难以理解,我们举个例子:
M = \left[ \begin{matrix} 3 & 0\\ 0 & 1 \\ \end{matrix} \right]
然后向量X=\left[ \begin{matrix} x\\ y \\ \end{matrix} \right]
所以M乘以一个向量(x,y)就表示为:
\left[ \begin{matrix} 3 & 0\\ 0 & 1\\ \end{matrix} \right] * \left[ \begin{matrix} x\\ y \\ \end{matrix} \right]=\left[ \begin{matrix} 3x\\ y \\ \end{matrix} \right]

我们在平面上化出这个变换过程如下:

在这里插入图片描述
同时,注意到上面的矩阵M是对称的,如果不是对称的,例如
那么M*X的变换就可以用下图来表示:
在这里插入图片描述
其中蓝色箭头指的是一个最主要的变换方向。如果要描述好一个变换,我们就需要描述好这个变换最主要的变化方向。

特征值分解

对于矩阵A,有一组特征向量v,将这组向量进行正交化单位化,就能得到一组正交单位向量。特征值分解,就是将矩阵A分解为如下式:
A=Q\sum Q^{-1}
其中,Q是矩阵A的特征向量组成的矩阵,\sum则是一个对角阵,对角线上的元素就是特征值。我们来分析一下特征值分解的式子,分解得到的\sum矩阵是一个对角矩阵,里面的特征值是由大到小排列的,这些特征值所对应的特征向量就是描述这个矩阵变换方向(从主要的变化到次要的变化排列)。
当矩阵是高维的情况下,那么这个矩阵就是高维空间下的一个线性变换,这个线性变换可能没法通过图片来表示,但是可以想象,这个变换也同样有很多的变化方向,我们通过特征值分解得到的前N个特征向量,就对应了这个矩阵最主要的N个变化方向。我们利用这前N个变化方向,就可以近似这个矩阵变换。也就是之前说的:提取这个矩阵最重要的特征

特征值分解举例

这里我们用一个简单的方阵来说明特征值分解的步骤。我们的方阵A定义为:

A = \begin{pmatrix} -1&1&0\\ -4&3&0\\ 1&0&2 \end{pmatrix}
首先,由方阵A的特征方程,求出特征值。
|A-\lambda E|= \begin{vmatrix} -1-\lambda & 1&0 \\ -4 & 3-\lambda &0 \\ 1&0&2-\lambda \end{vmatrix}=(2-\lambda)\begin{vmatrix} -1-\lambda & 1 \\ -4 & 3-\lambda \end{vmatrix}=(2-\lambda)(\lambda-1)^2=0
解方程得\lambda=2,1(重数为2)。
然后,把每个特征值\lambda代入到线性方程组(A-\lambda E)x=0里面,求出特征向量。
\lambda=2时,解线性方程组(A-2E)x=0
(A-2E)=\begin{pmatrix} -3&1&0\\ -4&1&0\\ 1&0&0 \end{pmatrix}->\begin{pmatrix} 1&0&0\\ 0&1&0\\ 0&0&0 \end{pmatrix},解得x_1=0,x_2=0。特征向量为:
p1=\begin{pmatrix} 0\\ 0\\ 1 \end{pmatrix}

同理,当\lambda=1时,解线性方程组(A-E)x=0
(A-E)=\begin{pmatrix} -2&1&0\\ -4&2&0\\ 1&0&1 \end{pmatrix}->\begin{pmatrix} 1&0&1\\ 0&1&2\\ 0&0&0 \end{pmatrix}
x_1+x_3=0,x_2+2x_3=0,特征向量为:p2=\begin{pmatrix} -1\\ -2\\ 1 \end{pmatrix}
最后,方阵A的特征值分解为:

A=Q\sum Q^{-1}=\begin{pmatrix} 0&-1&-1\\ 0&-2&-2\\ 1&1&1 \end{pmatrix}\begin{pmatrix} 2&0&0\\ 0&1&0\\ 0&0&1 \end{pmatrix}\begin{pmatrix} 0&-1&-1\\ 0&-2&-2\\ 1&1&1 \end{pmatrix}^{-1}

奇异值分解

上面讲解的特征值分解在实际应用的时候有一个最致命的缺点,就是只能用于方阵,也即是n*n的矩阵,而我们实际应用中要分解的矩阵大多数都不是方阵,所以SVD出现了。

什么是奇异值分解

奇异值分解是一个能适用于任意矩阵的一种分解的方法,对于任意矩阵A总是存在一个奇异值分解:
A=U\sum V^T
假设A是一个m*n的矩阵,那么得到的U是一个m*m的方阵,U里面的正交向量被称为左奇异向量。\sum是一个m*n的矩阵,\sum除了对角线其它元素都为0,对角线上的元素称为奇异值。V^T是V的转置矩阵,是一个n*n的矩阵,它里面的正交向量被称为右奇异值向量。而且一般来讲,我们会将\sum上的值按从大到小的顺序排列。
从上面的奇异值分解公式来看,我们是不知道如何拆分矩阵A的。但我们可以想办法把奇异值和特征值联系起来。首先,我们用矩阵A的转置乘以矩阵A,得到一个方阵,用这个方阵进行特征值分解,得到的特征值和特征向量满足下面的等式:
(A^TA)v_i=\lambda_iv_i
这里的v_i就是我们要求的右奇异向量。
其次,我们将A和A的转置做矩阵的乘法,得到一个方阵,用这样的方阵进行特征分解,得到的特征和特征向量满足下面的等式:
(AA^T) \mu _i = \lambda_i \mu_i
这里的\mu_i就是左奇异向量。

上面直接给出了\mu_i是左奇异向量,v_i是右奇异向量,根据是什么呢?我们尝试以V矩阵举例来证明一下:
A=U\sum V^T=>A^T=V\sum^TU^T=>A^TA=V\sum^TU^TU\sum V^T=V\sum^2V^T(*)
上式证明种使用了U^TU=I,\sum^T\sum=\sum^2。可以看出,A^TA的特征向量组成的矩阵就是我们SVD种的V矩阵,而AA^T得特征矩阵就是我们SVD中的U矩阵。
所以,我们可以得到奇异值的两种求法:

  • 第一种:
    A=U\sum V^T=>AV=U\sum V^TV=>AV=U\sum=>Av_i=\delta_iu_i=>\delta_i=\frac{A_vi}{u_i}
  • 第二种:
    通过上面*式的证明,我们还可以看出,特征值矩阵等于奇异值矩阵的平方,也就是说特征值和奇异值满足如下关系:

\delta_i=\sqrt{\lambda_i}
其中,\sigma_i就是奇异值,奇异值\sigma_i和特征值相似,在矩阵\sigma种也是从大到小排列。
上面了我们将了奇异值分解的原理,那么我们就可以直接按照上面的公式来计算了吗?
考虑一下直接使用上面的公式来计算,我们可以计算一下算法的复杂度,我们将m*n的矩阵分解为m\times m的矩阵U,m\times n的矩阵\sumn*n的矩阵V^T,这个矩阵的复杂度大概是O(N^3)级别的,所以在求奇异值的时候需要耗费非常多的时间,整个算法的开销就过大了。
幸运的是,当我们将分解矩阵\sum种的奇异值按照从大到小的顺序排列之后,奇异值\delta_i从大到小的顺序减小的特别快。在很多情况下,前10%甚至1%的奇异值的和就占了全部的奇异值之和的99%以上。也就是说,剩下的90%甚至99%的奇异值几乎没有什么作用。因此,我们可以用前面r个大的奇异值来近似描述矩阵,于是奇异值分解公式可以写成如下:

A_{m*n}\approx U_{m*r}\sum_{r*r}V_{r*n}^T

其中r是一个远远小于m和n的数,右边的三个矩阵的相乘结果将会得到一个接近A的矩阵。如果r的取值远远小于n,从计算机内存的角度来说,右边三个矩阵的存储内存要远远小于矩阵A的。所以在奇异值分解中r的取值很重要,就是在计算精度和时间空间之间做选择。

奇异值分解举例

定义矩阵A为:
A=\begin{pmatrix} 0&1\\ 1&1\\ 1&0 \end{pmatrix}
首先,我们先求出A^TAAA^T
A^TA=\begin{pmatrix} 0&1&1\\ 1&1&0\\ \end{pmatrix}*\begin{pmatrix} 0&1\\ 1&1\\ 1&0 \end{pmatrix}=\begin{pmatrix} 2&1\\ 1&2\\ \end{pmatrix}

AA^T=\begin{pmatrix} 0&1\\ 1&1\\ 1&0 \end{pmatrix}*\begin{pmatrix} 0&1&1\\ 1&1&0\\ \end{pmatrix}=\begin{pmatrix} 1&1&0\\ 1&2&1\\ 0&1&1 \end{pmatrix}
然后,求出A^TAAA^T的特征值和特征向量:
A^TA的特征向量和特征值:
\lambda_1=3, v_1=\begin{pmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}
\lambda_2=1, v_2=\begin{pmatrix} \frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}
AA^T的特征向量和特征值:

\lambda_1=3, u_1=\begin{pmatrix} \frac{1}{\sqrt{6}}\\ \frac{2}{\sqrt{6}}\\ \frac{1}{\sqrt{6}} \end{pmatrix}
\lambda_2=1, u_2=\begin{pmatrix} \frac{1}{\sqrt{2}}\\ 0\\ \frac{-1}{\sqrt{2}}\\ \end{pmatrix}
\lambda_3=0, u_3=\begin{pmatrix} \frac{1}{\sqrt{3}}\\ \frac{-1}{\sqrt{3}}\\ \frac{1}{\sqrt{3}}\\ \end{pmatrix}
然后,我们再利用Av_i=\sigma_i\mu_i,i=1,2求奇异值:
\begin{pmatrix} 0&1\\ 1&1\\ 1&0\\ \end{pmatrix}*\begin{pmatrix} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}=\sigma_1\begin{pmatrix} \frac{1}{\sqrt{6}}\\ \frac{2}{\sqrt{6}}\\ \frac{1}{\sqrt{6}} \end{pmatrix}=>\sigma_1=\sqrt{3}
\begin{pmatrix} 0&1\\ 1&1\\ 1&0\\ \end{pmatrix}*\begin{pmatrix} \frac{-1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{pmatrix}=\sigma_2\begin{pmatrix} \frac{1}{\sqrt{2}}\\ 0\\ \frac{-1}{\sqrt{2}} \end{pmatrix}=>\sigma_2=1

这一步也可以直接用\sigma_i=\sqrt{\lambda_i}求出奇异值\sqrt{3}和1。
最后,我们得到A的奇异值分解为:
A=U\sum V^T=\begin{pmatrix} \frac{1}{\sqrt{6}}&\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{3}}\\ \frac{2}{\sqrt{6}}&0&\frac{-1}{\sqrt{3}}\\ \frac{1}{\sqrt{6}}&\frac{-1}{\sqrt{2}}&\frac{-1}{\sqrt{3}}\\ \end{pmatrix}*\begin{pmatrix} \sqrt{3}&0\\ 0&1\\ 0&0 \end{pmatrix}*\begin{pmatrix} \frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\ \frac{-1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\ \end{pmatrix}

容易发现,经过SVD分解后,可以用前r个非零奇异值对应的奇异向量表示矩阵A的主要特征,这样就把矩阵A进行了降维。

协方差和散度矩阵

样本均值:
\hat{x}=\frac{1}{n}\sum_{i=1}^Nx_i
样本方差:
S^2=\frac{1}{n-1}\sum_{i=1}^n(x_i-\hat{x})^2
样本X和样本Y的协方差:
Cov(X,Y)=E[(X-E(X))(Y-E(Y))]=\frac{1}{n-1}\sum_{i=1}^n(x_i-\hat{x})(y_i-\hat{y})
由上面的公式,我们可以得到以下结论:

(1) 方差的计算公式是针对一维特征,即针对同一特征不同样本的取值来进行计算得到;而协方差则必须要求至少满足二维特征;方差是协方差的特殊情况。

(2) 方差和协方差的除数是n-1,这是为了得到方差和协方差的无偏估计。
协方差为正时,说明X和Y是正相关关系;协方差为负时,说明X和Y是负相关关系;协方差为0时,说明X和Y是相互独立。Cov(X,X)就是X的方差。当样本是n维数据时,它们的协方差实际上是协方差矩阵(对称方阵)。例如,对于3维数据(x,y,z),计算它的协方差就是:

在这里插入图片描述
散度矩阵定义为:
在这里插入图片描述
对于数据的散度矩阵为。协方差矩阵和散度矩阵关系密切,散度矩阵就是协方差矩阵乘以。因此它们的特征值和特征向量是一样的。这里值得注意的是散度矩阵就是上面讲的SVD分解的第一步,因此PCA和SVD也是有相关性的。

PCA算法

PCA即(Principal Component Analysis)主成分分析算法,是机器学习种应用得最广泛的数据降维算法。PCA的思想是将原始n维的数据映射到k维上(k<n),这k维是全新的正交特征,也叫主成分。PCA的工作就是在原始的数据空间种顺序的找一组相互正交的坐标轴,新的坐标轴和数据本身是密切相关的。其中第一个坐标轴是原始数据中方差最大的方向,第二个坐标轴是和第一个坐标轴相交的坐标轴种最大的,以此内推,k个坐标轴是完全正交的。研究发现,大部分方差都包含在k个坐标轴种,后面的坐标轴所含的方差几乎维0。所以可忽略不计,以实现对数据的降维处理。

PCA算法实现

基于特征值分解协方差矩阵实现PCA算法

输入数据集X={x_1,x_2,x_3,...,x_n},需要降维到k维。
1)去均值,即将每一维特征减掉各自的平均值。
2)计算协方差矩阵\frac{1}{n}XX^T,注:里除或不除样本数量n或n-1,其实对求出的特征向量没有影响。
3)用特征值分解方法求协方差矩阵\frac{1}{n}XX^T的特征值与特征向量。
4)对特征值从大到小排序,选择其中最大的k个。然后将其对应的k个特征向量分别作为行向量组成特征向量矩阵P。
5)将数据转换到k个特征向量构建的新空间中,即Y=PX。

举个例子:
我们要将之前提到的矩阵A=\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}降到1维。
因为这个矩阵的每行已经是零均值,这里我们直接求协方差矩阵:
C=\frac{1}{5}\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}\begin{pmatrix} -1 & -2 \\ -1 & 0 \\ 0 & 0 \\ 2 & 1 \\ 0 & 1 \end{pmatrix}=\begin{pmatrix} \frac{6}{5} & \frac{4}{5} \\ \frac{4}{5} & \frac{6}{5} \end{pmatrix}
然后求其特征值和特征向量,得到:
\lambda_1=2,\lambda_2=2/5
其对应的特征向量分别是:
c_1\begin{pmatrix} 1 \\ 1 \end{pmatrix},c_2\begin{pmatrix} -1 \\ 1 \end{pmatrix}
其中对应的特征向量分别是一个通解,c1c2可以取任意实数。那么标准化后的特征向量为:
\begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix},\begin{pmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix}
因此我们的矩阵P是:
P=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}
可以验证协方差矩阵C的对角化:
PCP^\mathsf{T}=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\begin{pmatrix} 6/5 & 4/5 \\ 4/5 & 6/5 \end{pmatrix}\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}=\begin{pmatrix} 2 & 0 \\ 0 & 2/5 \end{pmatrix}
最后我们用P的第一行乘以数据矩阵,就得到了降维后的表示:
Y=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}=\begin{pmatrix} -3/\sqrt{2} & -1/\sqrt{2} & 0 & 3/\sqrt{2} & -1/\sqrt{2} \end{pmatrix}
数据矩阵X降维投影结果为:

在这里插入图片描述

基于SVD实现PCA算法

输入数据集X={x_1,x_2,x_3,...,x_n},需要降维到k维。

  1. 去平均值,即每一位特征减去各自的平均值。
  2. 计算协方差矩阵。
  3. 通过SVD计算协方差矩阵的特征值与特征向量。
  4. 对特征值从大到小排序,选择其中最大的k个。然后将其对应的k个特征向量分别作为列向量组成特征向量矩阵。
  5. 将数据转换到k个特征向量构建的新空间中。
    和利用特征值分解实现PCA算法类似,我们需要找到样本协方差矩阵XX^T的最大k个特征向量,然后用这最大的k个特征向量组成的矩阵来做低维投影降维。可以看出,在这个过程中需要先求出协方差矩阵XX^T当样本数多、样本特征数也多的时候,这个计算还是很大的。当我们用到SVD分解协方差矩阵的时候,SVD有两个好处:
  • 有一些SVD的实现算法可以先不求出协方差矩阵XX^T也能求出我们的右奇异矩阵V。也就是说,我们的PCA算法可以不用做特征分解而是通过SVD来完成,这个方法在样本量很大的时候很有效。实际上,scikit-learn的PCA算法的背后真正的实现就是用的SVD,而不是特征值分解。
  • 注意到PCA仅仅使用了我们SVD的左奇异矩阵,没有使用到右奇异值矩阵,那么右奇异值矩阵有什么用呢?
    假设我们的样本是m\times n的矩阵X,如果我们通过SVD找到了矩阵X^TX最大的k个特征向量组成的k*n的矩阵V^T,则我们可以做如下处理:
    X'_{m*k}=X_{m*n}V_{n*k}^T
    可以得到一个m\times k的矩阵X',这个矩阵和我们原来m\times n的矩阵X相比,列数从n减到了k,可见对列数进行了压缩。也就是说,左奇异矩阵可以用于对行数的压缩;右奇异矩阵可以用于对列(即特征维度)的压缩。这就是我们用SVD分解协方差矩阵实现PCA可以得到两个方向的PCA降维(即行和列两个方向)。

代码实现

import numpy as np
# 零均值化,即中心化,是数据的预处理方法
def zero_centered(data):
    matrix_mean = np.mean(data, axis=0)
    return data - matrix_mean

def pca_eig(data, n):
    new_data = zero_centered(data)
    conv_mat = np.dot(new_data.T, new_data) #也可以用np.cov()方法
    eig_values, eig_vectors = np.linalg.eig(np.mat(conv_mat))
    # 求特征值和特征向量,特征向量是列向量
    value_indices = np.argsort(eig_values) #将特征值从小到大排序
    n_vectors = eig_vectors[:, value_indices[-1: -(n+1): -1]]
    # 最大的n个特征值对应的特征向量
    return new_data * n_vectors #返回低维特征空间的数据

def pca_svd(data, n):
    new_data = zero_centered(data)
    cov_mat = np.dot(new_data.T, new_data)
    U, s, V = np.linalg.svd(cov_mat) #将协方差矩阵奇异值分解
    pc = np.dot(new_data, U) #返回矩阵的第一个列向量即是降维后的结果
    return pc[:, 0]

def unit_test():
    data = np.array([[2.5, 2.4], [0.5, 0.7], [2.2, 2.9], [1.9, 2.2], [3.1, 3.0], 
                     [2.3, 2.7], [2, 1.6], [1, 1.1], [1.5, 1.6], [1.1, 0.9]])
    result_eig = pca_eig(data, 1)
    # 使用常规的特征值分解法,将二维数据降到一维
    print(result_eig)
    result_svd = pca_svd(data, 1)
    # 使用奇异值分解法将协方差矩阵分解,得到降维结果
    print(result_svd)
unit_test()

得到的降维结果如下:


在这里插入图片描述

可以看到两种方法的降维结果是一致的。

参考文章

https://mp.weixin.qq.com/s/Dv51K8JETakIKe5dPBAPVg
https://blog.csdn.net/program_developer/article/details/80632779
http://blog.codinglabs.org/articles/pca-tutorial.html

欢迎关注我的微信公众号GiantPadaCV,期待和你一起交流机器学习,深度学习,图像算法,优化技术,比赛及日常生活等。


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

推荐阅读更多精彩内容