Революція роботів

Що таке «лінійна алгебра»?

Існує чимало технічних визначень: у [Вікіпедії][linalg-wiki], у підручниках і на чудових сайтах, як-от [3Blue1Brown][3blue1brown], але для наших цілей можна бути менш формальними.

Супровідники Julia (і Python) обожнюють такі речі. А от керівництво Exercism - не дуже.

Ця концепція обмежиться простішими аспектами величезної теми, але застереження: це неминуче досить математична концепція.

Що таке вектор?

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

• У Julia Vector{T} - це лише псевдонім для Array{T, 1}, а eltype T може бути будь-чим.
• У (більшій частині) фізики вектор - це стрілка з length та direction у N-вимірному просторі, але без фіксованого положення.
• У лінійній алгебрі стрілка фізиків зафіксована в просторі, а її хвіст розміщено в початку координат. Через це стрижень стрілки стає зайвим, тож ми можемо подати вектор як point у N-вимірному просторі, де N елементів задають відстань від початку координат уздовж кожної з N осей.

Для наших поточних цілей знехтуймо векторами, складеними з рядків тексту (англ. string) чи символів. У цій концепції ми будемо працювати з числовими типами: Int, Float чи Complex.

Що таке матриця?

Тут теж відповіді бувають різні.

• Прямокутний двовимірний масив чисел (без рваних країв). Квадратна матриця - поширений окремий випадок.
• linear combination стовпцевих векторів, розміщених один біля одного.
• transformation, який можна застосувати до вектора, так само як function застосовують до інших типів вхідних даних.

Деякі типи квадратних матриць настільки поширені, що мають власні назви.

Діагональна матриця: усі ненульові елементи розташовані на main diagonal (від верхнього лівого кута до нижнього правого).

julia> [1 0 0; 0 2 0; 0 0 3]
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

# a convenient shortcut
julia> diagm(1:3)
3×3 Matrix{Int64}:
 1  0  0
 0  2  0
 0  0  3

Одинична матриця: діагональна матриця, на діагоналі якої стоять лише одиниці (з причин, які стануть зрозумілішими далі). Часто скорочується як I.

julia> Matrix{Float64}(I, 3, 3)
3×3 Matrix{Float64}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

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

Транспонування матриці

Якщо справді потрібно поміняти рядки місцями зі стовпцями, функція permutedims() зробить це і поверне нову матрицю.

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

Для цілей лінійної алгебри функція transpose() корисніша, бо швидко створює ліниву обгортку навколо вихідної матриці.

Ще корисніша функція adjoint(), яка до того ж змінює знак уявної частини комплексних чисел. Ця операція настільки поширена, що достатньо додати апостроф ' до назви змінної, щоб отримати adjoint.

julia> m
2×3 Matrix{Int64}:
 1  2  3
 4  5  6

# new, full copy
julia> permutedims(m)
3×2 Matrix{Int64}:
 1  4
 2  5
 3  6

# lazy version
julia> transpose(m)
3×2 transpose(::Matrix{Int64}) with eltype Int64:
 1  4
 2  5
 3  6

julia> mc = [1+2im 2+3im; 3+2im 1+2im]
2×2 Matrix{Complex{Int64}}:
 1+2im  2+3im
 3+2im  1+2im

# lazy conjugate transpose
julia> adjoint(mc)
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

# syntactic sugar for adjoint
julia> mc'
2×2 adjoint(::Matrix{Complex{Int64}}) with eltype Complex{Int64}:
 1-2im  3-2im
 2-3im  1-2im

Множення

Решта цього документа використовує можливості модуля LinearAlgebra. Рядок using LinearAlgebra додасть його до простору імен, але ми не будемо повторювати це в прикладах (занадто багато візуального шуму).

Вектори, поелементний добуток

Про це йшлося в концепції [операцій з векторами][vector-ops]. Оператор - це .*; він діє на вхідні вектори попарно й дає вихідні дані того самого розміру й типу, що й вхідні.

julia> [1, 2] .* [3, 4]
2-element Vector{Int64}:
 3
 8

Вектори, скалярний добуток

Цю надзвичайно поширену операцію в підручниках записують як u ⋅ v і називають скалярним, або «dot», добутком.

Скалярний добуток дорівнює сумі поелементного добутку.

Два вектори можна перемножити звичайним оператором *, але лише якщо лівий вектор є adjoint: зручно записувати це як u' * v. Це стане зрозумілішим у пізнішому розділі про множення матриць.

Підняту (середню) точку можна використати в Julia (вводиться як \cdot і Tab) як синтаксичний цукор для функції dot(). Вказувати adjoint не потрібно, бо ця деталь обробляється автоматично.

julia> using LinearAlgebra

# with element-wise syntax
julia> sum([1, 2] .* [3, 4])
11

# with u' * v syntax
julia> [1, 2]' * [3, 4]
11

# with dot()
julia> dot([1, 2], [3, 4])
11

# with \cdot syntax
julia> [1, 2] ⋅ [3, 4]
11

Для векторів із комплексними значеннями лівий вектор має бути спряженим (зі зміненим знаком уявної частини). Функція dot робить це автоматично, синтаксис u' робить це явно, але sum(u .* v) не спрацює, якщо u і v комплексні.

Норми

Наскільки «великий» вектор?

norm - це спроба це вловити, звівши вектор до відповідного скаляра.

Визначено цілу родину норм, але найпоширеніша з них - 2-норма, яка дорівнює √(v ⋅ v).

Ця середньоквадратична операція - це піфагорова відстань від початку координат (у N-вимірному просторі). Якщо уявити вектор як стрілку з хвостом у початку координат, то 2-норма - це довжина цієї стрілки.

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

# defaults to the 2-norm
julia> norm([1, 2, 3])
3.7416573867739413
# the 1-norm
julia> norm([1, -2, 3], 1)
6.0

Множення матриць

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

Це операція, у якій компʼютери дуже сильні:

• Вона надзвичайно повторювана.
• Її можна ефективно розпаралелити.
• Для її прискорення розробили багато спеціалізованого обладнання, зокрема графічний процесор (GPU), який, найімовірніше, вбудовано в ноутбук.

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

Розгляньмо матрицю A, яка множиться на вектор v і дає вихідні дані w (за домовленістю лінійної алгебри матриці позначаємо великими літерами, а вектори - малими).

Верхній рядок A скалярно множиться на v і дає перший елемент w, другий рядок дає другий елемент, і так далі вниз.

julia> A = [1 2; 3 4]
2×2 Matrix{Int64}:
 1  2
 3  4
julia> v = [5, 6]
2-element Vector{Int64}:
 5
 6
julia> A * v
2-element Vector{Int64}:
 17  # equals [1, 2] ⋅ [5, 6]
 39  # equals [3, 4] ⋅ [5, 6]

Множення матриці на матрицю поширює це на всі стовпці матриці праворуч.

Для C = A * B можна вважати, що A множиться на кожен стовпець B і дає відповідний стовпець C: це серія множень матриці на вектор.

Рівнозначно можна сказати, що перший рядок A скалярно множиться на кожен стовпець B і дає верхній рядок C, другий рядок дає другий рядок, і так далі вниз.

julia> A
2×2 Matrix{Int64}:
 1  2
 3  4

julia> B = [5 7; 6 8]
2×2 Matrix{Int64}:
 5  7
 6  8

# left column is the same as A*v previously
julia> A * B
2×2 Matrix{Int64}:
 17  23
 39  53
 
julia> B * A
2×2 Matrix{Int64}:
 26  38
 30  44

Як видно з прикладу вище, множення матриць не комутативне: між A*B і B*A немає простої залежності.

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

Розмірності

Скалярний добуток двох векторів можливий лише за однакової довжини.

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

Якщо подати розміри як кортежі (nrows, ncols), які повертає size(A), маємо (a, b) * (b, c) -> (a, c). «Внутрішні» розмірності, тут b і b, сумісні для обчислення скалярних добутків. «Зовнішні» розмірності, тут a і c, визначають розмірність вихідних даних.

Приклад із прямокутними матрицями:

julia> D = reshape(1:6, 2, 3)
2×3 reshape(::UnitRange{Int64}, 2, 3) with eltype Int64:
 1  3  5
 2  4  6

julia> E = reshape(1:12, 3, 4)
3×4 reshape(::UnitRange{Int64}, 3, 4) with eltype Int64:
 1  4  7  10
 2  5  8  11
 3  6  9  12

julia> D * E
2×4 Matrix{Int64}:
 22  49   76  103
 28  64  100  136

# (3, 4) * (2, 3) not possible
julia> E * D
ERROR: DimensionMismatch: matrix A has axes (Base.OneTo(3),Base.OneTo(4)), matrix B has axes (Base.OneTo(2),Base.OneTo(3))

Множення пар векторів у матричному стилі має два варіанти.

За домовленістю ми використали б u' * v як еквівалент скалярного добутку. Розмірності тут (1, 3) * (3, 1), і Julia спрощує вихід (1, 1) до скаляра (на відміну, наприклад, від R).

Або ж можна використати u * v' з розмірностями (3, 1) * (1, 3) -> (3, 3), щоб перемножити всі можливі пари елементів і розгорнути вектори в матрицю.

julia> u = [1, 2]
2-element Vector{Int64}:
 1
 2
julia> v = [3, 4]
2-element Vector{Int64}:
 3
 4
julia> u' * v
11
julia> u * v'
2×2 Matrix{Int64}:
 3  4
 6  8

u' * v іноді називають внутрішнім добутком.

Обертання

Один особливо поширений різновид множення матриць повʼязаний з матрицями обертання.

У двовимірному просторі є відносно проста матриця, яка обертає вектор проти годинникової стрілки на θ радіан.

julia> rot2d(θ, vec) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * vec
rot2d (generic function with 1 method)

# unit vector in the x direction
julia> i_hat = [1, 0]
2-element Vector{Int64}:
 1
 0

# rotate 45 degrees
julia> rot2d(π/4, i_hat)
2-element Vector{Float64}:
 0.7071067811865476
 0.7071067811865475

# rotate 90 degrees -> unit vector in the y direction, j_hat
julia> rot2d(π/2, i_hat)
2-element Vector{Float64}:
 6.123233995736766e-17  # zero, within numerical error
 1.0