zoukankan      html  css  js  c++  java
  • 主成分分析法

    主成分分析法


    主成分分析法:(Principle Component Analysis, PCA),是一个非监督机器学习算法,主要用于数据降维,通过降维,可以发现便于人们理解的特征,其他应用:可视化和去噪等。

    一、主成分分析的理解

    mark

    ​ 先假设用数据的两个特征画出散点图,如果我们只保留特征1或者只保留特征2。那么此时就有一个问题,留个哪个特征比较好呢?

    mark

    ​ 通过上面对两个特征的映射结果可以发现保留特征1比较好,因为保留特征1,当把所有的点映射到x轴上以后,点和点之间的距离相对较大,也就是说,拥有更高的可区分度,同时还保留着部分映射之前的空间信息。那么如果把点都映射到y轴上,发现点与点距离更近了,这不符合数据原来的空间分布。所以保留特征1相比保留特征2更加合适,但是这是最好的方案吗?

    1565281000485

    ​ 也就是说,我们需要找到让这个样本间距最大的轴?那么如何定义样本之间的间距呢?一般我们会使用方差(Variance),Var(x)=frac{1}{m}sum_{i=1}m(x_{i}-overline{x})2,找到一个轴,使得样本空间的所有点映射到这个轴的方差最大。

    1. 第一步:将样本进行均值归0,此时$Var(x)=frac{1}{m}sum_{i=1}m(x_{i})2$,此时计算过程就少一项,这就是为什么要进行样本均值归0。

    mark

    1. 第二步:需要定义一个轴的方向w=(w1, w2),使得我们的样本,映射到w以后,有:

    Var(X_{project}) = frac{1}{m}sum_{i=1}m(X_{project}(i) - overline{x}{project})^2最大。其实括号中的部分是一个向量,更加准确的描述应该是Var(X{project}) = frac{1}{m}sum_{i=1}mlVertX_{project}(i) - overline{x}{project}lVert^2,因为前面已经去均值,所以,这里只需Var(X{project}) = frac{1}{m}sum_{i=1}mlVertX_{project}(i) lVert^2最大

    1565281794841

    X^{(i)} cdot w = lVert X^{(i)} lVert cdot lVert w lVert cdot cos heta

    设w为单位向量,X^{(i)} cdot w = lVert X^{(i)} lVert cdot cos heta

    最终,X^{(i)} cdot w = lVert X_{project}^{(i)} lVert

    所以,最终只需Var(X_{project}) = frac{1}{m}sum_{i=1}^mlVert X^{(i)} cdot w lVert^2最大

    mark

    目标:求w,使得Var(X_{project}) = frac{1}{m}sum_{i=1}^m( X^{(i)} cdot w )^2最大

    Var(X_{project}) = frac{1}{m}sum_{i=1}^m( X_{1}^{(i)}w_{1} +X_{2}{(i)}w_{2}+dots+X_{n}{(i)}w_{n})^2

    Var(X_{project}) = frac{1}{m}sum_{i=1}^m( X^{(i)} cdot w )^2

    一个目标函数优化问题,使用梯度上升法解决。

    二、使用梯度上升法求解PCA

    mark

    mark mark

    首先生成一组数据:

    import numpy as np
    import matplotlib.pyplot as plt
    
    X = np.empty((100, 2))
    X[:, 0] = np.random.uniform(0., 100., size=100)
    X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0., 10., size=100) 
    
    plt.scatter(X[:, 0], X[:, 1])
    plt.show()
    
    def demean(X):
        # 不使用standardscale标准化数据,求解PCA需要方差,只去均值。
        return X - np.mean(X, axis=0) # axis=0表示最终求得的均值是一个行向量,也就是说将每一列求和 
    x_demean = demean(X)
    plt.scatter(x_demean[:,0], x_demean[:,1])
    plt.show()
    np.mean(x_demean[:,0])
    np.mean(x_demean[:,1])
    
    def f(w ,x):
        return np.sum((x.dot(w) ** 2)) / len(x)
    
    def df_math(w, x):
        return x.T.dot(x.dot(w)) * 2. / len(x)
    
    def df_denug(w, x, epsilon=0.0001):
        res = np.empty(len(w))
        for i in range(len(w)):
            w_1 = w.copy()
            w_1[i] += epsilon
            w_2 = w.copy()
            w_2[i] -= epsilon
            res[i] = (f(w_1, x) - f(w_2, x)) / (2 * epsilon)
        return res
    
    def direction(w):
        return w / np.linalg.norm(w)
    
    def gradient_ascent(df, x, init_w, eta, n_iters=1e-4, epsilon=1e-8):
        
        w = direction(init_w)
        cur_iter = 0
        while cur_iter < n_iters:
            gradient = df(w, x)
            last_w = w
            w = w + eta * gradient
            w = direction(w)
            if (abs(f(w, x) - f(last_w, x)) < epsilon):
                break
            cur_iter += 1
        return w
    

    测试:

    init_w = np.random.random(X.shape[1])    # 不能0向量开始
    init_w
    eta = 0.001
    gradient_ascent(df_denug, x_demean, init_w, eta)
    

    输出结果:array([0.84612666, 0.53298187])

    w = gradient_ascent(df_math, x_demean, init_w, eta)
    plt.scatter(x_demean[:, 0], x_demean[:,1])
    plt.plot([0, w[0]*30], [0, w[1]*30], color='r')
    plt.show()
    

    这只是一个从2维降到1维的情况,但在实际应用中往往都不这么简单。

    三、求数据的前n个主成分

    求出第一主成分以后,如何求出下一个主成分?数据进行改变,将数据在第一个主成分上的分量去掉。然后在新的数据上求第一主成分。

    mark

    第一步:先求出第一主成分

    import numpy as np
    import matplotlib.pyplot as plt
    
    X = np.empty((100, 2))
    X[:, 0] = np.random.uniform(0., 100., size=100)
    X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0., 10., size=100) 
    
    def demean(X):
        return X - np.mean(X, axis=0)
    
    X = demean(X)
    plt.scatter(X[:,0], X[:,1])
    plt.show()
    
    def f(w ,x):
        return np.sum((x.dot(w) ** 2)) / len(x)
    
    def df(w, x):
        return x.T.dot(x.dot(w)) * 2. / len(x)
    
    def direction(w):
        return w / np.linalg.norm(w)
    
    def first_component(df, x, init_w, eta, n_iters=1e-4, epsilon=1e-8):
        
        w = direction(init_w)
        cur_iter = 0
         
        while cur_iter < n_iters:
            gradient = df(w, x)
            last_w = w
            w = w + eta * gradient
            w = direction(w)
            if (abs(f(w, x) - f(last_w, x)) < epsilon):
                break
            cur_iter += 1
            
        return w
    
    init_w = np.random.random(X.shape[1])
    init_w
    eta = 0.01
    w = first_component(df, X, init_w, eta)
    w
    

    输出结果: array([0.71849226, 0.69553495])

    第二步:去掉第一主成分

    X2 = np.empty(X.shape)
    for i in range(len(X)):
        X2[i] = X[i] - np.dot(X[i], w) * w
    # X2 = X - X.dot(w).reshape(-1, 1) * w
    
    plt.scatter(X2[:,0], X2[:,1])
    plt.show()
    

    第三步:求第二主成分

    init_w2 = np.random.random(X2.shape[1])
    eta = 0.01
    w2 = first_component(df, X2, init_w2, eta)
    w2
    

    输出结果: array([ 0.98482183, -0.17356834])

    第四步:画出主成分

    plt.scatter(X2[:,0], X2[:,1])
    plt.plot([0, w[0]*30], [0, w[1]*30])
    plt.plot([0, w2[0]*30], [0, w2[1]*30])
    plt.show()
    

    mark

    import numpy as np
    import matplotlib.pyplot as plt
    
    X = np.empty((100, 2))
    X[:, 0] = np.random.uniform(0., 100., size=100)
    X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0., 10., size=100) 
    
    def demean(X):
        return X - np.mean(X, axis=0)
    
    X = demean(X)
    
    def direction(w):
        return w / np.linalg.norm(w)
    
    def df(w, x):
        return x.T.dot(x.dot(w)) * 2. / len(x)
    
    def f(w ,x):
        return np.sum((x.dot(w) ** 2)) / len(x)
    
    ef gradient_ascent(df, x, init_w, eta, n_iters=1e-4, epsilon=1e-8):
            
        w = direction(init_w)
        cur_iter = 0
        while cur_iter < n_iters:
            gradient = df(w, x)
            last_w = w
            w = w + eta * gradient
            w = direction(w)
            if (abs(f(w, x) - f(last_w, x)) < epsilon):
                break
            cur_iter += 1
        return w
    
    def fisrt_n_component(n, x, eta=0.01, n_iters=1e-4, epsilon=1e-8):
        x_pca = x.copy()
        x_pca = demean(x_pca)
        res = []
        for i in range(n):
            init_w = np.random.random(x_pca.shape[1])
            w = gradient_ascent(df, x_pca, init_w, eta)
            res.append(w)
            
            x_pca = x_pca - x_pca.dot(w).reshape(-1, 1) * w
            
        return res    
    
    fisrt_n_component(2, X)
    

    输出结果:[array([0.75978266, 0.65017714]), array([ 0.96416996, -0.26528531])]

    四、将高维数据向低维数据映射

    ​ 假设有数据矩阵X,通过主成分分析法得到前k个主成分对应的轴的单位方向向量之后,我们就可以将数据从高维向低维映射,同样也可以从低维在映射回高维空间。

    mark

    import numpy as np
    
    class PCA(object):
    
        def __init__(self, n_components):
            assert n_components >= 1, "n_component must be valid"
            self.n_components = n_components
            self.components_ = None
    
        def fit(self, x, eta=0.01, n_iters=1e-4):
            assert self.n_components <= x.shape[1], "n_components must be greater than the feature number of x"
    
            def demean(x):
                return x - np.mean(x, axis=0)
    
            def f(w, x):
                return np.sum((x.dot(w) ** 2)) / len(x)
    
            def df(w, x):
                return x.T.dot(x.dot(w)) * 2. / len(x)
    
            def direction(w):
                return w / np.linalg.norm(w)
    
            def gradient_ascent(df, x, init_w, eta, n_iters=1e-4, epsilon=1e-8):
    
                w = direction(init_w)
                cur_iter = 0
                while cur_iter < n_iters:
                    gradient = df(w, x)
                    last_w = w
                    w = w + eta * gradient
                    w = direction(w)
                    if (abs(f(w, x) - f(last_w, x)) < epsilon):
                        break
                    cur_iter += 1
                return w
    
            x_pca = demean(x)
            self.components_ = np.empty(shape=(self.n_components, X.shape[1]))
            for i in range(self.n_components):
                init_w = np.random.random(x_pca.shape[1])
                w = gradient_ascent(df, x_pca, init_w, eta, n_iters)
                self.components_[i, :] = w
    
                x_pca = x_pca - x_pca.dot(w).reshape(-1, 1) * w
    
            return self
    
        def transform(self, x):
            "将给定的x,映射到各个主成分分量中"
            assert x.shape[1] == self.components_.shape[1]
            return x.dot(self.components_.T)
    
        def inverse_transform(self, x):
            assert x.shape[1] == self.components_.shape[0]
            return x.dot(self.components_)
    
        def __repr__(self):
            return "PCA(n_components=%d)" % self.n_components
    
    if __name__ == '__main__':
        import numpy as np
        import matplotlib.pyplot as plt
    
        X = np.empty((100, 2))
        X[:, 0] = np.random.uniform(0., 100., size=100)
        X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0., 10., size=100)
    
        pca = PCA(n_components=1)
        pca.fit(X)
        print(pca.components_)
    
        x_reduction = pca.transform(X)
        print(x_reduction.shape)
    
        x_restore = pca.inverse_transform(x_reduction)
        print(x_restore.shape)
    
        plt.scatter(X[:, 0], X[:, 1], color='b', alpha=0.5)
        plt.scatter(x_restore[:, 0], x_restore[:, 1], color='r', alpha=0.5)
        plt.show()
    

    mark

    ​ 蓝色的点是随机初始生成,红色的点就是由pca降维之后,在从低维映射到高维的描述。其实完全可以只有一个轴来表示这些红色的点,也就完成了降维。

    五、scikit-learn中的PCA

    from sklearn.decomposition import PCA
    import numpy as np
    import matplotlib.pyplot as plt
    
    X = np.empty((100, 2))
    X[:, 0] = np.random.uniform(0., 100., size=100)
    X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0., 10., size=100) 
    
    def demean(X):
        return X - np.mean(X, axis=0)
    
    X = demean(X)
    
    pca = PCA(n_components=1)
    pca.fit(X)
    pca.components_
    

    输出结果:array([[-0.7896098, -0.6136093]])

    x_reduction = pca.transform(X)
    x_reduction.shape
    x_restore = pca.inverse_transform(x_reduction)
    x_restore.shape
    plt.scatter(X[:,0], X[:,1], color='b', alpha=0.5)
    plt.scatter(x_restore[:,0], x_restore[:,1], color='r', alpha=0.5)
    plt.show()
    

    mark

    接下来使用真实的数据集进行测试:

    import numpy as np
    import matplotlib.pyplot as plt
    from sklearn import datasets
    from sklearn.model_selection import train_test_split
    
    digits = datasets.load_digits()
    x = digits.data
    y = digits.target
    
    x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=0.2, random_state=666)
    x_train.shape
    

    接下来使用KNN先进行一个简单的测试:

    %%time
    from sklearn.neighbors import KNeighborsClassifier
    knn_clf = KNeighborsClassifier()
    knn_clf.fit(x_train, y_train)
    

    输出结果:Wall time: 4.88 ms

    knn_clf.score(x_test, y_test)
    

    输出结果:0.9888888888888889

    接下来使用降维:

    from sklearn.decomposition import PCA
    pca = PCA(n_components=2)
    pca.fit(x_train)
    #PCA(copy=True, iterated_power='auto', n_components=2, random_state=None,svd_solver='auto', tol=0.0, whiten=False)
    x_train_reduction = pca.transform(x_train)
    x_test_reduction = pca.transform(x_test)
    
    %%time
    knn_clf = KNeighborsClassifier()
    knn_clf.fit(x_train_reduction, y_train)
    

    输出结果:Wall time: 976 µs 相比降维之前速度快了很多。

    knn_clf.score(x_test_reduction, y_test)
    

    输出结果:0.6055555555555555,但是测试结果的准确率大打折扣。这是因为直接从64维降到2维,损失了太多信息。在sklearn中封装了一个函数,可以看出损失的详细信息。

    pca.explained_variance_ratio_
    

    输出结果:array([0.1450646 , 0.13714246]),可以看出这两个主成分分别损失了14.51%和13.71%,加和共损失28.22%。显然这不是想要的效果。

    那么我们查看一个64个主成分分别损失情况:

    pca = PCA(n_components=x_train.shape[1])
    pca.fit(x_train)
    pca.explained_variance_ratio_
    

    输出结果:

    array([1.45064600e-01, 1.37142456e-01, 1.19680004e-01, 8.43768923e-02,
           5.87005941e-02, 5.01797333e-02, 4.34065700e-02, 3.61375740e-02,
           3.39661991e-02, 3.00599249e-02, 2.38906921e-02, 2.29417581e-02,
           1.81335935e-02, 1.78403959e-02, 1.47411385e-02, 1.41290045e-02,
           1.29333094e-02, 1.25283166e-02, 1.01123057e-02, 9.08986879e-03,
           8.98365069e-03, 7.72299807e-03, 7.62541166e-03, 7.09954951e-03,
           6.96433125e-03, 5.84665284e-03, 5.77225779e-03, 5.07732970e-03,
           4.84364707e-03, 4.34595748e-03, 3.73352381e-03, 3.57655938e-03,
           3.30727680e-03, 3.18129431e-03, 3.06969704e-03, 2.89170006e-03,
           2.51205204e-03, 2.27743660e-03, 2.22760483e-03, 2.00065017e-03,
           1.89529684e-03, 1.56877138e-03, 1.42740894e-03, 1.39115781e-03,
           1.20896097e-03, 1.10149976e-03, 9.81702199e-04, 8.82376601e-04,
           5.69898729e-04, 4.10322729e-04, 2.32125043e-04, 8.49807543e-05,
           5.37426557e-05, 5.27990816e-05, 1.03398093e-05, 6.20749843e-06,
           5.03430895e-06, 1.15881302e-06, 1.09689390e-06, 5.51564753e-07,
           5.65215369e-08, 1.78867947e-33, 8.83490979e-34, 8.55393580e-34])
    

    从上面结果可以看出,64个主成分损失从大到小排列,最大14.5%,最小8.55393580e-34可以忽略不计。接下来使用曲线图绘制一下:

    plt.plot([i for i in range(x_train.shape[1])],
            [np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(x_train.shape[1])])
    plt.show()
    

    mark

    从图中,可以大概看出当降到30维的时候,保留信息大概在96%左右,因此可以考虑降到30为左右比较合适,既能保证精度又能保证速度。其实,sklearn中已经封装好了能够实现类似功能。

    pca = PCA(0.95)
    pca.fit(x_train)
    pca.n_components_
    

    输出结果:28,可以发现当损失信息为5%时,可以将数据从64维降到28维。

    x_train_reduction = pca.transform(x_train)
    x_test_reduction = pca.transform(x_test)
    
    %%time
    knn_clf = KNeighborsClassifier()
    knn_clf.fit(x_train_reduction, y_train)
    

    输出结果:Wall time: 2.93 ms

    knn_clf.score(x_test_reduction, y_test)
    

    输出结果:0.9833333333333333

    对比降维前后:

    1. 时间上,降维前:Wall time: 4.88 ms,降维后:Wall time: 2.93 ms
    2. 精度上,降维前:0.9888888888888889,降维后0.9833333333333333

    通过以上对比,可以发现在精度损失0.5%的情况下速度能提高一倍。这是一个比较的维度,如果维度很大,效果会更加优秀。

    在整个测试的过程中,曾尝试把数据降到2维,精度损失很大,但并不能说明这样做没有意义,通常情况下,我们会将高维数据降到3维或者2维进行可视化。

    pca = PCA(n_components=2)
    pca.fit(x)
    x_reduction = pca.transform(x)
    x_reduction.shape
    for i in range(10):
        plt.scatter(x_reduction[y==i, 0], x_reduction[y==i, 1], alpha=0.8)
    plt.show()
    

    mark

    通过对数据进行可视化,可以发现有蓝色、红色、紫色这些点在降到2维的情况下,依旧可以实现很高的识别度。这说明这些样本在原来的高维空间中差异就比较大。所以,一般会将数据进行可视化,进行一些简单的分析。

    六、对真实数据集MNIST使用PCA

    首先下载MNIST数据集的时候有个坑?

    import numpy as np
    from sklearn.datasets import fetch_mldata
    
    mnist = fetch_mldata("MNIST original", data_home='./datasets')
    

    如果直接运行会一直报错,反正就是下载不下来,但是会在datasets目录下生成一个文件夹,./datasets/mldata。所以需要下载数据集放在这个文件夹下,再次运行就可以了。

    数据集地址:

    链接:https://pan.baidu.com/s/15k6Sykf5TICN40KeEc6xgQ
    提取码:ujrq

    mnist
    
    {'DESCR': 'mldata.org dataset: mnist-original',
     'COL_NAMES': ['label', 'data'],
     'target': array([0., 0., 0., ..., 9., 9., 9.]),
     'data': array([[0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            ...,
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0],
            [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)}
    
    x, y = mnist['data'], mnist['target']
    x.shape
    y.shape
    

    输出结果:(70000, 784) (70000,)

    x_train = np.array(x[:60000], dtype=float)
    y_train = np.array(y[:60000], dtype=float)
    x_test = np.array(x[60000:], dtype=float)
    x_test = np.array(y[60000:], dtype=float)
    
    from sklearn.neighbors import KNeighborsClassifier
    
    knn_clf = KNeighborsClassifier()
    %time knn_clf.fit(x_train, y_train)
    

    输出结果:

    Wall time: 28 s
    
    KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
               metric_params=None, n_jobs=None, n_neighbors=5, p=2,
               weights='uniform')
    
    %time knn_clf.score(x_test, y_test)
    

    输出结果:

    Wall time: 11min 6s
    
    0.9688
    

    在这里,值得一提的是,之前讨论knn的时候,由于相比较样本之间距离,所以通常会对数据进行归一化。但是这里没用使用standardscale,因为这是图像。因为图像的像素点都在0-255,所以进行距离的比较是有意的,只有当尺度不一致时,才会考虑使用归一化。

    接下来使用PCA进行降维:

    from sklearn.decomposition import PCA
    
    pca = PCA(0.9)
    pca.fit(x_train)
    x_train_reduction = pca.transform(x_train)
    x_train.shape
    

    输出结果:(60000, 87) 可以发现通过降维实现了只有87维代替原来的784维,能够保留90的信息。

    from sklearn.neighbors import KNeighborsClassifier
    
    knn_clf = KNeighborsClassifier()
    %time knn_clf.fit(x_train_reduction, y_train)
    

    输出结果:Wall time: 436 ms 降维前使用了28s完成现在只需436ms。

    x_test_reduction = pca.transform(x_test)
    %time knn_clf.score(x_test_reduction, y_test)
    
    Wall time: 1min 17s
    0.9728
    

    将降维前后进行对比:

    1. 时间上:
      • 训练时间-降维前:28 s,降维后:436 ms
      • 测试时间-降维前:11min 6s,降维后:1min 17s
    2. 精度上:降维前:0.9688,降维后:0.9728

    经过上面的对比,发现经过降维之后速度提高了很多很多,而且准确率反而高了,这是为什么呢?这是因为其实降维还有另外一个作用,就是去噪

    七、使用PCA降噪

    import numpy as np
    import matplotlib.pyplot as plt
    
    X = np.empty((100, 2))
    X[:, 0] = np.random.uniform(0., 100., size=100)
    X[:, 1] = 0.75 * X[:, 0] + 3. + np.random.normal(0, 5, size=100)
    plt.scatter(X[:,0], X[:,1])
    plt.show()
    
    from sklearn.decomposition import PCA
    pca = PCA(n_components=1)
    pca.fit(X)
    x_reduction = pca.transform(X)
    x_restore = pca.inverse_transform(x_reduction)
    plt.scatter(x_restore[:, 0], x_restore[:, 1])
    plt.show()
    

    mark

    mark

    通过这个例子,发现pca还可以降噪,将高维数据降低到低维数据,这个过程中也可以理解成降低了维度,丢失了信息,同时也去除了噪音。为了更加直观,接下来使用手写数字识别数据集。

    from sklearn import datasets
    
    digits = datasets.load_digits()
    x = digits.data
    y = digits.target
    # 手动添加噪声
    noisy_digits = x + np.random.normal(0, 4, size=x.shape)
    # 取出部分数据做可视化
    example_digits = noisy_digits[y==0,:][:10]
    for num in range(1, 10):
        x_num = noisy_digits[y==num,:][:10]
        example_digits = np.vstack([example_digits, x_num])
    example_digits.shape
    
    def plot_digits(data):
        fig, axes = plt.subplots(10, 10, figsize=(10,10), 
                                subplot_kw={'xticks':[], 'yticks':[]}, 
                                gridspec_kw=dict(hspace=0.1, wspace=0.1))
        for i, ax in enumerate(axes.flat):
            ax.imshow(data[i].reshape(8, 8),
                     cmap='binary', interpolation='nearest',
                     clim=(0, 16))
        plt.show()
        
     plot_digits(example_digits)
    

    mark

    接下来使用pca去噪,

    from sklearn.decomposition import PCA
    
    pca = PCA(0.5)
    pca.fit(noisy_digits)
    pca.n_components_
    

    输出结果:12,可以看出保留50%的信息,只需要12维。

    components = pca.transform(example_digits)
    filterd_digits = pca.inverse_transform(components)
    plot_digits(filterd_digits)
    

    mark

    通过这个例子,可以对pca降维有一个感性的认识。

    八、PCA与人脸识别

    mark

    假设X一张人脸图像,W是钱k个主成分组成的矩阵,其实可以把W看做X的一些比较明显的特征。这就是特征脸(Eigenface)。

    首先下载数据集,sklearn中的lfw_people数据集:

    import numpy as np
    import matplotlib.pyplot as plt
    from sklearn.datasets import fetch_lfw_people
    
    faces = fetch_lfw_people()
    

    这里依旧有个小坑,就是数据下载不下来,不过没事,我们可以手动下载,复制这个网址 https://ndownloader.figshare.com/files/5976015 到迅雷中,下载到完整的lfwfunneded.tgz文件。并把这个文件放在C:UsersMogulscikit_learn_datalfw_home路径下,解压。然后再次运行就可以了。

    faces
    faces.keys()
    faces.data.shape
    #(13233, 2914),数据集中一共13233个样本,每个样本维度2914
    faces.images.shape
    # (13233, 62, 47),每个图片的大小(62, 47)
    random_indexes = np.random.permutation(len(faces.data))
    # 将数据集的顺序打乱
    X = faces.data[random_indexes]
    # 从中取除前36个
    example_faces = X[:36, :]
    example_faces.shape
    
    def plot_faces(faces):
        fig, axes = plt.subplots(6, 6, figsize=(10, 10),
                                subplot_kw={'xticks':[], 'yticks':[]},
                                gridspec_kw=dict(hspace=0.1, wspace=0.1))
        for i, ax in enumerate(axes.flat):
            ax.imshow(faces[i].reshape(62, 47), cmap='bone')
            
        plt.show()
    plot_faces(example_faces)
    

    mark

    faces.target_names
    # array(['AJ Cook', 'AJ Lamas', 'Aaron Eckhart', ..., 'Zumrati Juma',
           'Zurab Tsereteli', 'Zydrunas Ilgauskas'], dtype='<U35')
    len(faces.target_names)
    # 5749,数据集中一共有5749个人物,一共有13233张图片,差不多每个人2张。
    # 但实际情况并不是这样,因为样本并不均衡,有些可能只有一张,有些可能很多张。
    
    %%time
    from sklearn.decomposition import PCA
    pca = PCA(svd_solver='randomized')
    pca.fit(X)
    

    输出结果:Wall time: 18.4 s

    pca.components_.shape
    # (2914, 2914)
    # 取出前36个特征脸进行可视化
    plot_faces(pca.components_[:36, :])
    

    mark

    刚刚提到这个数据集的样本分布是不均衡的,那么接下来看看具体:

    faces2 = fetch_lfw_people(min_faces_per_person=60)
    faces2.data.shape
    # (1348, 2914)
    faces2.target_names
    

    输出结果:

    array(['Ariel Sharon', 'Colin Powell', 'Donald Rumsfeld', 'George W Bush',
           'Gerhard Schroeder', 'Hugo Chavez', 'Junichiro Koizumi',
           'Tony Blair'], dtype='<U17')
    

    根据输出结果可以发现,最少有60张图像的数据集一共有1348张图片,只有8个人的图像大于60,使用这样的数据进行人脸识别相对就比较合理了。

    我是尾巴:

    其实在PCA的背后是强大的数学支撑。在这里使用的梯度上升法对PCA进行求解,其实它是可以通过数学推导出相应的公式,具体的数学解释,以后会在写一篇。

    本次推荐:

    轻松玩转PDF

    你是如何强迫自己不断学习提升的?

    看过更大的世界,更优秀的人,就再也不甘心留在原地,不甘心就是动力!

    继续加油~

  • 相关阅读:
    AB(ApacheBench)工具 -- 压力测试
    请求错误 --405
    第三方支付 -----支付宝支付流程
    前后端分离 ---购物车
    canvas实现碰壁反弹(单个小方块)
    h5 的canvas绘制基本图形
    for in 循环
    js创建对象的三种方式和js工厂模式创建对象
    类的继承
    jq操作class类
  • 原文地址:https://www.cnblogs.com/zhangkanghui/p/11331276.html
Copyright © 2011-2022 走看看