Transformation Matrix: Eigenvalues & Eigenvectors "การคำนวณ"

พล่ามเรื่อง eigenvalue กับ eigenvector ไปเยอะละ แต่ยังไม่ได้บอกเลย ว่าเราจะหามันได้ยังไง :P คิดว่าถึงเวลาแล้วหละ มาลองดูกันดีกว่า

Analytical Solution

จำได้มั้ยว่า ใน Transformation Matrix: Eigenvalues & Eigenvectors "คุณสมบัติ" เรารู้มาอย่างนึงว่า

λ จะเป็น eigenvalue ของ A ก็ต่อเมื่อ det(A - λI) = 0

ดังนั้น เราก็ตั้งต้นด้วยสมการ

det(A - λI) = 0

สมการนี้ มักจะเรียกกันว่า "Characteristic Equation" ของปัญหาการหา eigenvalue และ eigenvector นี้ ส่วนพหุนาม P(λ) = det(A - λI) ก็จะเรียกว่า "Characteristic Polynomial"

ถ้า A มีมิติเป็น n × n เราจะรู้ว่า det(A - λI) เป็นพหุนามดีกรี n อย่างแน่นอน ดังนั้น เราจะสามารถหา eigenvalue λ ได้ n ตัว (นับตัวซ้ำด้วย)

คราวนี้ พอเราเอา λ แต่ละตัว ไปแทนใน

(A - λI)e = 0

จะได้ระบบสมการที่มี n สมการ และมีตัวแปร n ตัว ซึ่งมีคุณสมบัติ linear dependency (เพราะว่า det เป็น 0) ดังนั้น เซตของคำตอบจะมีจำนวนนับไม่ถ้วน ซึ่งก็ถูกต้อง เพราะว่า eigenvector มีจำนวนเป็นอนันต์

แต่ถ้าหาก λ ที่เลือกไปแทนในสมการ (A - λI)e = 0 นั้น ไม่ใช่รากซ้ำจากการแก้สมการ P(λ) = 0 เราจะรู้ว่า eigenvector ทั้งหมดที่ได้จะมีทิศขนานกัน

Eigenvalue ของเมตริกซ์สมมาตรที่มีสมาชิกเป็นจำนวนจริง เป็นจำนวนจริง

รู้สึกว่ากรณี A = A^T นี่จะพูดบ่อยจังเนอะ :D

วิธีพิสูจน์อันนี้ อาจจะดูแปลก ๆ หน่อย แต่ทุกขั้นตอนก็จำเป็นนะ เราจะเริ่มโดยสมมติว่า e เป็น eigenvector ที่อาจจะมีค่าที่ไม่ใช่จำนวนจริงก็ได้ เราจะพิสูจน์ว่า λ ซึ่งเป็น eigenvalue ของมัน ต้องเป็นจำนวนจริง

เราจะเขียน e^* เพื่อหมายถึง เมตริกซ์ conjugate transpose ของ e ซึ่งก็คือ เอาค่าในแต่ละมิติของ e^T ไปเปลี่ยนเป็น conjugate ของมัน

คราวนี้ ผลคูณภายในของเวกเตอร์เชิงซ้อน u กับ v ใด ๆ เราก็จะเขียนเป็น u^*v แทน u ⋅ v

เริ่มพิสูจน์โดยการหา e ⋅ Ae ...

e^*Ae = λ ||e||²

คราวนี้ ลองใส่ * เข้าไปทั้งสองข้าง

(e^*Ae)^* = (λ ||e||²)^*

ใช้คุณสมบัติที่ว่า (XY)^* = Y^*X^* กระจายเครื่องหมาย * ในฝั่งซ้าย ส่วนฝั่งขวา กระจายแค่ conjugate เข้าไปที่ λ เพราะว่า ||e||² เป็นจำนวนจริงอยู่แล้ว จะได้ว่า

e^*A^*e = λ^* ||e||²

แต่จากเงื่อนไขที่ว่า A เป็นเมตริกซ์จำนวนจริง และสมมาตร เราเลยรู้ว่า A = A^* ซึ่งทำให้

e^*Ae = λ^* ||e||²

ฝั่งซ้ายของสมการนี้ เท่ากับฝั่งซ้ายของสมการข้างบน (กลับไปดูสิ) พอจับมันมาเท่ากัน เราจะได้

λ ||e||² = λ^* ||e||²
λ = λ^*

สรุปได้แล้วว่า λ ต้องเป็นจำนวนจริง

ข้อจำกัดของ Analytical Solution

เนื่องจาก เราสามารถคำตอบเชิงวิเคราะห์ (Analytical Solution) ของสมการ P(λ) = 0 ได้ก็ต่อเมื่อ n ≤ 4 (ความจริงข้อนี้ พิสูจน์ยากมาก ๆ ... ขอยกมาเฉย ๆ ละกัน) ดังนั้น ในกรณีที่ n > 4 เราจะต้องหาวิธีทางตัวเลข (Numerical Method) เพื่อแก้สมการนี้

การหารากของสมการพหุนามแบบ numerical เนี่ย มันก็มีอยู่หลายวิธีแหละ แต่ ... จริง ๆ ถ้าเราจะหา numerical method อยู่แล้ว เราลองหาจากปัญหาต้นแบบได้มั้ย? แบบที่ไม่ต้องยุ่งกับ P(λ) หนะ

Numerical Solution

วิธีที่จะบอกเนี่ย เป็นหลักการเบื้องต้นนะ ในการใช้งานจริง ต้องมีการปรับปรุง เปลี่ยนแปลง อีกเยอะเลย

เริ่มละนะ ... สมมติว่าเรามีเวกเตอร์อะไรก็ไม่รู้ เส้นนึง ชื่อว่า u แล้วก็ สมมติว่า eigenvector ของเรา มีอยู่ n ตัว คือ e₁ e₂ e₃ ... e_n และ eigenvalue ก็คือ λ₁ λ₂ λ₃ ... λ_n ตามลำดับ

เนื่องจาก ตัวเลขห้อย e กับ λ มันจะเรียงกันยังไงก็ได้ เราก็สมมติซะว่า

λ₁ ≥ λ₂ ≥ λ₃ ≥ ... ≥ λ_n

ละกัน ... จากนั้น ก็สมมติว่า a_i เป็นผลจากการลองหาค่า Au แบบนี้

Au = a₁e₁ + a₂e₂ + a₃e₃ + ... + a_ne_n

สมมติว่า a₁ ≠ 0 ก่อนนะ คราวนี้ เราลองเอา A ใส่เข้าไปอีกที ทั้งสองข้างของสมการ ผลจะเป็น

A²u = λ₁a₁e₁ + λ₂a₂e₂ + λ₃a₃e₃ + ... + λ_na_ne_n

ถ้าเอา A ใส่ซ้ำไปเรื่อย ๆ มันจะเป็น ...

A^p+1u = λ₁^pa₁e₁ + λ₂^pa₂e₂ + λ₃^pa₃e₃ + ... + λ_n^pa_ne_n

คราวนี้ ถ้าหาก λ₁ > λ₂ ... ลองดึง λ₁ ออกมาจากสัมประสิทธิ์ทุกตัวสิ

A^p+1u = λ₁^p(a₁e₁ + (λ₂/λ₁)^pa₂e₂ + (λ₃/λ₁)^pa₃e₃ + ... + (λ_n/λ₁)^pa_ne_n)

พอเราให้ p → ∞ แล้ว เราจะเห็นว่า (λ_i/λ₁)^p เป็น 0 หมดทุกตัวยกเว้นเมื่อ i = 1 ดังนั้นข้อสรุปของเราก็คือ

A^p+1u ≈ λ₁^pa₁e₁

แสดงว่า การเอา A คูณซ้ำ ๆ ไปเรื่อย ๆ จะทำให้เราได้ e₁ ออกมาเอง

แต่ ถ้าหาก λ₁ = λ₂ = λ₃ = ... = λ_m หละ? (m ≤ n นะ) ... เราก็ดึงตัวร่วมออกจากทุกตัวได้อยู่ดีนะ มันจะเป็นงี้

จาก A^p+1u = λ₁^p(a₁e₁ + (λ₂/λ₁)^pa₂e₂ + (λ₃/λ₁)^pa₃e₃ + ... + (λ_n/λ₁)^pa_ne_n)

แต่ λ₁ = λ₂ = λ₃ = ... = λ_m > λ_m+1 (ถ้า m = n ก็ถือซะว่า λ_n+1 = 0 ไปเลย) ดังนั้น

A^p+1u = λ₁^p(a₁e₁ + a₂e₂ + a₃e₃ + ... + a_me_m
+ (λ_m+1/λ₁)^pa_m+1e_m+1 + (λ_m+2/λ₁)^pa_m+2e_m+2 + (λ_n/λ₁)^pa_ne_n)

ซึ่ง เมื่อ p → ∞ เราจะได้ว่า

A^p+1u ≈ λ₁^p(a₁e₁ + a₂e₂ + a₃e₃ + ... + a_me_m)

ผลลัพธ์ที่ได้เป็นเวกเตอร์ a₁e₁ + a₂e₂ + a₃e₃ + ... + a_me_m ซึ่งเรารู้ว่า มันก็เป็น eigenvector เหมือนกัน เพราะว่า eigenvalue มันเท่ากัน (พิสูจน์ไปแล้วในหัวข้อย่อย "จริง ๆ แล้ว จำนวน Eigenvector เป็นอนันต์" ในเรื่อง Transformation Matrix: Eigenvalues & Eigenvectors "คุณสมบัติ")

ดังนั้น เราจึงสรุปได้ว่า การเอาเวกเตอร์ที่เราสุ่มขึ้นมา คูณด้วย A ไปเรื่อย ๆ เราจะได้ eigenvector ตัวที่มีค่า eigenvalue สูงสุดออกมาเอง

มีอีกเรื่องที่ยังคาใจเราอยู่นะ คือ เมื่อตอนแรก เราสมมติว่า a₁ ≠ 0 ... แล้วถ้า a₁ = 0 มันจะเป็นยังไงน่ะเหรอ? เราก็คิดที่ a₂ แทนสิ ... ผลลัพธ์ก็คือ eigenvector ที่มี eigenvalue มากเป็นอันดับรองลงมาไงหละ

ดังนั้น ถ้าเรามั่วเวกเตอร์ u มาครั้งแรก แล้วเราหา e₁ ได้แล้ว เราก็มั่ว u มาอีกตัวที่ทำให้ a₁ = 0 ซะ คราวนี้ เราก็จะได้ e₂ แทน ... ทำซ้ำไปเรื่อย ๆ เราก็จะได้ e_i ครบทุกตัว

ในทางปฏิบัติ ทุกครั้งที่เราเอา A มาคูณ เราควรจะ normalize ผลลัพธ์ด้วย เพื่อให้ขนาดของเวกเตอร์ ไม่ใหญ่เกินไป จนเกิดปัญหา overflow

ผลพลอยได้จากปัญหา Eigenvector

มีคนเค้าพิสูจน์ไว้แล้วว่า ถ้าเราสมมติ P(λ) เป็นพหุนามดีกรี n ที่มีสัมประสิทธิ์เป็นจำนวนจริง ขึ้นมาอันนึง เราจะสามารถหา A ที่ทำให้ P(λ) = det(A - λI) ได้เสมอ

ดังนั้น ... แทนที่เราจะใช้ numerical method สำหรับหารากของสมการพหุนาม มาแก้ปัญหา eigenvector เราก็สามารถทำกลับกันได้ โดนเปลี่ยนปัญหาการหารากของสมการพหุนาม ไปเป็นปัญหา eigenvector แล้วใช้วิธีอย่างเมื่อกี๊

แปลว่า เราได้ numerical method สำหรับหารากของสมการพหุนาม เพิ่มขึ้นอีกวิธี ก็คือ วิธีการหา eigenvector นี่เอง

ทิ้งท้าย

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