01204472/การทดลองเกี่ยวกับเมตริกซ์

หน้านี้เป็นส่วนหนึ่งของวิชา 01204472

เราจะทดลองเกี่ยวกับเมตริกซ์ และการใช้งานโครงสร้างข้อมูล matrix ของ numpy

ตัวอย่างการใช้งาน

>>> a = matrix([[1,2,3],[4,5,6]])
>>> a
matrix([[1, 2, 3],
        [4, 5, 6]])
>>> a.T
matrix([[1, 4],
        [2, 5],
        [3, 6]])
>>> norm(a,2)                              # หา matrix norm
9.5080320006957209
>>> b = matrix([[7,8],[9,10],[11,12]])
>>> b
matrix([[ 7,  8],
        [ 9, 10],
        [11, 12]])

>>> a*b
matrix([[ 58,  64],
        [139, 154]])

Orthogonal matrices

เราจะทดลองสร้าง orthogonal matrices โดยการหาเวกเตอร์ที่ตั้งฉากกัน กระบวนการดังกล่าวเรียกว่า Gram-Schmidt process

สุ่มเวกเตอร์

ใน pylab เราสามารถสุ่มเวกเตอร์ได้โดยใช้คำสั่ง rand [1] ซึ่งจะสามารถสร้าง array ตามมิติที่เราระบุได้ เช่น

>>> rand(5)
array([ 0.46074869,  0.45697852,  0.72675971,  0.87655621,  0.59247653])

เราสามารถสร้างเมตริกซ์จาก array ดังกล่าวได้โดยสั่ง matrix แต่เมตริกซ์ที่ได้จะเป็นเมตริกซ์ที่มี 1 แถว ไม่ใช่คอลัมน์เวกเตอร์ที่เราต้องการ แต่เราสามารถ transpose ได้โดยเรียก attribute T ของผลลัพธ์ดังกล่าว เช่น

>>> matrix(rand(5)).T
matrix([[ 0.50004005],
        [ 0.41567827],
        [ 0.56018141],
        [ 0.37370744],
        [ 0.29102686]])

จงเขียนฟังก์ชัน runit(n) ที่สุ่มเวกเตอร์ที่มีขนาด 1 หน่วยที่มีขนาด n

เวกเตอร์ที่ตั้งฉากกัน

ในส่วนนี้ เราจะสนใจกรณีของเวกเตอร์ขนาด 5 เท่านั้น ใช้ฟังก์ชัน runit() สุ่มเวกเตอร์ v1 จากนั้นสุ่มเวกเตอร์ u2 ให้ทดสอบว่าเวกเตอร์ u2 นั้นขนานกับ v1 หรือไม่ (มีโอกาสน้อยมากที่จะเกิดขึ้น) จากนั้นให้คำนวณหา

• เวกเตอร์ u2a และ u2b ที่ u2 = u2a + u2b และ u2a นั้นตั้งฉากกับ v1
• ให้เวกเตอร์ v2 เท่ากับเวกเตอร์ที่มีทิศทางเดียวกับ u2a แต่มีขนาด 1 หน่วย
  • ทดลองได้ว่า v2.T * v1 เท่ากับ 0 หรือไม่

จากขั้นตอนดังกล่าว เราจะสามารถสุ่มเวกเตอร์ v3 ที่ตั้งฉากกับทั้ง v1 และ v2 ได้

ให้เขียนฟังก์ชัน randorth(vlist,n) ที่รับรายการของเวกเตอร์ vlist ที่ประกอบไปด้วยเวกเตอร์ขนาด n ที่ตั้งฉากกันทั้งหมด และคืน unit vector ที่ตั้งฉากกับเวกเตอร์ทุกเวกเตอร์ในรายการ

ให้ใช้ฟังก์ชันดังกล่าว หาเวกเตอร์ v3, v4, และ v5 ที่เป็น unit vector ที่ตั้งฉากกันทั้งหมด

หมายเหตุ: ทดลองหา vector v6 ที่ตั้งฉากกับทุก ๆ เวกเตอร์

เมื่อเราได้เวกเตอร์ที่มีคุณสมบัติดังกล่าวแล้ว เราสามารถสร้างเมตริกซ์ Q ที่เป็น orthogonal matrix ได้ไม่ยาก ทดลองสร้างและตรวจสอบว่า Q เป็น orthogonal matrix จริงหรือไม่

Matrix norm

ใน pylab มีฟังก์ชัน norm (จาก numpy) สำหรับหา norm ของเมตริกซ์ อ่าน reference โดยในการหา operator norm ของเมตริกซ์ M เราจะสั่ง norm(M,2) นอกจากนี้ยังมี norm ประเภทอื่น เช่น ถ้าเราสั่ง norm(M) จะเป็นการหา Frobenius norm

เราจะหา norm ของเมตริกซ์ $A$ โดยใช้อัลกอริทึมสำหรับคำนวณหา eigenvalue ที่เรียกว่า Power method ทั้งนี้เนื่องจาก Operator norm นั้น มีค่าเท่ากับรากที่สองของ eigenvalue ที่มากที่สุดของ $A^{T}A$ (ทำไม?)

เราจะหาค่าดังกล่าว ให้เมตริกซ์ $B=A^{T}A$ เราจะสุ่มเวกเตอร์ $u$ จากนั้นเราจะนำเมตริกซ์ $B$ ไปคูณกับเวกเตอร์ $u$ ซ้ำไปเรื่อย ๆ (พร้อมกับสเกลเวกเตอร์ให้มีขนาด 1) นั่นคือ ในแต่ละรอบเราจะปรับค่า

$u\leftarrow {\frac {Bu}{\|Bu\|}}$

ทำไปสักระยะ เวกเตอร์ $u$ จะเริ่มไม่เปลี่ยนแปลง (จะลู่เข้าสู่ค่า eigenvector ที่สอดคล้องกับ eigenvalue ที่มีค่ามากที่สุด)

เราจะหาค่า eigenvalue ที่สอดคล้องกับ $u$ ได้โดยการคำนวณค่า $\|Bu\|/\|u\|$ ดังนั้นค่า norm ที่คำนวณได้จะเท่ากับ ${\sqrt {\|Bu\|/\|u\|}}$

ให้สุ่มเมตริกซ์ $A$ ดังนี้

>>> A = matrix(rand(3,3))
>>> A
matrix([[ 0.08684718,  0.51124141,  0.59892317],
        [ 0.70633162,  0.01742227,  0.64947786],
        [ 0.85831143,  0.67425728,  0.38012037]])

จากนั้นให้ $B=A^{T}A$ และคำนวณหาค่า eigenvalue ที่มากที่สุดของ $B$ โดยวิธีที่กล่าวมาแล้ว

เปรียบเทียบค่าที่ได้กับค่า norm ของเมตริกซ์ A (คำนวณโดย norm(A,2))

Range และ Nullspace ของเมตริกซ์

พิจารณาเมตริกซ์ $A$ ต่อไปนี้

$A=\left[{\begin{array}{ccc}1&0&1\\1&2&3\\-1&2&1\end{array}}\right]$

Range ของ $A$. ให้สุ่มเวกเตอร์ $x\in {\mathbf {R} }^{3}$ มาจำนวนมาก ๆ (เช่นสุ่มโดยใช้ matrix(rand(3)).T) จากนั้น plot ค่าของ $Ax$ ที่ได้เป็นแบบ scatter ในสามมิติ เพื่อแสดง range ของ $A$

Nullspace ของ $A$. เราจะประมาณ nullspace ของ $A$ โดยขั้นตอนดังนี้

• สุ่มเวกเตอร์ $x\in {\mathbf {R} }^{3}$ มาจำนวนมาก ๆ (เช่นสุ่มโดยใช้ matrix(rand(3)-[0.5,0.5,0.5]).T) เช่น ประมาณ 10000 เวกเตอร์
• นำเวกเตอร์ $x$ ที่ $\|Ax\|<\epsilon$ (โดยที่ค่า $\epsilon$ เป็นค่าน้อย ๆ เช่น 0.1) มา plot แบบ scatter ในสามมิติ เพื่อแสดง บริเวณเวกเตอร์ที่ใกล้กับ nullspace ของ $A$

การแก้ Linear equations

ทดลองการแก้สมการ

ในปัญหา polynomial interpolation เราทราบจุด $(x_{0},f(x_{0})),(x_{1},f(x_{1})),\ldots ,(x_{d},f(x_{d}))$ ของ polynomial $f$ ที่มี degree $d$ เราต้องการหา $f$

สมมติให้ $f(x)=a_{d}x^{d}+a_{d-1}x^{d-1}+\ldots +a_{1}x+a_{0}$ เราสามารถเขียนระบบสมการเชิงเส้นได้ดังด้านล่าง

TBA

เขียนฟังก์ชัน geneq(d,plist) ที่รับรายการของคู่ลำดับ d+1 คู่ จากนั้นคืนค่าเป็น tuple (M,b) ที่ $Ma=b$ เป็น linear equations ของสัมประสิทธิ์ของ $f$

ตัวอย่างการทำงาน

>>> geneq(2,[(1,5),(2,2),(4,7)])
(matrix([[1, 1, 1],
        [4, 2, 1],
        [16, 4, 1]]),
 matrix([[5],
        [2],
        [7]))

ให้ทดลองใช้ฟังก์ชัน solve เพื่อแก้สมการที่เราสร้างขึ้น ในหาทดลองหา polynomial ที่ผ่านจุดต่อไปนี้

• ข้อ a. ดีกรี 2, ผ่านจุด (1,5),(2,2),(4,7)
• ข้อ b. ดีกรี 4, ผ่านจุด (1,0), (3,5), (6,2), (9,7), (10,8)

ทดลองเขียน gaussian elimination

เขียนฟังก์ชันที่ทำงานลักษณะเดียวกับ ฟังก์ชัน solve โดยใช้ gaussian elimination

อย่าลืมนำฟังก์ชันดังกล่าวทดลองกับ การทดลองในแลบ 1 ก่อน

ทดลองใช้ฟังก์ชันที่เขียนขึ้นในการหา polynomial ในส่วนก่อนหน้า