Лінійна алгебра для симуляцій: розріджені розв'язувачі, розклади та власні значення

1. Чому лінійні системи домінують у симуляції

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

МСЕ: глобальна матриця жорсткості

Метод скінченних елементів для механіки конструкцій дає:

N — кількість ступенів свободи (3 на вузол у 3D). Для сітки зі 100 000 вузлів N ≈ 300 000. K має ширину стрічки, визначену зв'язністю сітки — приблизно 20–30 ненульових на рядок — що дає близько 6 мільйонів ненульових елементів із загальних 90 мільярдів: розрідженість 99.993%.

Рідина: рівняння Пуассона для тиску

У симуляціях нестисливих Нав'є–Стокса (напр. рідина FLIP, SPH), забезпечення ∇·u = 0 вимагає розв'язування:

Дискретний лапласіан на сітці 256³ має N ≈ 16.7 мільйона невідомих. Пряма факторизація вимагала б O(N³) ≈ 4.6 × 10²¹ операцій — абсолютно нездійсненно. Саме тому існують ітераційні розв'язувачі.

Спільні властивості

• Розрідженість: O(N) ненульових із N² можливих елементів
• Симетрія: Для неорієнтованих графів або самоспряжених операторів A = Aᵀ
• Додатна визначеність (SPD): xᵀAx > 0 для всіх x ≠ 0 — уможливлює спеціальні алгоритми
• Структурована розрідженість: матриці МСЕ часто мають блочну структуру, що відповідає топології сітки

2. Формати розріджених матриць: CSR, CSC, COO

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

COO — координатний формат

Три масиви: індекси рядків, індекси стовпців, значення. Невпорядкований, допускає дублікати елементів (row, col) (підсумовуються при перетворенні). Найкращий для складання матриці.

CSR — стиснений за рядками

Найпоширеніший формат для арифметики. Замінює масив рядків на рядковий вказівник (довжиною N+1, де rowptr[i+1] − rowptr[i] — кількість ненульових у рядку i):

CSC — стиснений за стовпцями

Ролі рядків та стовпців помінялися. Переважний за угодою LAPACK та для стовпцево-орієнтованих операцій (напр. вибору головного елемента при факторизації). Власний формат MATLAB — CSC.

Порівняння форматів

3. Прямі розв'язувачі: LU, Холецький, QR

Прямі методи факторизують A один раз, а потім розв'язують через трикутну зворотну підстановку за O(N) операцій — але сама факторизація може бути дорогою, а заповнення (поява ненульових там, де A була нулем) обмежує масштабованість до N ≲ 10⁵ у 3D.

LU-розклад

Для загальної матриці n×n: A = P · L · U, де P — матриця перестановки (частковий вибір головного елемента), L — нижня трикутна з одиничною діагоналлю, U — верхня трикутна.

Розріджений LU через UMFPACK або SuperLU використовує порядки, що зменшують заповнення. Порядок вкладеного розсікання дає O(N^1.5) заповнення для 2D- сіток та O(N^2) для 3D — набагато краще за природний порядок, але все ще гірше за ітераційні розв'язувачі для дуже великих 3D-задач.

Розклад Холецького (SPD)

Для симетричних додатно визначених матриць: A = LLᵀ, де L — нижня трикутна з додатною діагоналлю. Вибір головного елемента не потрібен (стійкість гарантована для SPD). Приблизно вдвічі швидше за LU — зберігається й обробляється лише нижній трикутник.

Застосування: матриці жорсткості МСЕ, коваріаційні матриці в імовірнісних методах, множники SPD-передобумовлювачів (неповний Холецький IC(0)).

QR-розклад

Факторизує A = QR, де Q — ортогональна (Qᵀ = Q⁻¹), а R — верхня трикутна. Випадки застосування:

• Метод найменших квадратів: min‖Ax − b‖₂ → Rx = Qᵀb (n стовпців, m рядків, m > n)
• Власні значення: QR-ітерація для щільних матриць
• Ортогоналізація Грама-Шмідта: основа для методів підпросторів Крилова

Обчислюється через відбивачі Хаусхолдера або повороти Гівенса. QR Хаусхолдера числово кращий за класичний Грам-Шмідт — похибка ортогональності становить O(ε·κ(A)) проти O(ε·κ(A)²) для CGS.

4. Ітераційні розв'язувачі: CG, GMRES, BiCGSTAB

Ітераційні методи починають з початкового наближення x₀ й покращують його по одному «кроку» за раз. Кожен крок коштує O(nnz) для розрідженого множення матриці на вектор. Вони ніколи повністю не факторизують A — що робить їх застосовними до задач, де прямі методи нездійсненні.

Спряжені градієнти (CG) — лише SPD

Метод CG (Hestenes & Stiefel, 1952) мінімізує квадратичну форму f(x) = ½ xᵀAx − bᵀx за послідовністю A-спряжених напрямків пошуку. Для SPD-матриці він збігається щонайбільше за N кроків у точній арифметиці (N = розмір матриці), але на практиці збігається до машинної точності значно раніше.

GMRES — загальні матриці

Метод узагальненого мінімального залишку (Saad & Schultz, 1986) обробляє несиметричні та невизначені системи. Він будує ортонормований базис Крилова {v₁, v₂, ..., vₖ} через ітерацію Арнольді й мінімізує ‖b − Axₖ‖₂ над підпростором Крилова Kₖ = span{v₁, ..., vₖ}.

GMRES(k) з перезапуском перезапускає базис Крилова кожні k кроків, щоб контролювати пам'ять. Компроміс: менше k використовує менше пам'яті, але може стагнувати на складних системах. Гнучкий GMRES (FGMRES) дозволяє передобумовлювачу змінюватися між ітераціями.

BiCGSTAB

Стабілізований біспряжений градієнт (van der Vorst, 1992) — це варіант для несиметричних систем, що уникає зростання пам'яті GMRES. Він використовує дві послідовності залишків для досягнення плавнішої збіжності, ніж BiCG (який може поводитися хаотично). Пам'ять: O(N) — лише фіксована кількість векторів незалежно від кількості ітерацій. Переважний над GMRES, коли пам'яті обмаль.

5. Передобумовлювачі: Якобі, ILU, багатосіткові

Передобумовлювач M апроксимує A й перетворює систему на M⁻¹Ax = M⁻¹b. Збіжність CG/GMRES залежить від числа обумовленості κ(M⁻¹A) — хороший передобумовлювач зменшує κ за невелику ціну.

Передобумовлювач Якобі (діагональний)

Неповний LU: ILU(0) та ILU(k)

Неповний LU факторизує A ≈ L̃Ũ, де L̃ та Ũ мають той самий шаблон розрідженості, що й A (ILU(0)), або допускають k рівнів додаткового заповнення (ILU(k)). Наближені множники дешеві в застосуванні (розріджені трикутні розв'язування) та суттєво зменшують κ.

ILU(k) допускає заповнення до відстані k від наявних ненульових. ILU(1) чи ILU(2) часто зменшує кількість ітерацій у 5–10× порівняно з Якобі за ціну в 3–5× більше пам'яті на множник.

Багатосіткові методи: оптимальна збіжність O(N)

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

• Згладжування: кілька ітерацій Гаусса-Зейделя чи Якобі гасять високочастотні компоненти похибки
• Звуження: проєкція залишку на грубішу сітку (повне зважування або ін'єкція)
• Грубий розв'язок: точний розв'язок на найгрубішому рівні (тривіально малий)
• Продовження: інтерполяція корекції назад на дрібну сітку
• V-цикл: один прохід униз + один прохід угору; W-цикл: два рекурсивні проходи вниз

Геометричні багатосіткові методи вимагають структурованої ієрархії сіток. Алгебраїчні багатосіткові методи (AMG) будують ієрархію автоматично з матриці — BoomerAMG з HYPRE та ML/MueLu (Trilinos) є стандартними бібліотеками в науковому обчисленні. Для рівнянь Пуассона на однорідних сітках багатосіткові методи збігаються за O(1) ітерацій незалежно від розміру сітки.

6. Алгоритми власних значень: степеневий, QR, Ланцош

Власні значення з'являються всюди, де симуляція має власні частоти, аналіз стійкості чи модальний розклад. Знаходження всіх власних значень матриці N×N коштує O(N³) — здійсненно лише для N ≲ 10³. Великі розріджені задачі на власні значення вимагають методів, що дешево обчислюють кілька власних значень.

Степеневі ітерації

Знаходять власне значення з найбільшим модулем λ_max та його власний вектор:

Обернені ітерації (замінити A на (A − σI)⁻¹) знаходять власне значення, найближче до σ, використовуючи одне лінійне розв'язання за крок — ідеально, коли прямий розв'язувач дешевий, а цільове власне значення відоме приблизно. Зсув-обернення з CG/GMRES замінює щільне розв'язання розрідженою ітерацією.

Щільний QR-алгоритм

Для щільних матриць (N ≲ 10³) усі власні значення можна знайти через QR- ітерацію:

Алгоритм Ланцоша — великі розріджені

Алгоритм Ланцоша (1950) будує ортонормований базис Крилова й зводить A до симетричної тридіагональної матриці T розміру k×k (k ≪ N). Власні значення T (значення Рітца) апроксимують крайні власні значення A. Потрібно лише k SpMV — ідеально для знаходження 10–100 найбільших чи найменших власних значень розрідженої A.

На практиці накопичується втрата ортогональності — що вимагає або повної реортогоналізації (підхід IRLM з ARPACK), або вибіркової ортогоналізації (частковий Ланцош). ARPACK (Arnoldi Package) — стандартна бібліотека: її прив'язка для Python — це `scipy.sparse.linalg.eigsh` та `eigs` зі SciPy.

7. SVD, PCA та модальний розклад

Сингулярний розклад матриці A розміру m×n:

Теорема Еккарта-Янга стверджує, що найкраще наближення A рангу k (у будь-якій ортогонально інваріантній нормі) — це:

Власний ортогональний розклад (POD)

У гідродинаміці та динаміці конструкцій POD (SVD матриці знімків) виокремлює домінантні просторові моди симуляції. Якщо U = [u(t₁) | u(t₂) | ... | u(tₙ)] містить N знімків швидкості, то:

• Стовпці U (ліві сингулярні вектори) — це просторові моди POD, упорядковані за енергетичним вмістом σᵢ²
• k мод охоплюють частку (σ₁² + ... + σₖ²) / ‖U‖_F² від загальної дисперсії
• Проєкція РЧП на r мод POD дає модель зниженого порядку (ROM) розміру r ≪ N

Динамічний модальний розклад (DMD)

DMD (Schmid, 2010) знаходить найкраще наближений лінійний оператор, що відображає знімок Uₖ на Uₖ₊₁:

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

8. Розріджений метод спряжених градієнтів на JavaScript

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

// ============================================================
// Розріджена CSR-матриця + передобумовлені спряжені градієнти (PCG)
// Розв'язує A·x = b, де A — SPD і розріджена (формат CSR)
// ============================================================

class SparseCSR {
  constructor(n) {
    this.n = n;
    this.rowptr = new Int32Array(n + 1);
    this.cols   = [];    // буде перетворено на Int32Array
    this.vals   = [];    // Float64Array після побудови
  }

  // Складання з COO-трійок [{row, col, val}]
  static fromCOO(n, coo) {
    // Сортуємо за рядком, потім за стовпцем
    coo.sort((a, b) => a.row !== b.row ? a.row - b.row : a.col - b.col);
    const mat   = new SparseCSR(n);
    const cols  = [], vals = [];
    const rowptr = new Int32Array(n + 1);

    for (const {row, col, val} of coo) {
      rowptr[row + 1]++;
      cols.push(col);
      vals.push(val);
    }
    for (let i = 1; i <= n; i++) rowptr[i] += rowptr[i - 1];
    mat.rowptr = rowptr;
    mat.cols   = new Int32Array(cols);
    mat.vals   = new Float64Array(vals);
    return mat;
  }

  // y = A·x (розріджений добуток матриці на вектор)
  matvec(x, y) {
    const {n, rowptr, cols, vals} = this;
    for (let i = 0; i < n; i++) {
      let s = 0;
      for (let j = rowptr[i]; j < rowptr[i + 1]; j++) s += vals[j] * x[cols[j]];
      y[i] = s;
    }
  }

  // Діагональний передобумовлювач (Якобі)
  diagInv() {
    const d = new Float64Array(this.n);
    const {n, rowptr, cols, vals} = this;
    for (let i = 0; i < n; i++) {
      for (let j = rowptr[i]; j < rowptr[i + 1]; j++) {
        if (cols[j] === i) { d[i] = 1 / vals[j]; break; }
      }
    }
    return d;
  }
}

// Розв'язує A·x = b за допомогою передобумовленого CG
// Повертає {x, iter, residual}
function pcg(A, b, opts = {}) {
  const {tol = 1e-6, maxIter = 1000} = opts;
  const n = A.n;
  const diag  = A.diagInv();          // передобумовлювач Якобі M⁻¹
  const x     = new Float64Array(n);
  const r     = new Float64Array(n);
  const z     = new Float64Array(n);
  const p     = new Float64Array(n);
  const Ap    = new Float64Array(n);

  // r₀ = b − A·x₀  (x₀ = 0 → r₀ = b)
  for (let i = 0; i < n; i++) r[i] = b[i];

  // z₀ = M⁻¹·r₀  (Якобі: покомпонентне множення)
  let rz = 0;
  for (let i = 0; i < n; i++) { z[i] = r[i] * diag[i]; p[i] = z[i]; rz += r[i] * z[i]; }

  let iter = 0;
  while (iter < maxIter) {
    A.matvec(p, Ap);
    let pAp = 0;
    for (let i = 0; i < n; i++) pAp += p[i] * Ap[i];

    const alpha = rz / pAp;
    let rnorm2 = 0;
    for (let i = 0; i < n; i++) {
      x[i] += alpha * p[i];
      r[i] -= alpha * Ap[i];
      rnorm2 += r[i] * r[i];
    }

    if (Math.sqrt(rnorm2) < tol) break;

    let rzNew = 0;
    for (let i = 0; i < n; i++) { z[i] = r[i] * diag[i]; rzNew += r[i] * z[i]; }
    const beta = rzNew / rz;
    rz = rzNew;
    for (let i = 0; i < n; i++) p[i] = z[i] + beta * p[i];
    iter++;
  }
  return {x, iter, residual: Math.sqrt(rz)};
}

// ---- Приклад: 1D-Пуассон на однорідній сітці, 5 внутрішніх точок ----
// −u'' = f,  u(0)=u(1)=0,  f=1   → u_i = xi(1−xi)/2
const N  = 5;
const h  = 1 / (N + 1);
const h2 = h * h;

// Будуємо тридіагональну матрицю лапласіана у форматі COO
const coo = [];
for (let i = 0; i < N; i++) {
  coo.push({row: i, col: i, val:  2 / h2});
  if (i > 0)     coo.push({row: i, col: i-1, val: -1 / h2});
  if (i < N - 1) coo.push({row: i, col: i+1, val: -1 / h2});
}
const A = SparseCSR.fromCOO(N, coo);
const b = new Float64Array(N).fill(1);         // f=1

const {x, iter} = pcg(A, b, {tol: 1e-10});

console.log("PCG збігся за", iter, "ітерацій");
x.forEach((u, i) => {
  const xi  = (i + 1) * h;
  const ref = xi * (1 - xi) / 2;
  console.log(`u[${i}] = ${u.toFixed(8)}  (ref: ${ref.toFixed(8)})`);
});