วิธีกำลังสองน้อยที่สุดที่มีข้อจำกัด

ฉันมีสมการเชิงเส้น 37 ตัวพร้อมตัวแปร 36 ตัวในรูปเมทริกซ์: A x = b (A มี 37 แถวและ 36 คอลัมน์) สมการไม่มีวิธีแก้ปัญหาที่แน่นอนดังนั้นฉันจึงใช้ Matlab เพื่อค้นหาคำตอบที่ใกล้เคียงที่สุดโดยใช้ x = A \ b

ปัญหาคือฉันมีเงื่อนไข องค์ประกอบทั้งหมดของ x ควรเป็นบวก: xi > 0 สำหรับ i ทั้งหมด x = A \ b ให้ค่าลบสำหรับองค์ประกอบบางอย่าง ฉันจะใช้ข้อจำกัดนี้ได้อย่างไร


นี่คือค่าที่เป็นรูปธรรมของ A และ b ที่ฉันทำงานด้วย:

A = [0.83   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0.02   0.63    0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0.02    -0.37   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0.02    -0.2    0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0.15   0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0.15    0   0   0   0   0.02    -0.33   0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.18    0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0   0   0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.52   0.17    0   0   0   0.18    0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.35   0   0   0   0   0.018   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0   -0.32   0.17    0   0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0   0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0   0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17    0
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.34   0.17
     0  0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0.15    0   0   0   0   0.02    -0.17
     1  1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1];
b = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]';

person nasim    schedule 20.05.2015    source แหล่งที่มา
comment
ถ้าคุณมีสมการ 37 สมการสำหรับสมการที่ไม่ทราบ 36 สมการ นั่นเป็นระบบที่มีการกำหนดไว้มากเกินไป ซึ่งหมายความว่าไม่น่าจะมีวิธีแก้ปัญหาหลายข้อ (และจริงๆ แล้ว มีแนวโน้มสูงที่จะไม่มีวิธีแก้ปัญหา)   -  person    schedule 20.05.2015
comment
ฉันรู้ แต่อย่างน้อยฉันก็ไม่สามารถหาวิธีแก้ปัญหาโดยประมาณที่มีข้อผิดพลาดน้อยที่สุดได้ใช่ไหม   -  person nasim    schedule 20.05.2015
comment
หากคุณต้องการค้นหา X0 ที่ลดข้อผิดพลาดบางอย่างให้เหลือน้อยที่สุด คุณต้องกำหนดข้อผิดพลาดนั้นก่อน :-) คุณจะใช้บรรทัดฐานประเภทใดสำหรับข้อผิดพลาด E = B - A*X0 ยุคลิด, รถแท็กซี่, สูงสุด?   -  person    schedule 20.05.2015
comment
ฉันไม่รู้คำศัพท์พวกนั้นด้วยซ้ำ! (ฉันกำลังอ่านเกี่ยวกับพวกเขาตอนนี้) ตัวแปรของฉันคือความน่าจะเป็นในการเปลี่ยนแปลงในเครือข่ายที่เข้าคิวและโดยปกติจะอยู่ในช่วงระหว่าง 0.0001 ถึง 0.1 คุณแนะนำอันไหน?   -  person nasim    schedule 20.05.2015
comment
@CST-Link ชื่อบอกว่ากำลังสองน้อยที่สุด ซึ่งระบุฟังก์ชันต้นทุน   -  person A. Donda    schedule 20.05.2015
comment
@nasim คุณช่วยโพสต์ค่าที่แน่นอนของ A และ b ได้ไหม? นี่จะทำให้การดีบักคำตอบที่เป็นไปได้ง่ายขึ้น ในระหว่างนี้ โปรดดูวิธีแก้ปัญหาของฉัน   -  person A. Donda    schedule 21.05.2015
comment
A จะไม่พอดีกับความคิดเห็น :(   -  person nasim    schedule 21.05.2015
comment
@nasim ไม่ ไม่ใช่ความคิดเห็น แก้ไขคำถามและรวมเป็นรหัส   -  person A. Donda    schedule 21.05.2015
comment
@nasim ขอบคุณ แต่มันไม่สมเหตุสมผล A นี่คือ 36x36 และ b นี่คือ 1x38 โปรดแก้ไขปัญหานี้   -  person A. Donda    schedule 21.05.2015


คำตอบ (3)


คุณต้องการหาคำตอบโดยประมาณ x ถึง A x = b ในความหมายกำลังสองน้อยที่สุด เช่น คุณต้องการย่อให้เล็กสุด

||A x - b||2 = xT AT A x + bT b - 2 xT AT ข.

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

1/2 xT สูง x + fT x,

ถ้าเราเลือก H = AT A และ f = - AT b

การใช้ quadprog ที่สอดคล้องกันคือ:

H = A' * A;
f = - A' * b;
x = quadprog(H, f)

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

n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))

ข้อจำกัดด้านบวกไม่สมเหตุสมผล เพราะถ้าคำตอบที่ดีที่สุดเกี่ยวข้องกับองค์ประกอบของ x หนึ่งองค์ประกอบขึ้นไปให้เป็น 0 พอดี ดังนั้นคำตอบเชิงบวกอย่างเคร่งครัดจะดีกว่ายิ่งองค์ประกอบที่สอดคล้องกันมีขนาดเล็กลง: 0.01 จะดีกว่า 0.1, 0.001 จะดีกว่า จะดีกว่า 0.01 ฯลฯ เป็นต้น – ไม่มีขอบเขตตามธรรมชาติ หากคุณต้องการแน่ใจว่าผลเฉลยเป็นบวกทั้งหมด คุณต้องกำหนดขอบเขตจำกัดด้วยตัวเอง:

x = quadprog(H, f, -eye(n), zeros(n, 1) + 0.001)

ตอนนี้ค่าที่น้อยที่สุดที่เป็นไปได้ขององค์ประกอบของ x คือ 0.001


อัปเดต หลังจากเสริมคำถามด้วยข้อมูลจริงของ A และ b: การใช้โค้ด

H = A' * A;
f = - A' * b;
n = size(A, 2);
x = quadprog(H, f, -eye(n), zeros(n, 1))

ฉันได้รับผลลัพธ์:

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the function tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

x =
      0.000380906335150292
      3.90638261088393e-05
        0.0111196970167585
        0.0227055107206744
        0.0318402514628274
        0.0371743514880516
      0.000800900221354844
       0.00746652476710186
        0.0180511534370576
        0.0282423767946842
        0.0362606972021829
        0.0417582260990786
       0.00860220929402253
        0.0174105435824309
        0.0265771677458008
        0.0343071472371469
        0.0395176470725881
        0.0419494410289298
        0.0187719294637544
        0.0268976053211278
        0.0336818044612046
        0.0382365751296441
        0.0398823076542831
        0.0391016682549663
        0.0279383031707377
        0.0339393563379992
        0.0377917413001034
        0.0382731422972829
        0.0338557405807941
        0.0217568643500703
        0.0343698083354502
        0.0381554349806972
        0.0392353941260779
        0.0368010570888738
         0.031271868401718
        0.0258232230013864
person A. Donda    schedule 20.05.2015
comment
คำเตือน: อัลกอริทึม Trust-region-reflective ไม่สามารถแก้ปัญหาประเภทนี้ได้ โดยใช้อัลกอริทึมชุดที่ใช้งานอยู่ หากต้องการความช่วยเหลือเพิ่มเติม โปรดดูการเลือกอัลกอริทึมในเอกสารประกอบ - person nasim; 21.05.2015
comment
@nasim ฉันตรวจสอบรหัสนี้ด้วยข้อมูลที่สร้างขึ้นแบบสุ่ม ดังนั้นฉันจึงไม่สามารถตรวจสอบคำเตือนนี้ได้ด้วยตนเอง ตามที่ฉันแสดงความคิดเห็นกับคำถาม ดีกว่าถ้าคุณโพสต์เมทริกซ์ของคุณ - person A. Donda; 21.05.2015
comment
ฉันจะยินดี แต่เมทริกซ์ A ไม่พอดีกับความคิดเห็น ฉันสามารถโพสต์ได้ที่ไหน? - person nasim; 21.05.2015
comment
@nasim คลิกลิงก์แก้ไขด้านล่างคำถามของคุณ และคัดลอกข้อมูลลงในคำถาม - person A. Donda; 21.05.2015
comment
@nasim ขอบคุณสำหรับข้อมูลที่แก้ไข แต่ด้วยค่าเหล่านี้ ฉันจึงได้วิธีแก้ปัญหาที่เหมาะสมโดยไม่มีการเตือนล่วงหน้า วิธีแก้ปัญหาจะเป็นค่าบวกทั้งหมดแม้ว่าจะไม่ได้ระบุขอบเขตจำกัดก็ตาม x ที่เล็กที่สุดคือ 3.9e-5 - person A. Donda; 21.05.2015
comment
จริงหรือ?! ฉันขอดูค่า 36 ค่าได้ไหม และฉันขอรหัสที่คุณใช้เป็นคำตอบได้ไหม - person nasim; 21.05.2015
comment
ให้เราสนทนาต่อในการแชท - person nasim; 21.05.2015

นี่อาจเป็นปัญหาการหาค่าเหมาะที่สุดเชิงเส้นเมื่อคุณมีเงื่อนไข x>0 อัลกอริธึมที่ดีที่สุดในการแก้ปัญหาคืออัลกอริธึมแบบซิมเพล็กซ์ แนวคิดก็คือว่าสมการเชิงเส้นแต่ละสมการ aix=bi มีเส้นตรง และการรวมกันของเส้นเหล่านี้ทำให้เกิดรูปหลายเหลี่ยม/รูปทรงหลายเหลี่ยม และคำตอบคือหนึ่งในจุดยอดของรูปทรงหลายเหลี่ยม/หลายเหลี่ยมนี้ อัลกอริธึม Simplex ค่อนข้างเป็นมาตรฐานและมีฟังก์ชันและไลบรารีมากมายที่สามารถคำนวณได้

person Good Luck    schedule 20.05.2015

คุณสามารถใช้ฟังก์ชัน Matlab lsqlin หรือ lsqnonneg เพื่อแก้ไขปัญหาของคุณได้ เช่น:

x=lsqnonneg(A,b)

จะให้สิ่งที่คุณกำลังมองหา

person BayesianMonk    schedule 29.08.2018