8.3: Використання розріджених матриць

При роботі з розрідженими матрицями розміру\(10^3\times 10^3\) або більшого розміру слід використовувати розріджені формати матриць замість звичайних 2D масивів. (Для матриць розміром\(100\times 100\) або менше відмінності у продуктивності та використанні пам'яті зазвичай незначні.)

Зазвичай добре вибирати або формат CSR, або CSC, в залежності від того, які математичні операції ви маєте намір виконувати. Ви можете побудувати матрицю або (i) безпосередньо, за допомогою кортежу (data, idx), як описано вище, або (ii) шляхом створення матриці LIL, заповнення потрібних елементів, а потім перетворення її у формат CSR або CSC.

Якщо ви маєте справу з розрідженою матрицею, яка є «сильно діагональною» (тобто ненульові елементи займають дуже малу кількість діагоналей, таких як тридіагональна матриця), то можна розглянути можливість використання формату DIA. Основна перевага формату DIA полягає в тому, що його дуже легко побудувати, подаючи (дані, зсуви) вхід в функцію dia_matrix, як описано вище. Однак формат зазвичай не працює значно краще CSR/CSC; а якщо матриця не сильно діагональна, її продуктивність набагато гірше.

Іншим поширеним способом побудови розрідженої матриці є використання функцій scipy.sparse.diags або scipy.sparse.spdiags. Ці функції дозволяють вказати вміст матриці з точки зору її діагоналей, а також який розріджений формат використовувати. Ці дві функції мають дещо різні умови виклику; докладніше див. документацію.

8.3.1 Точковий метод

Кожна розріджена матриця має точковий метод, який обчислює добуток матриці з входом (у вигляді 1D або 2D масиву, або розрідженої матриці), і повертає результат. Для розріджених матричних добутків цей метод слід використовувати замість автономної функції, аналогічно названої точкою, яка використовується для нерозріджених матричних добутків. Метод розрідженої матриці використовує розрідженість матриці для прискорення обчислення виробу. Зазвичай формат CSR працює швидше, ніж інші розріджені формати.

Наприклад, припустимо, ми хочемо розрахувати товар\(\mathbf{A} \vec{x}\), де

\[\mathbf{A} = \begin{bmatrix}0 & 1.0 & 0 & 0 & 0 \\ 0 & 2.0 & -1.0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 6.6 & 0 & 0 & 0 & 1.4\end{bmatrix}, \quad \vec{x} = \begin{bmatrix}1\\1\\2\\3 \\5\end{bmatrix}.\]

Це може бути досягнуто за допомогою наступної програми:

from scipy import *
import scipy.sparse as sp

data = [1.0, 2.0, -1.0, 6.6, 1.4]
rows = [0, 1, 1, 3, 3]
cols = [1, 1, 2, 0, 4]
A = sp.csr_matrix((data, [rows, cols]), shape=(4,5))
x = array([1.,1.,2.,3.,5.])

y = A.dot(x)

print(y)

Запуск цієї програми дає очікуваний результат:

[  1.    0.    0.   13.6]

8.3.2 вирішити

Функція spolve, надана в модулі scipy.sparse.linalg, є вирішувачем для розрідженої лінійної системи рівнянь. Він приймає два входи,\(\mathbf{A}\) і\(\mathbf{b}\), де\(\mathbf{A}\) повинна бути розріджена матриця;\(\mathbf{b}\) може бути або розрідженим, або звичайним 1D або 2D масивом. Він повертає значення,\(\mathbf{x}\) яке вирішує лінійні рівняння\(\mathbf{x}\). Повертається значення може бути як звичайним масивом, так і розрідженою матрицею, залежно від того, чи\(\mathbf{b}\) є розрідженим.

Для розріджених проблем завжди слід використовувати spsolve замість scipy.linalg.solve (звичайний розв'язувач для нерозріджених проблем). Ось приклад програми, що показує spsolve в дії:

from scipy import *
import scipy.sparse as sp
import scipy.sparse.linalg as spl

## Make a sparse matrix A and a vector x
data = [1.0, 1.0, 1.0, 1.0, 1.0]
rows = [0, 1, 1, 3, 2]
cols = [1, 1, 2, 0, 3]
A = sp.csr_matrix((data, [rows, cols]), shape=(4,4))
b = array([1.0, 5.0, 3.0, 4.0])

## Solve Ax = b
x = spl.spsolve(A, b)
print(" x = ", x)

## Verify the solution:
print("Ax = ", A.dot(x))
print(" b = ", b)

Запуск програми дає:

 x =  [ 4.  1.  4.  3.]
Ax =  [ 1.  5.  3.  4.]
 b =  [ 1.  5.  3.  4.]

8.3.3 Цифри

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

Функція eigs, надана в модулі scipy.sparse.linalg, є власним вирішувачем для розріджених матриць. На відміну від власних розв'язувачів, які ми обговорювали раніше, таких як scipy.linalg.eig, функція eigs повертає лише вказану підмножину власних значень та власних векторів.

Функція eigsh аналогічна eigs, за винятком того, що вона спеціалізована для ермітієвих матриць. Обидві функції використовують низькорівневу числову бібліотеку власних розв'язувачів під назвою ARPACK, яка також використовується GNU Octave, Matlab та багатьма іншими числовими інструментами. Ми не будемо обговорювати, як працює алгоритм.

Перший вхід до eigs або eigsh - це матриця, для якої ми хочемо знайти власні значення. Також приймаються кілька інших необов'язкових входів. Ось найбільш часто використовувані з них:

• Необов'язковий параметр з ім'ям k визначає кількість власних значень, які потрібно знайти (типове значення — 6).
• Необов'язковий параметр з назвою M, якщо вказано, вказує праву матрицю для узагальненої задачі на власне значення.
• Необов'язковим параметром з назвою sigma, якщо вказано, має бути числом; це означає, що потрібно знайти k власних значень, найближчих до цього числа.
• Необов'язковий параметр з іменем, який визначає, які власні значення потрібно знайти, використовуючи критерії, відмінні від сигми: 'LM' означає знайти власні значення з найбільшими величинами, 'SM' означає знайти ті з найменшими величинами і т.д. одночасно вказуємо як сигму, так і яку. При знаходженні малих власних значень зазвичай краще використовувати сигму замість якої (див. Обговорення в наступному розділі).
• Необов'язковий параметр з назвою return_eigenvectors, якщо True (за замовчуванням), означає повернути як власні значення, так і власні вектори. Якщо False, функція повертає тільки власні значення.

Повний список входів наведено у повній документації щодо функцій eigs та eigsh.

За замовчуванням eigs і eigsh повертають два значення: 1D масив (який є складним для eigs і реальним для eigsh), що містить знайдені власні значення, і 2D масив, де кожен стовпець є одним з відповідних власних векторів. Якщо для необов'язкового параметра з ім'ям return_eigenvectors встановлено значення False, то повертається тільки 1D масив власних значень.

У наступному розділі ми побачимо приклад використання eigsh.