Розвʼязування лінійних рівнянь у

Про концепцію Розвʼязування лінійних рівнянь

Пригадаймо, як ми вперше вивчали алгебру ще в школі?

Зазвичай нам давали систему рівнянь і просили знайти x та y.

4x - 3y = -2
2x + 7y = 16

Два рівняння, дві невідомі: трохи підстановки, і ми швидко побачимо, що x = 1, y = 2.

Погляньмо ще раз на ліву частину попереднього прикладу. Вона має вигляд матриці [4 -3; 2 7], яка множить невідомий вектор [x, y] і дає вектор [-2, 16].

У позначеннях лінійної алгебри: A x = b, за звичною домовленістю писати матриці великими літерами, а вектори малими.

A x = 0 і нульовий простір

Важливий окремий випадок виникає тоді, коли в усіх наших рівняннях праворуч стоїть нуль.

4x - 3y = 0
2x + 7y = 0

У наведеному вище прикладі єдиний розвʼязок виникає за x = y = 0, а це зазвичай нецікаво.

Нетривіальні розвʼязки A x = 0 існують лише тоді, коли визначник A дорівнює нулю.

julia> using LinearAlgebra

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> det(A)  # not zero!
34.0

Натомість розгляньмо такі рівняння:

4x - 3y = 0
8x + 6y = 0

Друге рівняння просто вдвічі більше за перше і не додає нової інформації. Будь-які значення, де x = 0.75y, будуть розвʼязком.

Ці значення (x, y) лежать уздовж прямої у двовимірному просторі. Розвʼязків безліч, і вони утворюють нульовий простір матриці.

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

julia> A = [4 -3; 8 -6]
2×2 Matrix{Int64}:
 4  -3
 8  -6

julia> det(A)
0.0

julia> size(A)  # total number of rows and columns
(2, 2)

julia> rank(A)  # number of linearly-independent rows (or columns)
1

julia> nullspace(A)
2×1 Matrix{Float64}:
 -0.6
 -0.8

julia> nullspace(A) |> norm  # Julia returns unit vectors when possible
1.0

Функція nullspace() повертає одну точку, нормовану до одиничного вектора, але будь-яке скалярне кратне цього вектора теж належить нульовому простору.

Розвʼязування A x = b

А тепер уявімо, що у нас 1000 рівнянь із 1000 невідомими? Якщо це звучить безглуздо, згадаймо, що бездротові цифрові давачі тепер дешеві й універсальні (вимірюють деформацію, швидкість вітру, прискорення за трьома осями, та будь-що інше...). Цілком розумно мати 1000 таких давачів, що стежать за сучасною спорудою, як-от підвісний міст. Потрібна система керування, яка інтерпретує потік даних і готова подати сигнал тривоги, якщо ситуація стане загрозливою.

Знову A x = b, де в нас є A (з інженерного проєкту) і b (виміряні значення з давачів), але потрібно знайти x.

Перша думка, якщо це для нас нове, полягає в тому, щоб якось «поділити все на A» і перенести його в праву частину.

Насправді в A є обернена матриця (вона є в більшості квадратних матриць, але не в усіх: визначник має бути ненульовим), і обчислення працює (повільно!).

julia> A = [4 -3; 2 7]
2×2 Matrix{Int64}:
 4  -3
 2   7

julia> b = [-2, 16]
2-element Vector{Int64}:
 -2
 16

# the inverse matrix
julia> inv(A)
2×2 Matrix{Float64}:
  0.205882   0.0882353
 -0.0588235  0.117647

julia> inv(A) * b
2-element Vector{Float64}:
 0.9999999999999999
 2.0

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

На щастя, інші алгоритми значно швидші (у цьому випадку метод Гауса).

Julia (наслідуючи Matlab) просто використовує зворотну скісну риску для розвʼязувача (технічно це «ліве ділення»).

julia> x = A \ b
2-element Vector{Float64}:
 1.0
 2.0

Це тривіальний приклад, але є кілька деталей, за якими варто стежити.

• Коефіцієнти мають шикуватися в стовпці, тож додаймо нулі там, де бракує членів.
• Рівняння (а отже, і рядки матриці) мають бути лінійно незалежними. Якщо рядок 3 є просто сумою рядків 1 і 2, він не додає нової інформації, і задача недовизначена.

Команда rank() дасть кількість лінійно незалежних рядків/стовпців, тож прагнімо, щоб вона дорівнювала кількості невідомих змінних.

julia> rank(A)
2

Прямокутні матриці

Колись у дитинстві нам, мабуть, казали, що для N невідомих потрібна система з N рівнянь.

Як і в більшості речей у математиці, реальність трохи складніша (не лише деталь про лінійну незалежність, про яку йшлося в попередньому розділі).

«Замало» рядків

З N-1 рівняннями (рядками в матриці A) задача недовизначена. Невже опустимо руки з відчаю?

Як коли!

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

Це залишає нам безліч розвʼязків, але й безліч параметричного простору, де розвʼязків немає.

Чи проходить наша пряма розвʼязків через ділянку параметричного простору, яка мала б непокоїти порядного інженера? Або через таку, що виправдовує закриття споруди для відвідувачів (підвісні мости можуть добряче розхитуватися на сильному боковому вітрі)? Можливо, варто зробити ще кілька обчислень, щоб це зʼясувати!

julia> A3 = [4 -3 1; 2 2 2; 3 -4 -3]
3×3 Matrix{Int64}:
 4  -3   1
 2   2   2
 3  -4  -3

julia> b3 = [1, 12, -14]
3-element Vector{Int64}:
   1
  12
 -14

# fully-determined solution
julia> A3 \ b3
3-element Vector{Float64}:
 1.0
 2.0
 3.0

# remove row 3 from A3 and B3
julia> A2 = A3[1:2, :]
2×3 Matrix{Int64}:
 4  -3  1
 2   2  2

julia> b2 = b3[1:2]
2-element Vector{Int64}:
  1
 12

# an under-determined solution
julia> z1 = A2 \ b2
3-element Vector{Float64}:
 1.594594594594594
 2.445945945945947
 1.9594594594594592

Розвʼязувач зі зворотною скісною рискою дає нам розвʼязок! Неофіційно кажуть, що це розвʼязок із найменшою нормою (геометрично найближчий до початку координат).

Щоб отримати інші розвʼязки, нам потрібен nullspace для A2.

Додавання будь-якого скалярного кратного нульового простору до початкового розвʼязку дасть ще один дійсний розвʼязок.

# get the nummspace of A2
julia> N = nullspace(A2)
3×1 Matrix{Float64}:
 -0.46499055497527725
 -0.3487429162314578
  0.813733471206735

# add a random multiple of the nullspace to the earlier solution
julia> z = z1 + rand() * N
3×1 Matrix{Float64}:
 1.274331860650258
 2.205748895487695
 2.5199192438620472

# our random z is a valid solution, recovering b2 = [1, 12]
julia> A2 * z
2×1 Matrix{Float64}:
  0.9999999999999947
 12.0

Так само N-2 рівнянь обмежать (безліч) розвʼязків площиною, і діятимуть ті самі принципи.

«Забагато» рядків

Протилежна ситуація виникає, коли рівнянь більше, ніж невідомих, але вони все одно залишаються лінійно незалежними.

У реальній інженерії це цілком нормально, і це навіть добре!

Давачі мають обмежену точність, вектор b має обмежену точність, і в обчисленнях є шум. Тепер розвʼязком стає наближення методом найменших квадратів до зашумлених даних.

Найпростіший підхід використовує pseudoinverse матриці, який у Julia реалізовано як функцію pinv().

Його можна використовувати майже так само, як обернену матрицю (невиродженої) квадратної матриці, отримуючи оцінку змінних методом найменших квадратів, а не точний розвʼязок.

# create 5-row A and b for 2 variables
julia> A5 = [4 -3; 2 7; -1 2; 1 -1; 3 1]
5×2 Matrix{Int64}:
  4  -3
  2   7
 -1   2
  1  -1
  3   1

julia> rank(A5)
2

# add some random-normal noise to b5
julia> b5n = [-2, 16, 3, -1, 5] + randn(5) * 0.1
5-element Vector{Float64}:
 -1.9337349223581917
 16.024913008059457
  3.0087646177373992
 -0.8675748721176717
  4.914568047308805

# solve to get x ≈ 1, y ≈ 2, using the pseudoinverse
julia> pinv(A5) * b5n
2-element Vector{Float64}:
 1.006117942769543
 1.9962973764508272

Власні значення та власні вектори

У попередньому розділі йшлося про розвʼязки A x = b, що відповідають звичній алгебрі.

Цей розділ присвячено розвʼязкам A x = λ x, де λ позначає скаляр.

Це суто концепція лінійної алгебри, але її можна витлумачити геометрично:

Для відповідних значень x і λ квадратна матриця A масштабує ненульовий вектор x за довжиною у λ разів, не змінюючи його напрямку (хіба що відʼємне λ розвертає його).

Звучить вузько, але насправді це неймовірно корисно!

Термінологія: Допустимі значення λ називають власними значеннями A, а відповідні значення x називають власними векторами A. На жаль, доводиться миритися зі словами, які посередині змінюють мову.

Студентів зазвичай вчать обчислювати власні значення й вектори для матриць 2×2 вручну, але з компʼютером це значно простіше (хоча все ще досить повільно, і в загальному випадку масштабується як O(n^3) для матриці n×n).

julia> A = rand(-9:9, 2, 2)
2×2 Matrix{Int64}:
 9  5
 3  2

julia> F = eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274
vectors:
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

# eigenvalues
julia> F.values
2-element Vector{Float64}:
  0.27984674554472466
 10.720153254455274

# each column is an eigenvector, normalized to a unit vector
julia> F.vectors
2×2 Matrix{Float64}:
 -0.497417  0.945605
  0.867511  0.325317

Загалом матриця n×n має n власних значень, хоча вони не завжди різні. Уявляймо їх як корені многочлена n-го степеня (його називають характеристичним многочленом): вони можуть повторюватися й часто бувають комплексними навіть для матриці з дійсними значеннями.

Кожен власний вектор задає напрямок, і будь-яке його скалярне кратне теж є дійсним власним вектором. Для зручності подальших обчислень Julia повертає одиничні вектори з нормою, що дорівнює 1.

Застосування

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

На найвищому рівні власні вектори задають «найважливіші» осі (напрямки) в наборі даних. Власні значення вказують на «відносну важливість» кожної осі (за певних припущень про нормалізацію вхідних даних).

Що це насправді означає, залежить від застосування.

Метод головних компонент

У типовій задачі аналізу даних у нас може бути 100 або більше «ознак», збережених як стовпці даних. Майже неминуче будуть шум, надлишковість і небажані кореляції.

Нам потрібно зменшити розмірність, і PCA дає один зі способів упорядкувати цей хаос:

• Обчислімо коваріаційну матрицю даних.
• Тепер у нас є квадратна матриця, тож далі обчислімо власні значення та власні вектори.
• Упорядкуймо власні вектори за спаданням їхніх власних значень (точніше, за абсолютною величиною |λ|).
• Тепер перші k власних векторів утворюють principal components, де k значно менше за початкову кількість ознак.
• Спроєктуймо весь набір даних на ці k осей і почнімо шукати цікаві закономірності.

У такому вигляді PCA використовували найрізноманітніші групи: від медичних дослідників, які намагаються покращити властивості молекул у розробці ліків, до політичних агітаторів, які намагаються зрозуміти вподобання виборців за даними опитувань. Ймовірно, також маркетингові групи, які намагаються нам щось продати, але будь-яку технологію можна використовувати як на добро, так і на зло.

PCA не входить до мінімальної інсталяції Julia, але (поза Exercism) пакет MultivariateStats містить усе потрібне.

Обробка зображень

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

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

Дедалі частіше PCA є важливою складовою класифікаторів зображень, як-от розпізнавання облич. Наступного разу, коли ми проходимо через аеропорт, Старший Брат не просто спостерігає: він ще й використовує лінійну алгебру, щоб зрозуміти те, що бачить!

Машинобудування

Менш суперечливо: головні осі механічної деталі є власними векторами її тензора моменту інерції I (матриці під трохи іншою назвою).

Кожне колесо нашого автомобіля, ймовірно, має одну або кілька балансувальних вантажів, підібраних так, щоб обнулити позадіагональні елементи цього тензора I. Механік у гаражі, без сумніву, уникає цих обчислень (на відміну від авіаконструкторів та інженерів-ракетників), але «биття» лише побутова назва для перехресних членів тензора, які в цьому випадку роблять подорож менш комфортною та прискорюють механічний знос підшипників.