φυβλαςのβλογ
บล็อกของ phyblas



การวิเคราะห์องค์ประกอบหลัก (PCA) ด้วยการแยกค่าเอกฐาน (SVD)
เขียนเมื่อ 2019/09/16 20:50
ในบทความก่อนหน้านี้ได้เขียนถึงการวิเคราะห์องค์ประกอบหลัก (主成分分析, principle component Analysis, PCA) ไป https://phyblas.hinaboshi.com/20180727

โดยได้แสดงถึงว่าการวิเคราะห์องค์ประกอบหลักนั้นสามารถทำได้โดยหาเวกเตอร์ลักษณะเฉพาะ (本征向量, eigenvector) ของเมทริกซ์ความแปรปรวนร่วมเกี่ยว (协方差矩阵, covariance matrix)

แต่ว่าในทางปฏิบัติแล้ว การคำนวณเมทริกซ์ความแปรปรวนร่วมเกี่ยวเพื่อมาหากเวกเตอร์ลักษณะเฉพาะนั้นเป็นสิ่งที่กินเวลามาก โดยเฉพาะถ้าข้อมูลมีขนาดใหญ่

ในทางปฏิบัติแล้ว วิธีที่มักใช้แทนก็คือ การแยกค่าเอกฐาน (奇异值分解, singular value decomposition) หรือนิยมเรียกย่อๆว่า SVD

การทำ SVD นั้นง่ายกว่าการคำนวณเมทริกซ์ความแปรปรวนร่วมเกี่ยวโดยตรง แต่สามารถได้ผลลัพธ์แบบเดียวกันนี้

หากใครใช้ฟังก์ชัน sklearn.decomposition.PCA ใน sklearn ที่จริงแล้วภายในนั้นก็คำนวณด้วย SVD เหมือนกัน

SVD คือการนำเมทริกซ์มาแยกเป็น ๓ ส่วน

X = UΣV*

โดยที่ X เป็นเมทริกซ์เดิม สมมุติว่ามี n×m (n แถว m หลัก)

U เป็นเมทริกซ์ขนาด n×n
Σ เป็นเมทริกซ์ขนาด n×m
V เป็นเมทริกซ์ขนาด m×m

โดยที่สัญลักษณ์ * หมายถึงเมทริกซ์สลับเปลี่ยนสังยุค (共轭转置, conjugate transpose)

ในที่นี้ขอมองข้ามเรื่องจำนวนจินตภาพไปเพราะใช้แต่เมทริกซ์ที่มีค่าเป็นจำนวนจริงอยู่แล้ว ดังนั้นไม่ต้องสนเรื่องการสังยุค (conjugate) แบบนี้ * ก็เป็นแค่ T คือเมทริกซ์สับเปลี่ยน (transpose) ธรรมดา คือแค่สลับแถวกับหลัก

ทั้ง U และ V เป็นเมทริกซ์เชิงตั้งฉาก (正交矩阵, orthogonal matrix)

Σ เป็นเมทริกซ์ที่มีค่าเฉพาะตัวในแนวทแยง ที่เหลือเป็น 0

ค่าในแนวทแยงของ Σ เรียกว่าเป็นค่าเอกฐาน (奇异值, singular value)

ก่อนจะอธิบายรายละเอียดมากกว่านี้ลองมาทำให้เห็นภาพได้ชัดขึ้นโดยการเขียนโปรแกรมคำนวณด้วยไพธอนจริงๆสักหน่อย

numpy มีฟังก์ชัน np.linalg.svd สามารถคำนวน SVD ได้ทันที ฟังก์ชันนี้มีอยู่ใน scipy ด้วย คือ scipy.linalg.svd จะใช้อันไหนก็ได้

ลองสร้างเมทริกซ์ขนาด 6×4 ขึ้นมาแล้วใช้ฟังก์ชันนี้ทำการคำนวณ SVD
import numpy as np
X = np.arange(-4,20).reshape(6,4)/2
U,sv,W = np.linalg.svd(X)
Sigma = np.zeros_like(X)
Sigma[:len(sv),:len(sv)] = np.diag(sv)
print('X=\n%s'%X)
print('U=\n%s'%U)
print('Σ=\n%s'%Sigma)
print('V*=\n%s'%W)

ได้
X=
[[-0.9330098  -1.2272592  -1.1769271  -0.71982633]
 [-2.00659432  2.39830423  0.60754422 -0.90707789]
 [-0.00334835  0.71909839  0.4503522  -1.11424495]
 [-0.25129929  1.69495859  0.49591245  2.07755182]
 [-0.60816189 -0.84689852  0.96693315 -1.37090087]
 [-0.91821157 -1.25774444 -0.40389328 -1.07915852]]
U=
[[ 0.40403946 -0.13452204  0.65586452  0.11819427 -0.195199   -0.58002167]
 [-0.4338776  -0.81503693  0.21835986  0.09622186  0.29604682  0.05368089]
 [-0.03959412 -0.3152936  -0.42657657  0.41880222 -0.72824513 -0.106388  ]
 [-0.63333423  0.17706933  0.29491611 -0.48910041 -0.48374194 -0.08563415]
 [ 0.28073123 -0.35767022 -0.41518736 -0.70846344  0.02561177 -0.34395498]
 [ 0.40867147 -0.2427489   0.28431862 -0.24543126 -0.33053053  0.72369647]]
Σ=
[[4.05145648 0.         0.         0.        ]
 [0.         3.40844577 0.         0.        ]
 [0.         0.         1.78739278 0.        ]
 [0.         0.         0.         1.06649963]
 [0.         0.         0.         0.        ]
 [0.         0.         0.         0.        ]]
V*=
[[ 0.02639931 -0.83676923 -0.23809852 -0.49237123]
 [ 0.63311362 -0.3250715  -0.18742535  0.67702837]
 [-0.6329527  -0.05263584 -0.67214738  0.38054993]
 [ 0.44479334  0.43746456 -0.67557289 -0.39291846]]

np.linalg.svd นั้นจะให้ค่าออก ๓ ตัว โดยตัวแรกก็คือ U

ส่วนตัวที่ ๒ ที่ได้คือค่าเอกฐาน (SV) ซึ่งเป็นอาเรย์มิติเดียวซึ่งเป็นค่าแนวทแยงของ Σ สามารถนำมาแปลงเป็นเมทริกซ์เป็นเมทริกซ์ Σ ทีหลังได้ โดยตัวที่เหลืือก็มีแต่ 0

ส่วนตัวที่ ๓ คือ V* ซึ่งในกรณีที่เป็นจำนวนจริงก็แค่เอาเมทริกซ์มาสับเปลี่ยนแถวกับหลัก คือ VT ในที่นี้ใช้ตัวแปรเป็น W แทน

เพื่อให้เห็นภาพชัดไปอีกแบบ ลองเอามาวาดแสดงด้วยสี ให้สีอ่อนแสดงค่ามาก สีเข้มคือค่าน้อย
import matplotlib.pyplot as plt
for i,(a,s) in enumerate(zip([X,U,Sigma,W],['X','U','Σ','V*'])):
    plt.figure(figsize=a.T.shape,dpi=50)
    plt.axes([0,0,1,1])
    plt.imshow(a,cmap='copper')
    plt.axis('off')
    plt.savefig('SVD_%s.png'%s)
    plt.close()

X



U



Σ



V*





V ที่ได้มานี้ก็คือเมทริกซ์ที่รวมเวกเตอร์ลักษณะเฉพาะ (本征向量, eigenvector) ขององค์ประกอบ

ส่วน SV ก็คือตัวที่แสดงถึงอัตราความบ่งบอกความแปรปรวน (解释方差占比, variance explained ratio) ของแต่ละองค์ประกอบ

แล้ว UΣ ก็จะเทียบเท่ากับค่าของ X ในพิกัดขององค์ประกอบนั่นเอง เพราะ (UΣ)VT = X และ XV = UΣ

จากตรงนี้จะเห็นว่าผลที่ได้จาก SVD นั้นเทียบเคียงได้กับสิ่งที่เราต้องการจากการวิเคราะห์องค์ประกอบหลักเหมือนตอนที่คำนวณโดยหาเวกเตอร์ลักษณะเฉพาะของความปรวนร่วม

เพียงแต่ว่าวิธีการคำนวณเพื่อให้ได้มาซึ่งผลลัพท์นั้นต่างกัน ความเร็วก็ต่างกันด้วย

เมื่อลองเปรียบเทียบดูแล้ว ถ้าจำนวนมิติของข้อมูลมีมากวิธีนี้จะคำนวณได้เร็วกว่า เพราะการคำนวณเมทริกซ์ความปรวนร่วมนั้นถ้ามิติมากจะต้องคำนวณสร้างเมทริกซ์ขนาดใหญ่เป็นจำนวนมิติกำลังสอง จำนวนมิติจึงมีผลต่อความเร็วอย่างมาก

ในขณะที่การคำนวณ SVD นั้นไม่ว่าจะนวนมิติมาก หรือจำนวนข้อมูลมาก ก็มีผลต่อความเร็วไม่ต่างกัน

ซึ่งที่จริงถ้าจำนวนมิติมีน้อยแล้วจำนวนของข้อมูลมีมาก วิธีการนี้กลับจะช้ากว่า



ตัวอย่างการใช้กับข้อมูลที่มีมิติมาก

ลองทดสอบกับข้อมูลที่มีมิติมากถึงหมื่นมิติดู โดยลดให้เหลือแค่ ๒ มิติแล้ววาดออกมาดู

ในตัวอย่างนี้จะสร้างคลาสในลักษณะเช่นเดียวกับที่สร้างไปใน https://phyblas.hinaboshi.com/20180727

เพียงแต่ว่าเปลี่ยนมาคำนวณด้วย SVD แทน
from sklearn import datasets

class WikhroOngprakopLak:
    def rianru(self,X):
        _,sv,W = np.linalg.svd(X)
        self.ver = sv/sv.sum()
        self.V = W.T

    def plaeng(self,X):
        return X.dot(self.V)

    def rianru_plaeng(self,X):
        self.rianru(X)
        return self.plaeng(X)

np.random.seed(2)
X,z = datasets.make_blobs(200,n_features=10000,centers=4,cluster_std=50)
wol = WikhroOngprakopLak()
wol.rianru(X)
Xi = wol.plaeng(X)[:,:2]
plt.gca(aspect=1)
plt.scatter(*Xi.T,c=z,cmap='rainbow',edgecolor='k')
plt.show()



จะเห็นว่าแม้จะมีมากถึงหมื่นมิติแต่ก็คำนวณออกมาได้ในเวลาไม่นานนัก แต่ถ้าลองเปลี่ยนไปใช้วิธีสร้างเมทริกซ์ความแปรปรวนร่วมจะต้องคำนวณสร้างเมทริกซ์ขนาด "หมื่นคูณหมื่น" ซึ่งนอกจากจะช้าแล้วอาจ RAM ไม่พอได้ง่าย



การคำนวณ SVD ของแบบเมทริกซ์มากเลขศูนย์

บ่อยครั้งที่เวลาคำนวณหาองค์ประกอบหลักนั้นเราต้องการลดข้อมูลจากที่มีมิติมากให้เหลือมิติน้อย

ปกติแล้วเวลาคำนวณ SVD เราจะได้เมทริกซ์ V ขนาดจำนวนแถวและหลักเท่ากับจำนวนมิติของข้อมูล เช่นมี m มิติก็จะได้เมทริกซ์ m×m

แต่ว่าสุดท้ายก็จะคัดทิ้งเหลือองค์ประกอบเท่ากับจำนวนมิติที่ต้องการให้เหลือ คืือ k ตัว นั่นคือใช้แค่เมทริกซ์ขนาด m×k เท่านั้น

หากรู้แน่อยู่แล้วว่าจะทิ้งเหลือกี่มิติ แบบนั้นถ้าใช้วิธีการคำนวณที่ไม่ต้องไปคำนวณในส่วนที่ยังไงเดี๋ยวก็จะต้องทิ้งอยู่แล้วเสียตั้งแต่แรกจะดีกว่า

ใน scipy มีวิธีที่จะทำแบบนั้นได้อยู่ คือใช้ฟังก์ชัน scipy.sparse.linalg.svds ในมอดูลย่อย scipy.sparse ซึ่งเป็นมอดูลย่อยที่เอาไว้จัดการเมทริกซ์ที่โล่งหร็อมแหร็ม

เมทริกซ์แบบนั้นเรียกว่า sparse matrix ส่วนศัพท์บัญญัติของไทยใช้คำว่า "เมทริกซ์มากเลขศูนย์" แม้จะไม่ได้แปลตรงตัว แต่ก็เป็นการแปลที่สื่อความหมายจริงๆ

sparse จริงๆถ้าแปลเป็นไทยก็คือ "หร็อมแหร็ม, เบาบาง, กระจัดกระจาย" ประมาณนั้น แต่ในที่นี้ก็คือหมายถึงเมทริกซ์ที่องค์ประกอบส่วนใหญ่เป็นเลข 0

รายละเอียดจะไม่พูดถึงมากตรงนี้ แต่สรุปง่ายๆก็คือ เมื่อใช้หลักการของเมทริกซ์มากเลขศูนย์ หรือเมทริกซ์หร็อมแหร็ม แล้วก็จะทำให้คำนวณ SVD ได้ด้วยความเร็วสูงขึ้นมาก

เราสามารถใช้ฟังก์ชัน scipy.sparse.linalg.svds แทน np.linalg.svd ได้

วิธีการใช้คล้ายๆกับ np.linalg.svd แต่ต่างกันตรงที่เราต้องกำหนดจำนวนมิติที่ต้องการเหลือไว้ตั้งแต่แรกเลย โดยใส่เป็นค่าอาร์กิวเมนต์ตัวที่ ๒ หรือจะใส่คีย์เวิร์ดเป็น k= ก็ได้

ตัวอย่าง ลองสร้างข้อมูลแสนมิติ แล้วคัดเหลือแค่ ๒ มิติ
from scipy.sparse.linalg import svds

np.random.seed(137)
X,z = datasets.make_blobs(200,n_features=100000,centers=5,cluster_std=80)

U,sv,W = svds(X,2)
Xi = X.dot(W.T)
plt.gca(aspect=1)
plt.scatter(*Xi.T,c=z,cmap='rainbow',edgecolor='k')
plt.show()



ถ้ามิติมากจริงๆยังไงก็ยิ่งช้า แต่ก็ถือว่าเป็นวิธีที่เร็วที่สุดแล้วหากต้องการวิเคราะห์องค์ประกอบของข้อมูลที่มีจำนวนมิติมากขนาดนั้น

เพียงแต่ถ้าข้อมูลไม่ได้มีขนาดใหญ่มาก การใช้ np.linalg.svd ก็ยังกลับจะเร็วกว่า

ถ้าใช้ฟังก์ชัน pca ของ sklearn นั้นภายในจะตรวจดูขนาดของข้อมูลให้อัตโนมัติ ถ้าข้อมูลใหญ่ก็จะใช้ sparse.linalg.svds ถ้าข้อมูลเล็กก็จะใช้ linalg.svd

ดังนั้นในการใช้งานจริงนั้นใช้ sklearn ไปเลยจะสะดวกที่สุด

สำหรับวิธีการใช้ sklearn เพื่อทำการวิเคราะห์องค์ประกอบหลักนั้นอ่านได้ใน https://phyblas.hinaboshi.com/20190921


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

囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧囧

ดูสถิติของหน้านี้

หมวดหมู่

-- คอมพิวเตอร์ >> ปัญญาประดิษฐ์
-- คอมพิวเตอร์ >> เขียนโปรแกรม >> python >> sklearn
-- คอมพิวเตอร์ >> เขียนโปรแกรม >> python >> matplotlib
-- คอมพิวเตอร์ >> เขียนโปรแกรม >> python >> numpy
-- คอมพิวเตอร์ >> เขียนโปรแกรม >> python >> scipy

ไม่อนุญาตให้นำเนื้อหาของบทความไปลงที่อื่นโดยไม่ได้ขออนุญาตโดยเด็ดขาด หากต้องการนำบางส่วนไปลงสามารถทำได้โดยต้องไม่ใช่การก๊อปแปะแต่ให้เปลี่ยนคำพูดเป็นของตัวเอง หรือไม่ก็เขียนในลักษณะการยกข้อความอ้างอิง และไม่ว่ากรณีไหนก็ตาม ต้องให้เครดิตพร้อมใส่ลิงก์ของทุกบทความที่มีการใช้เนื้อหาเสมอ

สารบัญ

รวมคำแปลวลีเด็ดจากญี่ปุ่น
python
-- numpy
-- matplotlib

-- pandas
-- pytorch
maya
การเรียนรู้ของเครื่อง
-- โครงข่าย
     ประสาทเทียม
javascript
บันทึกในญี่ปุ่น
บันทึกในจีน
-- บันทึกในปักกิ่ง
-- บันทึกในฮ่องกง
-- บันทึกในมาเก๊า
บันทึกในไต้หวัน
บันทึกในยุโรปเหนือ
บันทึกในประเทศอื่นๆ
เรียนภาษาจีน
qiita
บทความอื่นๆ

บทความแบ่งตามหมวด



ติดตามอัปเดตของบล็อกได้ที่แฟนเพจ

  ค้นหาบทความ

  บทความแนะนำ

ทำความเข้าใจระบอบประชาธิปไตยจากประวัติศาสตร์ความเป็นมา
เรียนรู้วิธีการใช้ regular expression (regex)
หลักการเขียนทับศัพท์ภาษาจีนกวางตุ้ง
การใช้ unix shell เบื้องต้น ใน linux และ mac
หลักการเขียนทับศัพท์ภาษาจีนกลาง
g ในภาษาญี่ปุ่นออกเสียง "ก" หรือ "ง" กันแน่
ทำความรู้จักกับปัญญาประดิษฐ์และการเรียนรู้ของเครื่อง
ค้นพบระบบดาวเคราะห์ ๘ ดวง เบื้องหลังความสำเร็จคือปัญญาประดิษฐ์ (AI)
หอดูดาวโบราณปักกิ่ง ตอนที่ ๑: แท่นสังเกตการณ์และสวนดอกไม้
พิพิธภัณฑ์สถาปัตยกรรมโบราณปักกิ่ง
เที่ยวเมืองตานตง ล่องเรือในน่านน้ำเกาหลีเหนือ
บันทึกการเที่ยวสวีเดน 1-12 พ.ค. 2014
แนะนำองค์การวิจัยและพัฒนาการสำรวจอวกาศญี่ปุ่น (JAXA)
เล่าประสบการณ์ค่ายอบรมวิชาการทางดาราศาสตร์โดยโซวเคนได 10 - 16 พ.ย. 2013
ตระเวนเที่ยวตามรอยฉากของอนิเมะในญี่ปุ่น
เที่ยวชมหอดูดาวที่ฐานสังเกตการณ์ซิงหลง
บันทึกการเที่ยวญี่ปุ่นครั้งแรกในชีวิต - ทุกอย่างเริ่มต้นที่สนามบินนานาชาติคันไซ
หลักการเขียนทับศัพท์ภาษาญี่ปุ่น
ทำไมจึงไม่ควรเขียนวรรณยุกต์เวลาทับศัพท์ภาษาต่างประเทศ
ทำไมถึงอยากมาเรียนต่อนอก
เหตุผลอะไรที่ต้องใช้ภาษาวิบัติ?

บทความแต่ละเดือน

2020年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2019年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2018年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2017年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

2016年

1月 2月 3月 4月
5月 6月 7月 8月
9月 10月 11月 12月

ค้นบทความเก่ากว่านั้น

ไทย

日本語

中文