Вступ до лінійних рявнянь для програмістів

Сучасному програмістові не доводиться самому писати алгоритми для розв'язку систем лінійних рівнянь. Цей вступ, власне, і не про деталі реалізацій. Втім, проблеми, які можна розв'язати у формі лінійних рівнянь, сучасному програмістові зустрічаються. Цей вступ, відповідно, саме про концептуальні основи, які допомогають такі проблеми розпізнати і вирішити.

Насправді, навіть якщо вас взагалі не збуджує лінійна алгебра, вам все одно можуть бути цікаві деякі ідеї про які тут йдеться. Сходимість, похибка обчислення, складність алгоритмів, — всі ці речі будуть наочно проілюстровані на прикладах простих лінійних систем.

Лінійна система — це система лінійних рівнянь. Рівняння називається лінійним, якщо воно не містить нічого окрім зваженої суми змінних що дорівнює константі. У двовимірному просторі лінійне рівняння працює ще як рівняння прямої. У тривимірному — площини. У чьотирьох і більше — гіперплощини.

Систему можна записати просто як множину рівнянь.

Цей графік інтерактивний. На ньому можна покерувати кривими, а системи нагорі автоматично перебудуються.

Графічне представлення — це така корисна ментальна модель. Про прямі думати легше ніж про числа і змінні. Наприклад, розв'язком системи рівнянь є перетин прямих. Дуже просто.

А буває навпаки, буває що кількість розв'язків нескінченна.

Чисельні методи зазвичай розраховані на пошук єдиного розв'язку і такі системи можуть бути для них проблемою. Втім, деякі імплементації проблемою це не вважають і повертають перший ліпший розв'язок. Що само по собі вже проблема.

Якщо помножити кожен коефіцієнт рівняння разом із константою на одне і те саме число, всі факти про точки, які це рівняння описує, лишаться незмінними. Це буде те саме рівняння. Воно буде рівнянням тієї ж самої прямої. Рівняння, які мають різні коефіцієнти, але описують одну і ту ж пряму, називаються лінійно залежними.

Система рівнянь може мати більше рівнянь ніж змінних. Такі системи називаються перевизначеними.

В загальному випадку, такі системи не мають розв'язку. Але в окремих випадках можуть і мати.

Може бути і навпаки. Система може мати більше змінних аніж рівнянь. Такі системи називаються недовизначеними.

Звісно, у двовимірному випадку, недовизначена система складатиметься з одного рівняння. Кожна точка прямої цього рівняння і буде розв'язком. Та чи система це взагалі?

Але у багатовимірних випадках, недоспецифіковані системи розв'язків можуть і не мати.

Є два класи алгоритмів які розв'язують системи лінійних рівнянь: ітеративні та точні.

Ітеративні алгоритми концептуально прості. Вони всі зводяться до трьох кроків.

Давайте зробимо з цього конкретний алгоритм.

Той крок, що наближує точку до розв'язку, буде простою проекцією на пряму. Точка, проекція, і розв'язок утворюватимуть прямокутний трикутник. Отже проекція буде принаймні не далі від розв'язку, аніж сама точка, тому що гіпотенуза не може бути коротшою за катет.

Ми будемо по черзі проецювати точку з однієї прямої на іншу, і рано чи піздно доберемося до якогось задовільно малого околу розв'язку.

Щодо умови виходу, ми зробимо ще простіше. Ми мірятимо те, скільки точка подолала за свій останній перехід щоб спроецюватися. Очікувано, із тим як ми наближатимемося до розв'язку, і переходи ставатимуть коротші. Отже, коли перехід стане зовсім маленьким, ми вважатемемо, що розв'язок достатньо близько аби ближче вже не підходити.

А от на питання «де починати» відповіді не буде. Чудово було б почати одразу біля розв'язку, але наш алгоритм має працювати починаючи з будь-якої точки, тож і багато думати про вибір першої не варто.

Ось інтерактивна ілюстрація такого алгоритму.

Якщо подібний до цього алгоритм і справді наближається до розв'язку із кожною ітерацією, то кажуть, що він сходиться. Якщо навпаки, віддаляється — розходиться. Сходимість — це друга за важливістю характеристика ітеративного алгоритму. Фактично це його швидкість.

А найголовніша властивість ітеративного алгоритму — це стабільність. Кожна операція на числах з плаваючою точкою може додавати якусь нехай невелику похибку. Якщо ми виконуватимемо такі операції багато разів підряд, похибка накопичуватиметься і зрештою може стати проблемою.

Стабільні алгоритми вміють не накопичуюти похибку. Нехай на кожній окремій ітерації якась похибка і присутня, але вона не впливатиме на похибку на наступній ітерації.

Цей конкретний алгоритм стабільний і геометрично не може розходитись. Більш того, його сходимість наочно пов'язана із кутом між прямими.

Втім, стабільність і сходимість — це ще не все. Цей конкретний алгоритм, наприклад, непридатний до практичного застосування з іншої причини. Якщо прямі майже паралельні, умова виходу яку ми обрали буде давати завелику похибку.

Точні алгоритми, на відміну від ітеративних, сходимості не мають. Натомість, вони виконуть свою работу за певну наперед визначену кількість операцій.

Але при деяких умовах, ітеративні алгоритми поводяться як точні. Теж виконуть роботу за наперед визначену кількість ітерацій.

Точні алгоритми працюють так незалежно від системи. Що найцікавіше, ми теж можемо перетворити наш простий ітеративний алгоритм на точний додавши всього одну невелику деталь

Ми робимо переходи послідовно проецюючи точку з прямої на пряму. Проекція відбувається ортогонально до наступної прямої. Давайте тепер замість цього проецювати точку на наступну пряму вздовж прямої поточної.

Раз алгоритм тепер не ітеративний, то і умова виходу йому не потрібна. А ще, так як ми маємо тепер «абсолютну сходимість», вибір першої точки теж ні на що не впливає. Хіба це не чудово?! Навіщо взагалі видумали ці ітеративні алгоритми?

Питання, звісно, саркастично-риторичне. Так само як практична ефективність ітеративного алгоритма залежить від сходимості, практична ефективність алгоритму точного залежить від його складності. Наприклад, цей конкретний алгоритм як раз непрактично складний.

Складність — це не те саме що просто кількість операцій. Як і сходимість, складність алгоритму може залежати від багатьох обставин. Інакше кажучи, змінних. Отже, це не число, а функція. Зазвичай, для того щоб записати функцію складності, ми використовуємо О-нотацію. Хитрість в тому, що у цій нотації ми описуємо вид трохи іншої функції — такої, яка обмежує складність алгоритму в залежності від певних параметрів.

Якщо, наприклад, треба оцінити цикл, який робить щось щонайбільше n раз, ми кажемо що його складність у О-нотації лінійна тобто O(n). Насправді, навіть якщо в деяких випадках цикл може закінчити роботу не доходячи до кінця, на оцінку обмежуючої функції це не впливає. Ми описуємо обмеження влгоритму а не його справжні числа. Це дозволяє нам мати уявленя про нього не занурюючись у тонкощі імплементації

Якщо цикл має в середині ще один цикл, то оцінка складності стає квадратичною: O(n²). Якщо це цикл всередині цикла всередині цикла: O(n³). Оцінка подібного роду називається поліноміальною, і часто зустрічається на практиці. Наприклад, наш новий алгоритм теж має поліноміальну оціночну складність.

def cross_of(A):
  ''' n-dimensional cross product on (n-1) vectors in list A '''
  D = len(A[0])
  N = len(A)

  v_res = [0.] * D
  for i in xrange(0, D):
    for jk in xrange(0, D ** N):
      v_ijk = [i] + [(jk/(D ** (N-j-1))) % D for j in xrange(0, N)]
      t_res = v__E(v_ijk)
      if t_res != 0:
        for k in xrange(0, N):
          t_res *= A[k][v_ijk[k + 1]]
        v_res[i] += t_res
  return v_res


def solution_for(A, B):
  p = [0. for each in B]
  for i in xrange(len(A)):
    plane_n = A[i]
    plane_d = -B[i]
    other_planes_ns = A[:i] + A[i+1:]
    projection_vector = cross_of(other_planes_ns)
    p = project_by_vector(p, projection_vector, plane_n, plane_d)
  return p

Як ви можете собі уявити, така оцінка робить алгоритм неефективним на великих системах. Для розв'язку системи із тисячею рівнянь нашому алгоритму знадобиться порядка триліона операцій. Це немало. Якщо збільшити систему ще на кілька порядків, то Сонце згасне до того як він взагалі закінчить свою работу.

От ітеративні алгоритми для того і придумали, щоб не чекати тисячі років там де можна почекати кілька хвилин і мати не такий точний, але все-таки розв'язок.

Це не означає, що вони завжди працюють швидше за точні. Просто вони надають здатність балансувати між точністью і швидкістю.

Отже, наші попередні алгоритми годяться тільки для того, щоб навчати на них сходимості і складності. На щастя, ліпші варіанти вже давно знайдені, описані, і реалізовані в коді. Наприклад, алгоритм Гауса, який має складність O(n³). При тому, що теоретична найліпша складність LU-розкладення взагалі O(n^2.376).

Головна ідея розкладання полягає в тому, що ми можемо змінювати систему як нам заманеться, поки ми робитимемо це за правилами, які не порушують істиності кожного окремого рівняння.

Ми вже бачили, що якщо помножити рівняння на якесь ненульове число, її відповідна пряма не змінюється. Рівняння містить ту саму інформацію, що і до множення.

А ще ми можемо додавати бідь-яке рівняння до будь-якого іншого. Звісно, це змінить його пряму, але не точку розв'язку, бо обидва рівняня мають однакову точку розв'язку за визначенням. Не неї додавання сусіднього рівняння матиме такий самий ефект як і множення. Ніякого.

Цих двох операцій цілком достатньо, щоб перетворити якусь довільну систему рівнянь на тривіальну.

P. S.

Навчальні алгоритми — це просто забавки. Як обрати практичний алгоритм у кожному конкретному випадку? Ось декілька підказок.

{	a11x1 + a12x2 = b1
	a21x1 + a22x2 = b2

{	a11x1 + a12x2 = b1
	a21x1 + a22x2 = b2

Вступ до лінійних рівнянь для програмістів

P. S.

Вступ до лінійних рівнянь
для програмістів