การคำนวณเชิงตัวเลข (๑๒) LU decomposition ร่วมกับ Iterative improvement MO Memoir : Saturday 22 August 2558

Memoir ฉบับนี้จัดทำขึ้นเพื่อขยายความในเอกสารประกอบคำสอนวิชาคณิตศาสตร์สำหรับวิศวกรรมเคมี ในส่วนของการคำนวณเชิงตัวเลข เรื่องการแก้ปัญหาระบบสมการพีชคณิตเชิงเส้น

ในการแก้ปัญหาระบบสมการพีชคณิต A.x = b ด้วยเทคนิคการกำจัดของเกาส์ (Gauss elimination) นั้น จำนวนครั้งของการคูณในขั้นนตอนของการกำจัดไปข้างหน้า (Forward elimination) จะแปรผันตามจำนวนสมการยกกำลัง 2 ในขณะที่จำนวนครั้งของการคูณในขั้นตอนการแทนค่าย้อนกลับ (Back substitution) แปรผันตามจำนวนสมการยกกำลัง 1 หรือกล่าวอีกนัยหนึ่งคืองานส่วนใหญ่ไปหนักอยู่ที่ขั้นตอนการกำจัดไปข้างหน้า
  

ในกรณีที่เราต้องแก้ปัญหาระบบสมการพีชคณิต A.x = b หลายครั้ง โดยแต่ละครั้งนั้นมีเมทริกซ์ A ที่เหมือนกัน แต่แตกต่างกันตรงเวกเตอร์ b แม้ว่าเราจะสามารถใช้เทคนิคการกำจัดของเกาส์ในการแก้ปัญหาได้ แต่จะพบว่ามีการทำงานที่ซ้ำซ้อนกันอยู่ คือขั้นตอนการกำจัดไปข้างหน้าที่เราต้องทำการแปลงเมทริกซ์ A เพื่อให้ค่าตัวเลขในซีกซ้ายล่างของเมทริกซ์ A (พวกที่อยู่ใต้แนวเส้นทแยงมุม) มีค่าเป็น 0 และที่สำคัญคือขั้นตอนนี้เป็นขั้นตอนที่ใช้เวลาคำนวณมากที่สุด
 

ด้วยเหตุนี้จึงได้มีการหาวิธีที่จะทำการเก็บค่าตัวเลขที่ใช้ในการแปลงตัวเลขในตำแหน่งต่าง ๆ ที่อยู่ใต้เส้นทแยงมุมของเมทริกซ์ A ให้มีค่าเป็น 0 เพื่อที่จะได้สามารถนำไปใช้งานในครั้งต่อไปได้เลยโดยไม่ต้องทำการคำนวณใหม่ นั่นคือการแตกเมทริกซ์ A ให้กลายเป็นผลคูณของเมทริกซ์ 2 เมทริกซ์คือ L กับ U หรือ A = L.U

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

ดังนั้นจากโจทย์ของเรา A.x = b

จะแปลงเป็น L.(U.x) = b

กำหนดให้ U.x = y

ทำการแก้ระบบสมการ L.y = b เพื่อหาค่า y ก่อน

จากนั้นแก้ระบบสมการ U.x = y เพื่อให้ได้ค่า x

เพื่อให้ได้คำตอบที่ออกมาถูกต้อง (เท่าที่ความอุปกรณ์ที่ใช้ในการคำนวณจะให้ได้) การแก้ปัญหาด้วยเทคนิคการกำจัดของเกาส์นั้นจะเริ่มจากการทำการกำจัดไปข้างหน้าพร้อมกับการหาตัวหลัก (pivoting) ไปด้วยพร้อมกัน และในขณะเดียวกันก็จะทำการแตกเมทริกซ์ A ให้กลายเป็นผลคูณของเมทริกซ์ 2 เมทริกซ์คือ L กับ U เมื่อได้คำตอบ (คือ x) แล้วจะนำคำตอบที่ได้ไปแทนค่าในโจทย์เริ่มต้นเพื่อตรวจสอบว่าผลการคำนวณมีความคลาดเคลื่อนสะสมมากน้อยเท่าใด (ดูผลต่างระหว่างค่า b ที่โจทย์กำหนดมาให้กับค่าที่คำนวณได้จากการแทน x เข้าไปในสมการ) ถ้าพบว่าค่า b ที่ได้จากการแทนค่า x ที่คำนวณได้แตกต่างจากค่า b ที่โจทย์กำหนดมาให้ นั่นแสดงว่าค่า x ที่คำนวณได้นั้นมีความคลาดเคลื่อนสะสมอยู่ ก็จะทำการปรับแก้ด้วยการใช้เทคนิค iterative improvement
 

การปรับแก้คำตอบด้วยเทคนิค iterative improvement นั้นเป็นการแก้ปัญหาระบบสมการ A.dx = db เพื่อหาว่าค่าความคลาดเคลื่อนของคำตอบ (dx) นั้นมีค่าเท่าใด เพื่อที่จะได้ไปหักค่าความคลาดเคลื่อนนั้นออกจากค่าที่คำนวณได้ครั้งแรก กระบวนการปรับแก้ดังกล่าวเป็นกระบวนการที่ทำซ้ำไปเรื่อย ๆ จนกว่าจะพบว่าไม่สามารถปรับลดความคลาดเคลื่อนให้ต่ำลงไปได้อีกแล้ว (เกินความสามารถของเครื่องคำนวณ) และกระบวนการนี้ก็เป็นการทำให้การแก้ปัญหาระบบสมการ A.x = b ของเราเพียงข้อเดียว กลายเป็นการแก้ปัญหาระบบสมการ A.x = b หลายครั้ง โดยที่แต่ละครั้งนั้นมีค่าเวกเตอร์ b ที่แตกต่างกัน
  

เพื่อให้เห็นภาพเรามาลองดูตัวอย่างโดยใช้ระบบสมการพีชคณิต 5 สมการ 5 ตัวแปรที่เคยยกมาก่อนหน้านี้

0.001x₁ + 200x₂ + x₃ - 10x₄ - 200x₅ = 5     สมการที่ (1)

0.6x₁ + 5x₂ + 11.5x₃ + 4x₄ + 3x₅ = -3      สมการที่ (2)

-11x₁ + 3x₂ + 0.8x₃ + x₄ - 15.8x₅ = 7      สมการที่ (3)

2x₁ + 2x₂ + 2x₃ + 0.000067x₄ + 7x₅ = 6     สมการที่ (4)

-3x₁ - x₂ + 2x₃ + 8x₄ + 1.3x₅ = -11     สมการที่ (5)

จะทำการหาคำตอบของระบบสมการ (ค่า x₁, x₂, x₃, x₄ และ x₅) เปรียบเทียบกันระหว่างการใช้ Gauss elimination โดยไม่มีการหาตัวหลัก จากนั้นจะใช้เทคนิค Iterative improvement เพื่อทำการปรับแก้คำตอบ การคำนวณกระทำโดยใช้โปรแกรม Spreadsheet ของ OpenOffice 3.3.0
  

รูปที่ ๑ เป็นผลการคำนวณโดยไม่มีการหาตัวหลัก เมื่อนำค่า x₁, x₂, x₃, x₄ และ x₅ ที่คำนวณได้แทนค่ากลับเข้าไปในสมการที่ (1)-(5) พบว่าค่า b นั้นมีความคลาดเคลื่อนอยู่ (ตัวเลขสีแดงที่มุมขวาล่างของรูป) ดังนั้นจึงทำการแก้ระบบสมการ L.U.dx = db เพื่อหาค่า dx (ณ ขณะนี้เราไม่จำเป็นต้องทำการแตกเมทริกซ์ A ให้กลายเป็นผลคูณของเมทริกซ์ L กับ U เพราะเราได้ทำไปแล้วในระหว่างการหาคำตอบ x ชุดแรก) รูปที่ ๒ แสดงขั้นตอนคำนวณเพื่อคำนวณหาค่า dx ครั้งแรกเพื่อนำไปปรับคำตอบ x ที่ได้จากการคำนวณในครั้งแรก และเมื่อได้ค่า x ชุดใหม่แล้วก็ทำการทดสอบว่ายังมีความคลาดเคลื่อนอยู่หรือไม่ คือดูว่าค่า db เป็นศูนย์หรือไม่ ซึ่งก็พบว่ายังมีความคลาดเคลื่อนอยู่ (ค่า db ยังไม่เป็นศูนย์)
  

รูปที่ ๓ แสดงการทำ Iterative improvement ซ้ำเป็นครั้งที่สอง เพื่อหาว่าค่า x ที่ได้มาจากการปรับแก้ครั้งแรกนั้นยังมีความคลาดเคลื่อนหลงเหลืออยู่เท่าใด ซึ่งก็พบว่ามีค่าความคลาดเคลื่อนลดลง (จาก 10^-11-10^-12 ลงเหลือ 10^-16-10^-17) และเมื่อนำค่าคลาดคลาดเคลื่อนที่คำนวณได้ใหม่นี้ไปหักออกจากค่า x ที่ได้จากการปรับแก้ครั้งแรก ก็จะได้ค่า x ที่ผ่านการปรับแก้สองครั้ง แต่เมื่อนำค่า x ที่ได้จากการปรับแก้ครั้งที่สองนี้แทนค่ากลับเข้าไปใหม่กลับพบว่าค่าความคลาดเคลื่อน db นั้นไม่ลดต่ำลง ทั้งนี้เป็นเพราะค่า dx ที่ได้จากการคำนวณครั้งที่สองมีขนาดประมาณ 10^-16 เท่าของค่า x ที่ได้จากการปรับแก้ครั้งแรก และการคำนวณนั้นกระทำที่นัยสำคัญ 16 ตำแหน่ง ดังนั้นเมื่อนำค่า dx ที่ได้จากการคำนวณครั้งที่สองไปหักออกจากค่า x ที่ได้จาก การคำนวณครั้งแรก ผลลัพธ์ที่ได้จึงเป็นค่า x ชุดเดิม กล่าวอีกนัยหนึ่งคือเราไม่สามารถปรับคำตอบให้ถูกต้องได้มากกว่านี้อีกแล้วเนื่องด้วยข้อจำกัดของ machine precision ดังนั้นการคำนวณจะยุติลงเพียงแค่นี้

ไฟล์ที่แนบมาด้วยเป็นไฟล์ Spreadsheet ที่ใช้ในการคำนวณ