เขียน PCA บน MATLAB

บล็อกนี้จะมานำเสนอการทำ Principal Components Analysis (PCA) บน MATLAB กันเองนะครับ ผมไม่อธิบายนะครับว่า PCA คืออะไร ใครอยากรู้มากขึ้นเกี่ยวกับ PCA แนะนำให้ไปอ่านที่วิกิ หรือที่บล็อกของ อ.ปืน นะครับ เราลองมาดูกันว่าเราสามารถเขียน PCA ได้เองอย่างไรบน MATLAB

เริ่มต้นเรามีข้อมูลอยู่ชุดหนึ่ง มี 10 ค่า และมี 3 มิติ (dimension) ในที่นี้ผมมั่วตัวเลขขึ้นมาเองนะครับ เราก็หาจำนวนแถวและคอลัมภ์ของข้อมูลไว้ด้วย เพราะต่อไปต้องใช้

x = [ 3 2 3; 
      6 5 6;
      9 8 10;
      1 8 9;
      4 5 2;
      1 9 3;
      2 10 9;
      8 1 1;
      3 4 6;
      7 2 5 ];
[Rows, Columns] = size(x);

เสร็จแล้วเราก็หาค่าเฉลี่ย (mean) และค่าความแปรปรวนร่วมเกี่ยว (covariance) ของข้อมูลชุดนี้โดย

m = mean(x);
y = x - ones(size(x, 1), 1) * m;
c = cov(y);

หลังจากนี้ MATLAB ก็จะมีฟังก์ชั่นแจ่มๆ ให้ใช้ เพื่อหา eigenvector และ eigenvalue ดังนี้

[V, D] = eig(c);

โดย V คือ eigenvector และ D คือ eigenvalue จะง่ายไปไหนครับ? เสร็จแล้วเหรอ? ยังครับ เรายังต้องเรียงค่าของ eigenvalue จากมากไปน้อยด้วย ทำได้ดังนี้

[D, idx] = sort(diag(D));
D = D(end:-1:1)';

ตัวแปร idx จะเก็บค่า index ของ eigenvalue นั้นไว้ แล้วเราก็มาเรียง eigenvector ตาม eigenvalue ได้ดังนี้

V = V(:, idx(end:-1:1));

เท่านี้เราก็จะได้สิ่งที่เราต้องการ นั่นก็คือ eigenvector กับ eigenvalue ที่เรียงลำดับค่าจากมากไปน้อยแล้ว เราสามารถนำค่าพวกนี้ไป

1. ลดจำนวนมิติ (dimension reduction) หรือ
2. จำแนกพวก outliers ได้

ถ้าต้องการทำ dimension reduction ก็ต้องเลือกจำนวน dimension ที่เราต้องการมา ในที่นี้ผมเลือกมา 2 จาก 3 นะครับ แล้วเราก็เลือก eigenvector แค่ 2 อันแรกมา

OutputSize = 2;
V2d = V(:, 1:OutputSize);

จากนั้นก็ project ข้อมูลเดิมลงไปบน subspace ใหม่ (การคูณกันของเมตริกซ์) ได้เลย

prefinal = V2d'*y';
final = prefinal'

เสร็จสิ้นการทำ dimension reduction

ส่วนการจำแนก outliers ทำได้ดังนี้ สมมุติว่าผมมีข้อมูลที่อยากตรวจสอบ

new_data = [3 2 5];

ผมก็เอาข้อมูลนั้นมาลบด้วยค่าเฉลี่ยที่ผมได้หาว่าตอนแรก แล้วก็ project ลงไปบน subspace

new_data_subtracted = new_data - ones(size(new_data, 1), 1) * m;
output = V2d' * new_data_subtracted';

เราก็จำแนกข้อมูลนี้โดยดูค่า Mahalanobis distance ถ้าค่านี้มากกว่าค่า threshold ที่เราตั้งไว้ ข้อมูลนี้ถือว่าเป็น outliers ครับ 🙂

distance = sum( ( output'./sqrt( D(1:OutputSize) ) ).^2 );

เนื่องจาก eigenvector (แกนของ PCA) นั้นตั้งฉากกัน Mahalanobis distance ก็คือค่า Euclidean distance นั่นเอง หลังจากที่ normalize ด้วยค่ารากที่สองของ eigenvalue (เนื่องจาก eigenvalue แทนค่า variance แต่เราต้องหารด้วย standard deviation)

ใครอยากลองเล่น ก็โหลดโค้ดกันไปได้เลยครับ pca.zip
โค้ดข้างต้นนี้ ผมอ้างอิงจาก บล็อกของ James Rossiter ครับ

จริงๆ แล้วมีอีกวิธีหนึ่งที่สามารถใช้หา eigenvector และ eigenvalue ได้เช่นกัน นั่นก็คือใช้ Singular Value Decomposition (SVD) นั่นเอง ซึ่ง MATLAB ก็มีฟังก์ชั่นให้ใช้ด้วย โดยเราแค่หา SVD จาก covariance ของข้อมูลเราดังนี้

[U,S,V] = svd(c);

โดย eigenvalue ก็คือค่าบนแกน diagonal ของเมตริกซ์ S และ eigenvector ก็คือคอลัมภ์ของเมตริกซ์ V นั่นเอง! คอลัมภ์ของ U ก็เป็น eigenvector เหมือนกันแต่เป็นของคนละเมตริกซ์ รายละเอียดสามารถหาอ่านได้ตามเว็บครับ

การใช้ svd() จะดีกว่า eig() ตรงที่ eig() ต้องใส่ input เป็นเมตริกซ์จัตุรัส แต่ svd() นั้นไม่จำเป็น

Author: zkan

Soon to be a newbie data scientist. I ♥ machine learning, computer vision, robotics, image processing, data visualization, and data analytics. View all posts by zkan