A tutorial on Principal Components Analysis | 主成分分析（PCA）教程

A tutorial on Principal Components Analysis

原著：Lindsay I Smith, A tutorial on Principal Components Analysis, February 26, 2002.

翻譯：houchaoqun.
時間：2017/01/18.
出處：http://blog.csdn.net/houchaoqun_xmu  |  http://blog.csdn.net/Houchaoqun_XMU/article/details/54601984

說明：本文在翻譯原文的過程中，在保證不改變原文意思的前提下，對一些知識點進行了擴充並附上參考網址。本人剛開始嘗試翻譯文獻，經驗不足，有些理解可能不到位，有些翻譯可能不夠准確。若讀者們有發現有誤之處，還望海涵並加以指正，本人會及時加以指正。

• 本文中涉及到公式的編輯，使用的是“在線LaTeX公式編輯器”，將編輯完成后的公式（圖片的形式）復制到博文中即可。
• 原文可能缺少一些圖，本人通過Matlab作圖並補充。
• 本人在原文的基礎上，額外提供了一些參考網址，有興趣的讀者們可以詳細閱讀。
• 為了更讀者更深刻的理解PCA，本人還另外提供了案例，並用Matlab實現，僅供讀者們參考。
• 為了更好的將理論運用到實際，本人打算在介紹完這篇PCA教程后，整理一篇基於PCA進行人臉識別的具體案例，由於水平有限，需要改正和完善的地方希望讀者們提供寶貴的意見。

--------------------------------

主成分分析（PCA）教程

第1章：簡介（Introduction）

    本教程的目的是為了讓讀者理解主成分分析（PCA）。PCA是一項很有用的統計學技術，已經應用在人臉識別（Face Recognition）和圖像壓縮（Image Compression）等領域，並且作為一項通用技術在高維數據中發現數據模型（降維技術）。

    開始介紹PCA之前，本教程首先介紹一些學習PCA相關的數學概念。包括，標准差，協方差，特征向量及其對應的特征值。這部分背景知識旨在使得PCA相關的章節更加簡潔明了，如果讀者對這部分的概念很熟悉，可以選擇跳過這部分內容。

    本教程的一些相關案例旨在說明這些數學概念，讓讀者對這部分內容理解得更加深刻。如果讀者想要更進一步的學習，“初等線性代數 第五版（Elementary Linear Algebra 5e）”這本數學工具書是一個很好的關於這些數學背景的資源。這本書的作者是Howard Anton，由John Wiley和Sons股份有限公司出版（國際標准圖書編號[ISBN] 0-471-85223-6）。

 

第2章：數學背景（Background Mathematics）

    本章將給出一些對於理解PCA原理（涉及相關的推導過程）所需的基礎數學技能。這部分知識都是相互獨立的，並且都有相關的例子。理解為什么可能使用到這個技術以及在數據方面，某個操作的結果告訴了我們什么，比記住准確的數學技術更重要。雖然並不是這些技術都運用到PCA，但是最重要的技術（運用到PCA的技術）正是基於那些沒有明確需要的技術（沒有明確運用到PCA），所以對於初學者，還是有必要好好理解一下各個知識點。

    本章包含兩大部分，第一部分介紹了關於統計學的知識，側重於介紹數據（分布測量）是如何分布的（標准差，方差，協方差以及協方差矩陣）。第二部分介紹了矩陣代數，側重於特征向量和特征值，矩陣的一些重要性質是理解PCA的基礎。

2.1 統計學（Statistics）

    關於統計學的整個主題是圍繞着你有一個大的數據集，並且你想要分析這個數據集中點跟點之間的關系。本章將側重介紹一些讀者能夠親自操作（手算驗證）的測量指標（標准差，方差，協方差以及協方差矩陣），並且對於數據本身，這些測量指標告訴了我們什么（即，測量值反映數據的關系）。

2.1.1 標准差（Standard Deviation）

    為了理解標准差，我們需要一個數據集，統計學家通常使用來自總體（population）中的一個樣本（sample）。此處使用選舉投票（election polls）為例，“總體”是這個國家的所有人，而統計學家測量（measure）的樣本是這個總體的一個子集。統計學的偉大之處在於僅僅通過測量（可以通過電話調查或者類似的方法）總體中的一個樣本，就能夠算出與總體最相似的測量結果（即，通過測量樣本進而解決總體）。在統計學這個章節，作者假設本文使用的數據集是來自更大總體的樣本（總體的數據集比樣本的數據集更大）。（本章）下文有一個關於樣本和總體更詳細信息的參考網址。

    例如，對於數據集X：

    我們簡單的使用符號X表示整個數據集。如果想表示數據集中的某個元素，可以使用下標表示指定的元素（值），如   表示數據集X的第三個元素，對應的值是4；  指的是這個序列中的第一個元素，值為1；但是這里並沒有像我們在一些教科書上提到的 （因為此處的下標從1開始）。此外，符號n表示數據集X的元素個數。

    我們可以對一組數據集進行很多種計算。例如，我們可以計算這個樣本（數據集）的平均值（mean）。作者假設讀者們理解樣本X的均值的含義，此處僅給出計算公式，如下：

    注：符號   表示數據集X的均值。上述公式表示，“累加數據集X中的所有元素的值，然后除以元素總數n”。

   不幸的是，均值僅僅告訴我們數據集的中間點。例如，以下兩組數據集具有相關的均值（mean = 10），但是這是兩組完全不同的數據。

    因此，這兩組數據集有什么不同呢？數據的分布不同。數據集的標准差反映了數據是如何分布的（個體間的離散程度）。

    我們要怎樣計算標准差（Standard Deviation）呢？標准差（SD）的英文定義是，“The average distance from the mean of the data set to a point（數據集中每個點到均值的平均距離）”。標准差的計算方式是，首先計算每個點到樣本均值的距離的平方，然后進行累加，最后將累加的結果除以 n-1 再開根號，如下公式所示：

    其中，通常使用s表示一個樣本的標准差。讀者可能會問，“為什么是除以 (n-1) ，而不是 n？”。答案有些復雜，但總的來說，如果你的數據集是樣本數據集，也就是說如果你選擇的是一個真實世界中的子集（例如，對樣本數為500個人的選舉投票案例進行測量），你就需要使用 (n-1) 。因為事實證明，對於樣本，使用(n-1)得到的答案更接近於標准差的結果（如果你嘗試使用Matlab提供的cov函數，也是除以n-1）。如果你處理的數據集是總體，那你應該使用n。（ 其他解釋參考[ 點擊進入鏈接 ]：如是總體，標准差公式根號內除以n，如是樣本,標准差公式根號內除以（n-1)；這樣能使我們以較小的樣本集更好的逼近總體的標准差，即統計上所謂的“無偏估計”）。不管怎樣，如果你計算的不是樣本的標准差，而且總體的標准差，那么公式根號內應該除以n，而不是(n-1)。對於這部分內容，有意願更進一步了解的讀者可以參考網址： http://mathcentral.uregina.ca/RR/database/RR.09.95/weston2.html。這個網址用相似的方式介紹了標准差，還提供了一個實驗案例，解釋了標准差公式根號內的兩種分母（即，n和n-1）的區別，並對“樣本”和“總體”的區別進行了討論。

    如下表2.1所示的兩組數據，分別計算了它們的標准差。

   正如預期的，表2.1的第一部分的標准差（8.3266）更大是因為這組樣本數據相對於均值的分布更離散。正如另一個例子所示，樣本數據為[ 10 10 10 10 ]，均值同樣為10，但是它的標准差為0。因為所有元素的值都一樣，不存在偏離均值的元素。

2.1.2 方差（Variance）

    方差是另一種測量數據集離散程度的方式。實際上，方差和標准差幾乎是相同的。公式如下：

    讀者們可以發現方差僅僅是對標准差等式左右兩邊進行平方（方差的公式中沒有對其開根號）。對於樣本中的方差   是常用的符號。這兩個衡量值（方差和標准差）都是衡量數據的離散程度。標准差是最常用的衡量值，但方差同樣適用。本文在介紹了標准差之外，還介紹方差的原因是為下一個章節的協方差提供基礎（鋪墊，拋磚引玉）。

    練習：如下三組數據集所示，分別求出均值，標准差和方差。

• [ 12 23 34 44 59 70 98 ]
• [ 12 15 25 27 32 88 89 ]
• [ 15 35 78 82 90 96 97 ]

2.1.3 協方差（Covariance）

    上述我們關注的兩種（標准差和方差）測量指標的對象僅僅是一維數據集（局限性），比如房間里所有人的身高或者教師批改試卷等。然而，許多數據集都不只是一維，而且對於這些數據集，統計分析的目標通常是尋找這些維度之間是否存在任何關系。例如，我們將班級上所有學生的身高以及學生考試的分數作為數據集，然后使用統計分析方法去判斷學生的身高是否對他們的分數有任何影響。

    標准差和方差僅僅處理一維數據，因此你只能分別對數據集中的每一維計算出標准差（每一維對應一個標准差）。然而，找到一種與上述相似，並且能夠度量各個維度偏離其均值程度的測量指標是非常有用的。

    【協方差】就是這樣的一種測量指標。協方差通常用來測量2維的數據集。如果你計算的協方差中的2個維度都是同一個變量，此時得到的就是方差。（協方差正常處理的是二維數據，維數大於2的時候，就用協方差矩陣表示）因此，如果有一組三維數據集（x, y, z），那么我們需要分別計算x維度和y維度，x維度和z維度，y維度和z維度的協方差。計算x維度和x維度，y維度和y維度，z維度和z維度的協方差就相當於分別計算x，y，z維度的方差（這就是兩個維度都是同一個變量的情況）。

    協方差的公式和方差的公式很相似，方差的公式（度量各個維度偏離其均值的程度）如下所示：

    其中，作者僅僅是拓展了平方項變成兩個部分（如，  表示成  *  ）。因此，協方差的公式如下：

    上述兩個公式（方差和標准差）除了分子中的第二個括號內將X用Y代替以為，其他都一樣。用英文定義就是，“For each data item, multiply the different between the x value and the mean of x, by the difference between the y value and the mean of y. Add all these up, and divide by (n -1)”（對於X和Y的每個數據項，  與 （變量X的均值）的差乘以   與  （變量Y的均值）的差，然后將它們累加起來后除以 (n-1) ）。

    那么，協方差是如何運作的呢？讓我們使用一些實例數據。想象我們走進世界中，收集到一些二維數據。我們統計一群學生用在學習  COSC241 這個課程的時間（小時）以及他們獲得的分數。現在我們有了二維的數據，第一維是H維，表示學生學習的時間（以小時為單位），第二維是M維，表示學生獲得的分數。如下圖所示，包含了我們剛剛想象的數據以及H和M之間的協方差  的計算結果。

圖示：學生學習的時間H，獲得的分數M

    那么，這些數據告訴我們什么？協方差的符號（正好或者負號）比數值更加重要。如果協方差的值是正數，表明這兩維（H和M）是正相關，也就是說，隨着學生學習的時間（H）的增加，學生獲得的分數（M）也增加。如果值是負數，一個維度增加則另一個維度減小（數據呈負相關）。對於這個例子，如果協方差最終得到的結果是負值，也就是負相關，意味着隨着學生學習時間的增加，學生獲得的分數則減少。最后還有一種情況是，當協方差的值為0，意味着這兩維（如，變量H和M）是相互獨立的（不相關）。

    實驗結果（上述數據）可以很直觀的通過對原始數據作圖，如下圖所示，學生所得的分數M隨着學生學習的時間H的增加而增加。然而，數據的可視化（作圖）也僅僅局限於2維和3維的數據。因為對於任意2維的數據都可以計算協方差的值，因此這項技術（協方差）經常用於發現高維數據之間的關系（這些高維數據是很難進行可視化的）。

    讀者們或許會問，“   等於   嗎？”。通過協方差公式告訴我們，答案是“   等於   ”，（展開協方差公式）這兩者唯一的區別就是式子   用式子    替代。根據“乘法交換律”可知，兩個因數相乘，交換因數的位置，積不變，也就是說 和 得到相同的結果。

2.1.4 協方差矩陣（The Covariance Matrix）

    回顧協方差的相關知識，它通常測量的是二維數據。假如我們有一組大於2維的數據，那么將有多個協方差可以計算（有多種組合形式的協方差）。例如，對於一個三維的數據（x, y, z），我們能夠計算 ,   和  。實際上，對於一組n維數據，我們可以計算   種不同組合的協方差。一種有用的方式是，“計算不同維度之間所有可能情況的協方差值，然后將它們放入一個矩陣”，構成協方差矩陣。（本教程假設讀者們對矩陣很熟悉並且知道如何去定義矩陣）。因此，一組n維數組的協方差矩陣的定義如下：

其中， 表示一個n行n列的矩陣，  表示第x維。這一串復雜的公式（ugly looking formula）的意思是，如果你有一組n維的數據集，那么得到的是一個n行n列的協方差矩陣，矩陣中的每一個元素值是通過計算兩個不同維度的協方差得到的（矩陣的位置跟協方差的組合有對應關系）。例如，第2行，第3列位置的元素值是通過計算第2維和第3維之間的協方差得到的。

    舉一個例子，我們將虛構一組的3維數據集的協方差矩陣（使用x, y, z表示這三維）。得到一個3行3列的協方差矩陣，對應的元素值如下所示：

    一些關鍵的注意點：一個注意點是，沿着協方差矩陣的主對角線方向，你會發現元素值是某一維和它本身的協方差，也就是這個維的方差。另一個注意點是，由 可知，協方差矩陣是一個對稱陣。（協方差還有另一種定義形式，參考 百度百科）

   練習：

    如下所示的一組2維數據集，計算出x和y維之間的協方差，並根據計算結果說明一下數據。

    計算如下所示的一組3維數據集的協方差矩陣。

2.2 矩陣代數（Matrix Algebra）

    本章節提供矩陣代數的背景知識有助於下文對PCA的理解。作者將側重的介紹對於給定一個矩陣（方陣）的特征向量和特征值。同樣的，作者假設讀者們已經了解了一些矩陣的基礎知識。

2.2.1 特征向量（Eigenvectors）

   如讀者們所知，假如兩個矩陣的行數和列數滿足矩陣相乘的條件（例如，C=AB，當矩陣A的列數等於矩陣B的行數時，A與B可以相乘），就可以將他們相乘。特征向量就是其中的一個特例，考慮如下2個“矩陣乘以向量”的式子。

式1：沒有特征向量

式2：有一個特征向量

    在第一個乘法式子（式1）里，矩陣乘以向量的結果並不是形如一個整數乘以原始向量（沒有特征向量）。而在第二個乘法式子（式2）里相乘得到的結果是4乘以原始向量 。為什么是這樣子呢？原因如下：原始向量 是一個2維空間的向量，（在x-y坐標系下）向量 表示一個從原點（0, 0）指向（3, 2）的一個箭頭。式2左邊的另一個乘法因子是一個方陣（行數等於列數的矩陣）可以表示為一個變換矩陣。如果將這個變換矩陣左乘以（左乘和右乘有區別）一個向量，所得的結果是由原始狀態（向量）通過變換矩陣得到的另一個向量。

    特征向量來自於變換（矩陣）的性質（ 線性變換的特征向量是一個非簡並的向量，其方向在該變換下不變。該向量在此變換下縮放的比例稱為其特征值）。（A * V）想象有這么一個變換矩陣A作為左乘因子，將向量V映射到直線 y=x 上。然后你可以注意到是否有一個向量坐落於直線方程 y=x 上，這就是它自身的映射。此時，向量V（可以任意擴大它的倍數，因為向量的長度並不影響結果）就是變換矩陣A的一個特征向量。

    特征向量還有其他什么性質呢？讀者們首先應該知道特征向量必須是由一個方陣求解得到。但也不是所有的方陣都有特征向量。對於給定一個 n * n 的矩陣，並且該矩陣存在特征向量，那么這個矩陣有n個特征向量。類似的，對於給定一個 3 * 3 的矩陣（假設存在特征向量），那么這個矩陣有3個特征向量。

   特征向量的另一個性質是，即使變換矩陣A乘以原始向量V =  之前，先將V擴大2倍，再將矩陣與擴大后的向量 V' 相乘，得到的結果仍然是變換矩陣的特征向量，且同樣擴大了2倍，如下圖所示。

向量與矩陣相乘前，乘以倍數2

變換矩陣乘以擴大了2倍的原始向量

    這是因為，如果你對原始向量擴大n倍，意味着將向量變得更長但不改變向量的方向。最后，由變換矩陣得到的所有特征向量都是正交的，不管數據集有多少維，兩兩都相互垂直。這部分知識非常重要，因為這意味着你可以通過這些正交的特征向量來數據表達，而不用根據 x-y 坐標系。我們將在下文的PCA章節進行討論。（參考： 基與正交基,特征值與特征向量）

    另一個讀者們應該知道的重點是，“數學家們在求解特征向量時，他們更願意去找到這些模正好為1的特征向量”。這是因為，如讀者們所知，向量的大小並不影響這個向量是否為特征向量，但是方向就有影響。因此，為了標准化特征向量，當我們求得一個特征向量時，通常將它的模縮放為1，使得所有特征向量都有相同的模（模為1）。如下，通過一個例子加以驗證。（此部分更詳細的內容可參考： PCA數學原理）

向量   是一個特征向量，這個向量的模為： ，然后我們將原始向量   除以模 ，使得特征向量的模變為1，如下：

    那么，我們如何着手去找到這些神秘的特征向量呢？不幸的是，能夠簡單的進行求解的條件是你處理的對象是一個很小的矩陣，就像一個規模不大於 3*3 的方陣。對於求解規模較大方陣的特征向量，通常的方法是通過一些復雜的迭代方法進行求解（這些方法超出了本教程的范圍，作者沒做進一步的介紹）。如果讀者們需要通過程序求解方陣的特征向量，可以找一個數學庫（例如，Matlab，你只需要調用相關的函數即可，具體的實現過程這個庫已經幫你完成了）。此外，有一個叫做“newmat”的數學工具包僅供讀者參考： http://webnz.com/robert/。

    想要進一步了解“關於如何求解一般情況下的特征向量，正交性”的讀者們，可以參考由Howard Anton編著，John Wiley & Sons股份有限公司出版的教科書“基礎線性代數 第五版”（ISBN 0-471-85223-6）。

2.2.2 特征值

    特征值跟特征向量緊密相關，實際上，讓我們回顧一下“2.2.1 特征向量”這節的2個例子，我們可以注意到這2個例子中，“相同的方陣A，分別乘以原始向量   以及擴大2倍后的向量  ，最后得到的特征向量前面的系數是一樣的”。在這2個例子中，得到的系數都是4，4就是這個這個特征向量對應的特征值。在方陣A乘以向量V之前，無論我們先對向量V乘上一個多大的系數（非零：線性代數中規定特征向量不可以為零向量）得到的V'，最后得到的都會是4乘以 V' 。（ 特征值的定義：設A為n階矩陣，若存在常數λ及n維非零向量x，使得Ax=λx，則稱λ是矩陣A的特征值，x是A屬於特征值λ的特征向量）

    讀者們可以發現特征向量和特征值是成對出現的。當你使用一個理想的程序庫去計算特征向量時，同時你也會得到相應的特征值。

   練習：

    對於下列矩陣：

    分別對下列的5個向量，判斷是否為上述方陣的特征向量，如果是，求出該特征向量向量對應的特征值。

，  ，  ，  ， 

    【補充】可以通過Matlab提供的“eig()”函數求解特征向量及其對應的特征值，如下所示：

close all;
clc;
clear;

A = [3 0 1; -4 1 2; -6 0 -2];
[d, v] = eig(A)

A * d - v * d

    如上實驗結果所示，方陣A求得2個特征向量 [0; 1; 0] 和[0; -1; 0] 對應的特征值都是1。而上述例子中所給的5個向量中，只有[0;1; 0]是方陣A的特征向量。值得注意的是：對於上述例子的求解方法，應該根據特征值和特征向量的定義進行求解，即根據是否滿足“A*Vector = egienvalue*Vector”的形式，如果滿足，則向量Vector是方陣A的特征向量，反之不是。

第3章：主成分分析（PCA）

    最后，我們來到了主成分分析（PCA）這節。主成分分析是個什么東西呢？它是一種“識別數據中的模式，通過強調數據的異同（similarities and differences）的方式來表達數據”的方法。由於很難發現高維數據中的模式，對於圖形表示法是不可用的，而PCA正是一種強大的數據分析工具。

    PCA的另一個主要優勢是，一旦你發現數據中的這些模式，你就可以通過減少維數來壓縮數據（保留主成分的數據），而且不會損失太多的信息。這項技術（優勢）運用在圖像壓縮上，后面章節會具體介紹。

    本章將為讀者們介紹“在一組數據集上實現主成分分析”所需的相關步驟。本文並不會詳細介紹PCA這項技術為什么能起作用（下文中，作者沒有具體給出計算過程，是直接給出結果），但是作者會嘗試着提供“每個要點（步驟）會發生什么（有什么作用）”的解釋，因此當讀者們嘗試着使用到PCA這項技術時，能夠做出明確的決定。

3.1 理論方法（PCA的六大步驟）

步驟一：獲取數據集

    為了簡單起見，本文將使用作者自己構造的數據集。這是一組2維的數據集，選擇這樣一組數據集是因為能夠對其作圖，更直觀的顯示PCA分析過程中的每一個步驟的效果。如下左圖所示是本章的原始數據集Data和減去均值后的數據集DataAdjust，右圖是對DataAdjust進行作圖的效果圖。

         

步驟二：減去均值

    為了讓PCA正常工作，必須將原始數據集中的每一維的所有元素值減去該維的均值（數據的標准化）。（此處的數據集是2維數據，即x-y）也就是說，所有的 值都減去 （ 就是x維上所有元素值的均值），所有的   值都減去  。這個過程使得我們得到一組均值為0的數據集。

步驟三：計算協方差矩陣

    這個步驟正如2.1.4節討論過的計算協方差矩陣。因為這是一組2維數據集，所以將會得到一個 2 X 2 的協方差矩陣。在這里也是一樣（無意外），因此作者直接給出結果，如下：

    根據上述結果，由於協方差矩陣   非對角線上的元素值是正數，我們可以推出變量 x 和 y 是正相關（同增同減）。

步驟四：計算協方差矩陣的特征向量及其對應的特征值

    因為協方差矩陣是一個方陣，所以我們可以計算協方差矩陣的特征向量及其對應的特征值，這一步驟是非常重要的。特征向量及其對應的特征值告訴我們關於數據有用的信息，待會兒將展示給讀者們看。在此期間，計算所得的特征向量和特征值如下所示：

    ，     

注意：原文中給的特征向量的符號跟Matlab求得的不一樣（這里給出Matlab求得的解），基於Matlab求解特征值和特征向量的代碼如下所示：

cov = [0.616555556, 0.61544444; 0.616555556, 0.716555556];  % 協方差矩陣
[eigenvectors, eigenvalues] = eig(cov);  % 求解特征向量及其對應的特征值

值得注意的是，此處求解得到的這些特征向量都是單位特征向量，即 ，這兩個特征向量的模都為1。這一點對於PCA是非常重要的，幸運的是，大多數數學工具包求解特征向量時，得到的結果都是單位特征向量。
    那么，這又是什么意思呢？如下圖所示為標准化后的數據集（各維變量減去對應的均值），讀者們可以發現這些數據中存在一個明顯的模式（quite a strong pattern）。正如協方差矩陣所預期的一樣，x和y兩組變量的值呈正相關（同增同減）。在圖的上方，作者也畫出了對應的兩個特征向量所表示的直線方程，它們以虛線的形式出現於坐標系的對角線。正如前面提到的一樣，這一組特征向量相互垂直（正交），分別是圖中斜率大於0以及斜率小於0的這兩條虛線。但更重要的是，特征向量為我們提供了數據模式（此處的模式就是直線方程）的相關信息。首先讓我們先關注這條經過數據集（標准化后）中心的特征向量對應的直線（斜率小於0那條直線），看起來像畫了一條最擬合的直線嗎（很明顯不是）？這個特征向量告訴我們這兩組數據集（應該指的是這兩維的數據，x-y）是如何關聯到這條直線。（這就是我們要找的直線，數據模式）而第二個特征向量為了我們提供了另一條直線（另一種數據模式），斜率大於0的那條線，坐標上所有的點沿着這條主線分布，但都以一定程度的距離偏離在這條主線的兩側。

    所以，通過求解協方差矩陣的特征向量，我們已經能夠提取出描述數據集的直線。現在，我們剩下的步驟包括，對數據進行轉化，使其能夠表達在特征向量對應的直線上。

步驟五：選擇主成分，形成特征向量

    接下來就是“數據壓縮和降維”的到來。如果讀者們仔細觀察前幾章講到的特征向量和特征值，注意到求解得到的特征值是不一樣的。實際上，最大的特征值對應的特征向量就是數據集的主成分。在本文的例子里，最大特征值對應的特征向量就是x-y坐標系中，斜率大於0的那條直線（此處論文好像有點問題，讀者們可查看原文進一步驗證）。這就是表達數據集維度之間的關系最有意義，最具代表性的模式。

    通常來說，一旦從協方差矩陣中求解得到特征向量，下一步就是從大到小對特征值進行排序（此處要對應好特征向量和特征值），也就是對數據成分主次的排序（特征值越大，數據在該維就越有意義，信息量越大）。現在，如果讀者們有需求的話，就可以忽略那些成分較小的維度。雖然你丟失了一些數據信息，但是忽略的這些維度的特征值很小，丟失的信息很少。如果你忽略一些成分（維度），最終得到的數據集的維數將比原始數據集的維數少（這就是數據降維）。（To be precise）確切地講就是，如果你處理的原始數據集是n維，你將計算得到n個特征向量和特征值（有對應關系），然后你選擇特征值前p（p < n）大的特征向量，忽略掉其他維度的數據，最后得到的數據集就是一組p維的數據集。

    至此，讀者們現在需要做的事情就是構造一個特征向量（fearue vector），其實就是由向量構成的一個矩陣。  就是取前k大的特征值對應的eigenvectors，然后根據如下形式，將這k個向量以列的形式構成一個矩陣。

    回到本文例子的數據集，通過求解，我們有2個特征向量   ，  ，意味着我們有2種選擇。如下所示，我們可以同時選擇這2個特征向量作為 ：

  

    我們也可以選擇忽略特征值較小的那維數據，僅僅保留特征值較大的那維數據，如下所示：

    作者將在下一節展示以上兩種結果。

（注：以上的 FeatureVector ，下文運用的時候，需要對其先進行轉置）

步驟六：產生新的數據集（降維后）

    這是PCA的最后一個步驟，也是最簡單的。一旦我們選好了主成分（eigenvectors），並構造了特征向量（ ）后。我們首先對  進行轉置（得到1行2列或者2行2列的 ），然后左乘以原始數據（已標准化且轉置后，即2行n列），如下所示：

其中，  表示提取前k維主成分后構成的   並進行轉置后的矩陣（  是按列進行存儲，經過轉置后的  變成按行進行存儲）；  表示原始數據減去均值后再進行轉置得到的結果，  矩陣的每一列表示一組具體的數據項，如 ，每一行表示獨立的一個維（如，x維或者y維）。讀者可能會因為這個步驟中出現的突然“轉置”感到困惑，作者為此感到很抱歉，但如果讀者們先對 和 進行轉置，這里的給出的公式將更加簡單。而不是在 和 的右上方表上一個轉置符號T（其實，翻譯這里的時候，我不是很清楚這句話的用意）；  表示經過PCA算法后最終得到的數據集，其中每一列表示數據項 ，每一行表示具體的維（x維或者y維）。

    行文至此，（降維后的   ）我們將獲得什么呢？僅僅根據我們選擇的前k維特征向量，我們就可以得到原始數據。案例中的數據集有2維，包括x維和y維。讀者們可以選擇自己喜歡的兩個維度來表示這些數據。如果這些維度是相互垂直（perpendicular）的，那么以這種方式表達數據會更有效。這就是為什么特征向量需要兩兩互相垂直會這么重要。我們已經改變了數據的表示形式（原來是基於x-y坐標系），現在是表示在2個特征向量上。當新的數據集進了降維的情況下，我們已經是忽略掉一些特征向量（ ），只保留了我們選擇的前k個特征向量（ ）。

    為了在數據上展示這些，作者已經對每一種可能的特征向量（ ）完成了最后的變換（  左乘 ）。同時作者也對每一種情況的結果進行了轉置，將數據轉換回原始漂亮的表格格式（就是將更漂亮，更直觀的表示格式）。作者還對最后得到的新數據進行作圖，體現這些數據是怎樣關聯到成分（提取的前k個主成分）。

• 第一種情況（保留了所有的特征向量）：如下左圖表示經過 變換后的 ，右圖表示在新的坐標系（）下繪制 這些數據坐標點，可以理解為原始數據和坐標系經過旋轉一定的角度得到的。不難理解，這種情況下，沒有損失任何數據。

• 第二種情況（只選擇最大特征值對應的特征向量）：如下圖所示為降維后的數據集（從2維降到1維，和預期的一樣）。讀者們如果將這組數據與第一種情況的數據作比較，你會發現這組數據就是第一種情況那組數據的第一列，即對應 x 維的數據。

    所以，如果你想對這組數據作圖，如下圖所示，得到的就是一維空間，而且圖上 y=0 直線附近的這些坐標點恰好是 x 維空間上的位置（此時得到的坐標系，相當於y=0，主成分都在x軸上的數據）。我們已經拋掉另一維空間信息（y維），也就是另一個特征向量（值較小的那個特征向量被拋棄了）。

    那么，到這里我們完成了什么呢？我們基本上已經對數據進行了變換，使得數據表達成它們之間的模式，本例中這些模式是“描述這些數據之間的關系的最擬合的直線”。這是非常有用的，因為我們現在已經將數據“根據這些直線（也就是特征向量eigenvector）對原始數據貢獻信息的程度”進行組合分類。最初，我們有 x 和 y 坐標軸（完整的數據信息），但每個數據點的 x 值和 y 值實際上並不能確切告訴我們這個數據跟其他數據之間的關系（y值提供的信息量很少）。現在（降維后），數據點的坐標值（只有 x 值）告訴我們這個數據點適合於哪條趨勢線（trend lines）。本例中，兩個eigenvector都使用（轉換）的情況下，我們僅僅簡單的改變了數據，用求解得到的兩個特征向量（eigenvectors）代替原本的 x-y 坐標系而已。但是，將特征值較小（y維，貢獻小）的那一維去掉，只保留與特征值較大的那一維（x維）相關的數據的這種情況，就達到很好的降維效果，而且在丟失少量數據信息的情況下，用一維就能夠很好的表示原始數據。

3.1.1 恢復舊數據

    如果讀者們使用 PCA 進行數據壓縮，一定都會想要“恢復原始數據”，這個問題想必受到很多讀者的關注，畢竟數據經過PCA變換后，仍是丟失了一些信息（下文有提供相關的案例）。這部分內容摘抄於 http://www.vision.auc.dk/ sig/Teaching/Flerdim/Current/hotelling/hotelling.html。（這是原文提供的網址，但是好像進不去了- -）

    那么，我們該如何恢復原始數據呢？在這之前，讀者們應該知道只有在“將所有特征向量（eigenvectors）共同構成最終的變換矩陣”這種情況下，我們才能精確的恢復原始數據。如果最終的變換矩陣已經是經過降維（丟掉一些特征向量，比如上個例子中的 y 軸被丟棄），那么原始數據已經是丟失掉一些信息了，盡管是少量的信息。

    回顧一下之前的“最終變換矩陣”，如下：

    為了恢復原始數據，我們對上述公式做一下變形，得到如下公式：

其中，  是   的逆（inverse）。當我們將所有的eigenvectors作為feature vector時（因為eigenvector和feature vector都翻譯成“特征向量”，所以此處，我直接使用英文表示，加以區分），你會發現，原來矩陣 feature vector 的逆恰好等於feature vector的轉置（當且僅當滿足“feature vector矩陣的元素是數據集對應的單位特征向量，eigenvectors”這個條件的時候，成立。證明過程可以參考， 這里）。這使得恢復數據變得更加容易，因為上述的表達式得到了改善，如下所示：

但是，為了恢復到最原始數據，我們還需要對每個維度的元素都加上相應的均值（因為，我們一開始的時候就對原始數據做了標准化，也就是各維度的元素都減去了均值）。所以，如下所示，是更完整的表達式：

這個公式對“並沒有把所有的 eigenvectors 作為feature vector”這種情況也同樣適用。所以即使當你丟棄一些eigenvector時，上述等式仍然可以得到正確的轉換。

    作者在本文中並不打算使用完整的feature vector（使用所有的eigenvector）演示數據恢復，因為這樣得到的結果恰好就是我們一開始處理的數據集（最原始的數據）。然而，作者將做如下實驗：“僅保留特征值較大的維度空間的情況，展示數據是如何丟失的”，如下圖所示（降維后的數據），將圖中的數據與最原始的數據作比較，你會發現只有沿着主成分eigenvector方向上的變量（x維）被保留了下來，沿着另一個eigenvector方向上的變量（y維）已經不見了。

    練習：

• 求解協方差矩陣得到的特征向量能夠帶給我們什么（有什么作用）？
• 在PCA流程中，我們能夠在“哪個步驟”決定壓縮數據？會有什么影響？
• 找一個關於“PCA運用到面部識別”的案例並用圖示法表示主要的eigenvectors（前k大的特征值對應的eigenvectors），可以通過“Eigenfaces”關鍵字進行搜索。

第4章：計算機視覺領域的應用（Application to Computer Vision）

    本章將對“PCA在計算機視覺中的應用”進行概述。首先介紹圖像通常是怎樣表示的，然后介紹PCA能夠幫助我們對這些圖像做什么樣的處理。在本章中，關於“面部識別（facial recognition）”的相關信息是來自於“Face Recognition: Eigenface, Elastic Matching, and Neural Nets”, Jun Zhang et al. Proceedings of the IEEE, Vol. 85, No. 9, September 1997。 表征信息來自於“Digital Image Processing”, Rafael C. Gonzalez and Paul Wintz, Addison-Wesley Publishing Company, 1987（這同樣是關於“一般情況下， K-L變換更進一步的知識”一篇的很好的參考文獻）。圖像壓縮的相關資料來自於“ http://www.vision.auc.dk/ sig/Teaching/Flerdim/Current/hotelling/hotelling.html”（又進不去0- -，原文是2002年的，可能有些網址已經發生了改變，后續有找到的話，再更新分享給大家），這個參考網址也為我們提供一些關於“使用不同數量的eigenvector的圖像重建”的案例。

4.1 （圖像）表示方法（Representation）

    在計算機視覺領域，當我們使用矩陣方法時，一定要考慮圖像的表示形式。一個 N*N 的方陣可以被表示成一個  維的向量，如下所示：

其中，圖像中的像素以行為單位，一行接着一行被放置，形成一個一維的圖像。例如，前N個元素（  - ）將被作為圖像的第一行，接下來的N個元素（  - ）作為第二行，以此類推。向量中的值反應圖像的強度信息，或許就是一個簡單的灰度值（對於灰度圖像，就是對應灰度值）。

補充，如下矩陣所示，表示一個圖像進行數字化的矩陣表示形式：

4.2 基於PCA尋找模式（模型）（PCA to find patterns）

    假設我們有20張圖像，每張圖像由N個像素的高和N個像素的寬組成（N*N的矩陣）。對於每張圖像，我們可以使用上一節的方法將其表示為一個圖像向量。然后，我們可以將所有的圖像（現在一張圖像對應一個向量）放到一個大矩陣里面，形如：

    以此作為我們使用PCA算法的第一步（現在處理的對象是一個由20張圖像構成的大矩陣）。一旦使用PCA，我們要做的就是求解協方差矩陣得到特征向量（eigenvectors）。這（PCA）為什么有用呢？假設我們想實現面部識別，原始的數據集是人臉圖像。問題是，給定一張新的人臉圖像，識別出這張新圖像對應原始數據圖像中的哪一類人臉圖像，也就是根據人臉圖像信息進行分類（注意，這張新圖像並不是來自於我們一開始所給的那20張圖像）？這個問題在計算機視覺的解法是，基於PCA分析得到的新坐標系下，測量新的人臉圖像與已知的20張人臉圖像之間的區別，而不是在原來的坐標系。

    事實證明經過PCA算法得到的新的坐標系更有利於識別人臉，這是因為PCA算法告訴我們原始圖像（數據）之間的差異（differences and similarities）。PCA算法確定了數據中的統計模型。

    因為所有的向量都是  維，所以我們將得到 個特征向量（eigenvector）。實際上，我們也可以丟棄一些意義不大的eigenvectors（只保留特征值前k大對應的eigenvectors），識別的效果同樣不錯。

4.3 基於PCA的圖像壓縮（PCA for image compression）

    使用PCA算法進行圖像壓縮又稱為Hotelling 變換或者K-L變換。假如我們有20張圖像，每張 個像素。我們可以構造 個向量，每個向量20維，每一維對應這20張圖像中同一個像素的強度值，下文我（博主）將補充說明。這一點與上一個例子的大不同，上一個例子是構成的向量vector中的每個元素都是對應不同的像素，而現在這個例子構成的每個向量（ 個）的元素對應20張圖像相同的一個像素值。

---------------

補充說明： 個20維的向量

• 第1個向量形如：
• 第2個向量形如：
• 第個向量形如：

其中，每個向量中有20個像素值，比如第1個向量  的20個元素值分別對應從第一張圖像到第20張圖像的第一個像素值 ，以此類推。

---------------

    現在開始對這組數據使用PCA算法。我們將得到20個特征向量（eigenvectors），因為每個vector都是20維。為了壓縮數據，我們可以使用特征值前15大對應的15個eigenvector作為變換矩陣。最后得到的數據集只剩下15維，減少了1/4的空間。然而，但需要重構原始數據集時，得到的圖像將丟失一些信息。PCA是一項有損信息的壓縮技術，因為解壓后的圖像並不跟原始數據完全一樣，通常更糟（差異更大）。

5. 博主補充：PCA實例

• 原始數據集：，使用PCA算法將這組二維數據降維成一維數據。
• 因為這個矩陣的每行已經是零均值，這里我們直接求協方差矩陣，過程如下：

• 然后求其特征值和特征向量，具體求解方法不再詳述，可以參考相關資料。求解后特征值如下所示：

• 其對應的特征向量分別是：

•  其中對應的特征向量分別是一個通解， 和  可取任意實數。那么標准化后的特征向量（是一個單位特征向量）為 ：

• 因此我們的矩陣P是：

• 可以驗證協方差矩陣C的對角化：

• 最后我們用P的第一行乘以數據矩陣，就得到了降維后的表示：

• 降維投影結果如下圖：

注：本例（來自於文章“ PCA數學原理”）在求解協方差矩陣的過程中，協方差的分母是用n表示。但是我們在本教程中提到，處理樣本而非總體的時候，我們更傾向於選擇 (n-1) ，即 。在此，我們保留以上例子，但為了讓例子更貼近本教程，博主使用Matlab對該例子進行求解，協方差的分母使用的是 (n-1)，代碼和實驗結果如下所示：

• 數據集還是上個例子的數據集：
• 求解協方差矩陣，代碼如下：

close all;
clc;
clear;

%% PCA案例
X = [-1,-1,0,2,0];
Y = [-2,0,0,1,1];
mean_X = mean(X);
mean_Y = mean(Y);
n = 5;
sum_xy = 0;
sum_yx = 0;
sum_xx = 0;
sum_yy = 0;

for i = 1:5
    sum_xy = sum_xy + (X(i) - mean_X) * (Y(i) - mean_Y);
    sum_yx = sum_yx + (Y(i) - mean_Y) * (X(i) - mean_X);
    sum_xx = sum_xx + (X(i) - mean_X) * (X(i) - mean_X);
    sum_yy = sum_yy + (Y(i) - mean_Y) * (Y(i) - mean_Y);
end

cov_xx = sum_xx / (n-1)
cov_xy = sum_xy / (n-1)
cov_yy = sum_yy / (n-1)
cov_yx = sum_yx / (n-1)

    -- 也可以使用Matlab提供的cov函數直接得到協方差矩陣：

 

       -- 協方差矩陣：

• 根據協方差矩陣求解特征向量和特征值：

cov_test = cov(X,Y);    % 協方差矩陣
[eigenvectors eigenvalues] = eig(cov_test);  %特征向量和特征值

       -- 特征向量（通過Matlab的eig求得的特征向量已經是標准化的）：

       -- 特征值：0.5 和 2.5

• 選擇最大特征值（此處為2.5）對應的特征向量作為feature vector：
• 根據公式 ，求解降維后的新數據集：（從2維降到1維）

X = [-1,-1,0,2,0];
Y = [-2,0,0,1,1];

% 特征向量
featurevector = [0.7071;0.7071];     % 直接提取最大特征值對應的特征向量
RowFeatureVector = featurevector';   % 轉置
RowAdjustData = [X;Y];
FinalData = RowFeatureVector * RowAdjustData

附錄A

實現代碼：

    本文代碼是基於Scilab（一款開源的軟件）實現的。與Matlab類似， SCILAB也是一種科學工程計算軟件，其數據類型豐富，可以很方便地實現各種矩陣運算與圖形顯示，能應用於科學計算、數學建模、信號處理、決策優化、線性/非線性控制等各個方面。作者使用以下代碼生成了本文中的所有案例。代碼中除了第一部分的宏（macro）以外，其余部分都是作者自己完成的。

    注：這部分代碼，本人尚未進一步驗證，后續驗證了再更新。

%  This macro taken from
%  http://www.cs.montana.edu/˜harkin/courses/cs530/scilab/macros/cov.sci
%  No alterations made
%  Return the covariance matrix of the data in x, where each column of x
%  is one dimension of an n-dimensional data set. That is, x has x columns
%  and m rows, and each row is one sample.
% 
%  For example, if x is three dimensional and there are 4 samples.
%  x = [1 2 3;4 5 6;7 8 9;10 11 12]
%  c = cov (x)
function [c]=cov (x)
    % Get the size of the array
    sizex=size(x);
    % Get the mean of each column
    meanx = mean (x, 'r');
    % For each pair of variables, x1, x2, calculate
    % sum ((x1 - meanx1)(x2-meanx2))/(m-1)
    for var = 1:sizex(2),
        x1 = x(:,var);
        mx1 = meanx (var);
        for ct = var:sizex (2),
            x2 = x(:,ct);
            mx2 = meanx (ct);
            v = ((x1 - mx1)' * (x2 - mx2))/(sizex(1) - 1);
            cv(var,ct) = v;
            cv(ct,var) = v;
            % do the lower part of c also.
        end
    end
    c=cv;
    
% This a simple wrapper function to get just the eigenvectors
% since the system call returns 3 matrices
function [x]=justeigs (x)
% This just returns the eigenvectors of the matrix
    [a, eig, b] = bdiag(x);
    x= eig;
%  this function makes the transformation to the eigenspace for PCA
%  parameters:
%  adjusteddata = mean-adjusted data set
%  eigenvectors = SORTED eigenvectors (by eigenvalue)
%  dimensions = how many eigenvectors you wish to keep
% 
%  The first two parameters can come from the result of calling
%  PCAprepare on your data.
%  The last is up to you.
function [finaldata] = PCAtransform(adjusteddata,eigenvectors,dimensions)
    finaleigs = eigenvectors(:,1:dimensions);
    prefinaldata = finaleigs'*adjusteddata';
    finaldata = prefinaldata';
%  This function does the preparation for PCA analysis
%  It adjusts the data to subtract the mean, finds the covariance matrix,
%  and finds normal eigenvectors of that covariance matrix.
%  It returns 4 matrices
%  meanadjust = the mean-adjust data set
%  covmat = the covariance matrix of the data
%  eigvalues = the eigenvalues of the covariance matrix, IN SORTED ORDER
%  normaleigs = the normalised eigenvectors of the covariance matrix,
%  IN SORTED ORDER WITH RESPECT TO
%  THEIR EIGENVALUES, for selection for the feature vector.
% 
%  NOTE: This function cannot handle data sets that have any eigenvalues
%  equal to zero. It’s got something to do with the way that scilab treats
%  the empty matrix and zeros.

function [meanadjusted,covmat,sorteigvalues,sortnormaleigs] = PCAprepare (data)
% Calculates the mean adjusted matrix, only for 2 dimensional data
means = mean(data,'r');
meanadjusted = meanadjust(data);
covmat = cov(meanadjusted);
eigvalues = spec(covmat);
normaleigs = justeigs(covmat);
sorteigvalues = sorteigvectors(eigvalues', eigvalues');
sortnormaleigs = sorteigvectors(eigvalues', normaleigs);


%  This removes a specified column from a matrix
%  A = the matrix
%  n = the column number you wish to remove
function [columnremoved] = removecolumn(A,n)
inputsize = size(A);
numcols = inputsize(2);
temp = A(:,1:(n-1));
for var = 1:(numcols - n)
temp(:,(n+var)-1) = A(:,(n+var));
end,
columnremoved = temp;


%  This finds the column number that has the
%  highest value in it’s first row.
function [column] = highestvalcolumn(A)
inputsize = size(A);
numcols = inputsize(2);
maxval = A(1,1);
maxcol = 1;
for var = 2:numcols
if A(1,var) > maxval
maxval = A(1,var);
maxcol = var;
end,
end,
column = maxcol

%  This sorts a matrix of vectors, based on the values of
%  another matrix
% 
%  values = the list of eigenvalues (1 per column)
%  vectors = The list of eigenvectors (1 per column)
% 
%  NOTE: The values should correspond to the vectors
%  so that the value in column x corresponds to the vector
%  in column x.
function [sortedvecs] = sorteigvectors(values,vectors)
    inputsize = size(values);
    numcols = inputsize(2);
    highcol = highestvalcolumn(values);
    sorted = vectors(:,highcol);
    remainvec = removecolumn(vectors,highcol);
    remainval = removecolumn(values,highcol);
    for var = 2:numcols
        highcol = highestvalcolumn(remainval);
        sorted(:,var) = remainvec(:,highcol);
        remainvec = removecolumn(remainvec,highcol);
        remainval = removecolumn(remainval,highcol);
    end
    sortedvecs = sorted;


%  This takes a set of data, and subtracts
%  the column mean from each column.
function [meanadjusted] = meanadjust(Data)
    inputsize = size(Data);
    numcols = inputsize(2);
    means = mean(Data,'r');
    tmpmeanadjusted = Data(:,1) - means(:,1);
    for var = 2:numcols
        tmpmeanadjusted(:,var) = Data(:,var) - means(:,var);
    end
    meanadjusted = tmpmeanadjusted

附錄B：本文相關的Matlab作圖代碼

1. 【協方差部分】H和M的關系圖

%% ctrl + R 注釋
%% ctrl + T 取消注釋

close all;
clc;
clear;
fontsize = 11;
figure;
% x = 0:0.00025:1;
% y = x.*sin(10*3.14.*x)+2;
% plot(x,y)

x = [9 15 25 14 10 18 0 16 5 19 16 20];
y = [39 56 93 61 50 75 32 85 42 70 66 80];
plot(x, y, 'o');
xlabel({'H / 小時'}, 'FontSize',fontsize);
ylabel({'M / 分數'}, 'FontSize',fontsize);

h = legend(['H(學習時間)和M(獲得分數)的關系', sprintf('\n')], 'location', 'best');
grid on;

set(h, 'Box', 'off');

參考文獻

1. PCA數學原理：http://www.360doc.com/content/13/1124/02/9482_331688889.shtml
2. 杜洪，夏欣，琚生根，王能。基於PCA圖像壓縮算法研究與實現：http://www.cnki.com.cn/Article/CJFDTotal-SCDX201405009.htm
3. 本文對應的原文：【PDF】http://pan.baidu.com/s/1nvNPJ4H   |   【文章出處】https://arxiv.org/abs/1404.1100
4. 博主的CSDN主頁：http://blog.csdn.net/houchaoqun_xmu
5. 【Latex】制作演示文檔或者課程報告-常見用法（一）
6. 【Latex】制作演示文檔或者課程報告 - 制作演示文檔（二）
7. 【Latex】制作演示文檔或者課程報告 - 制作課程報告（三）
8. TSP_旅行商問題 - 蠻力法DFS（一）
9. TSP_旅行商問題 - 模擬退火算法（三）
10. TSP_旅行商問題 - 遺傳算法（四）