Основи лінійної алгебри у

Про концепцію Основи лінійної алгебри

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

Існує безліч технічних означень: у Вікіпедії, у підручниках і на чудових сайтах, як-от 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(), яка до того ж змінює знак уявної частини в комплексних числах (якщо цікаво, чому це так важливо, коротка відповідь - квантова механіка). Ця операція настільки поширена, що достатньо додати апостроф ' до назви змінної, щоб отримати спряження.

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 додасть його до простору імен, але ми не будемо повторювати це в прикладах (забагато візуального безладу).

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

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

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

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

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

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

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

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

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 комплексні.

Вектори, векторний добуток

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

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

У геометричному поданні векторного простору новий вектор перпендикулярний до площини, що містить обидва вхідні вектори. Якщо вхідні вектори паралельні (з точністю до знака), вони не задають площини, тож результатом буде [0, 0, 0].

Порядок має значення: u × v == -(v × u). Це відоме правило правої руки, через яке багато з нас витріщалися на власний великий палець і два пальці, крутячи їх у просторі (часто зі спантеличеним виразом).

# with \times syntax
julia> [1, 2, 3] × [3, 4, 5]
3-element Vector{Int64}:
 -2
  4
 -2
# with cross()
julia> cross([1, 2, 3], [3, 4, 5])
3-element Vector{Int64}:
 -2
  4
 -2

Зауважмо, що ця операція обмежена векторами довжини 3 (що відповідає евклідовому простору з ортогональними осями x, y, z). Математики радо визначають цю операцію й для інших вимірів, отримуючи в результаті незрозумілий набір слів (тензори вищого порядку, клиновий добуток, зовнішній добуток, мультивектори...). У цей момент більшість із нас тікає або швидко змінює тему!

Норми

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

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 років на те, щоб перетворити кожне обчислення на послідовність множень матриць.

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

• Вона дуже повторювана.
• Її можна ефективно розпаралелити.
• Було розроблено багато спеціалізованого обладнання, щоб прискорити її: від Cray-1 вартістю 80 мільйонів доларів у 1970-х до графічного процесора, який, імовірно, вбудовано в ноутбук.

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

Розгляньмо матрицю 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, другий рядок дає другий рядок, і так далі вниз.

Існує кілька таких ментальних уявлень, і зацікавлені студенти можуть подивитися цілу лекцію MIT, присвячену їм.

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 є багато відео, які показують його графічно, тож пошукайте «matrix multiplication» і виберіть те, що відповідає бажаному стилю, рівню деталізації та мові.

Розмірності

Скалярний добуток двох векторів можливий за умови, що вони мають однакову довжину.

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

Якщо записати розміри як кортежі (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 іноді називають внутрішнім добутком.

(Менш поширений) u * v' відповідно називають зовнішнім добутком. Він повʼязаний із тензорним добутком, хоч це вже далеко за межами наших розглядів.

Обертання

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

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

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

Обертання в 3D потребують складнішої матриці з двома кутами. Це важко показати в документі Markdown без великої кількості LaTeX, тож формулу можна знайти у Вікіпедії.

У 3D-графіці, як-от OpenGL та його наступники, прийнято працювати з однорідними координатами.

Кожну вершину подають як 4-вектор (або, що рівнозначно, стовпець матриці): [x, y, z, 1.0] для точки в (x, y, z).

Це дає змогу будувати складніші матриці перетворень. Обертання й далі міститься в A[1:3, 1:3], зсуви - у A[1:3, 4] як [Δx, Δy, Δz], масштабування - на діагоналі, а для скосу, перспективи тощо є й інші можливості.

Швидкодія

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

Швидкодія - велика тема, у яку ми не зможемо заглибитися, але для ілюстрації спробуймо перемножити випадкові матриці різних розмірів.

Наведений нижче код - це швидка й приблизна оцінка (для точнішого підходу скористайтеся BenchmarkTools.jl).

Використано невеликий компʼютер за 430 доларів (США): процесор Ryzen 9, 32 ГБ оперативної памʼяті, Linux Mint 22.1, Julia 1.11.6.

julia> mmul(A) = A * A
mmul (generic function with 1 method)
julia> A = rand(Float64, 1_000, 1_000);
julia> @time mmul(A);
  0.011055 seconds (3 allocations: 7.629 MiB)
julia> A = rand(Float64, 10_000, 10_000);
julia> @time mmul(A);
  7.236284 seconds (1.61 k allocations: 763.027 MiB, 17.92% gc time, 0.15% compilation time)

Приблизно: пара матриць з мільйоном елементів множилася кілька мілісекунд, а матриці зі 100 мільйонами елементів - кілька секунд. Щоб додати ще кілька порядків, знадобиться краще обладнання...