Редактирование: Разреженные матрицы

{{Понятие
|Description=Способы хранения разреженных матриц
Варианты представления матриц
Базовые операции с разреженными матрицами
|Field_of_knowledge=Информатика, Математика
|Возрастная категория=16
}}
Способы хранения разреженных матриц
	Для хранения разреженных матриц целесообразно хранить только ненулевые элементы, чтобы обеспечить экономию памяти и числа операций, необходимых для преобразования матрицы в процессе решения линейных систем, а также простоту доступа к любому элементу матрицы при формировании системы.
Координатный формат
	Используется 3 массива: AA, JR, JC. В первом массиве хранятся значения ненулевых элементов. В двух последних хранятся их координаты внутри матрицы (индексы строк и столбцов соответственно). Данный способ хранения является полным, поскольку представлена вся матрица, и неупорядоченным, поскольку элементы матрицы могут храниться в произвольном порядке.
	Данный формат обеспечивает медленный доступ к элементам матрицы, и является затратным по используемой памяти. Для хранения данных используется 3a чисел (a – кол-во ненулевых элементов).
 
Сжатый разреженный строчный формат
	Если в координатном формате элементы были бы упорядочены построчно, то массив JC содержал бы излишнюю информацию. В таком случае он может быть заменен на массив, который указывает на начало каждой строки, что обеспечит экономию памяти.
	Новая структура данных содержит 3 массива, выполняющих следующие функции:
	AA содержит ненулевые значения матрицы, записанные построчно с 1 до n строки
	JA содержит индексы столбцов всех элементов, представленных в матрице
	IA содержит указатели на начало каждой строки в массивах AA и JA. Таким образом, содержимое IA(i) – позиция в массивах AA и JA, где начинается i-я строка. Длина IA равна n+1, IA(n+1) равен количеству ненулевых элементов + 1.
 
	Данная схема хранения предпочтительнее координатного формата, поскольку часто оказывается удобнее для выполнения типовых вычислений.
	Существуют вариации сжатого разреженного строчного формата. Наиболее очевидный вариант – хранение столбцов вместо строк. Данная схема известна как Сжатый разреженный столбцовый формат.
Модифицированный разреженный строчный формат
	Другой вариант основывается на факте, что диагональные элементы большинства матриц обычно ненулевые и/или что они требуются чаще, чем остальные элементы. Как результат, они могут храниться отдельно. Таким образом, используется два массива:
	AA, в котором первые n позиций содержат диагональные элементы матрицы по порядку. Неиспользованная позиция n+1 может иногда содержать информацию о матрице. Начиная с позиции n+2, идут ненулевые элементы матрицы, исключая диагональные элементы, упорядоченные построчно.
	JA. JA(k) представляет значение индекса столбца для каждого элемента AA(k). n+1 первых позиций JA содержат указатель на начало каждой строки в AA и JA. Начиная с n+2, JA указывает на индексы столбцов элементов с соответствующими индексами в AA.
 
Диагонально структурированный формат
	Диагонально структурированные матрицы – это матрицы, ненулевые элементы которой расположены вдоль небольшого числа диагоналей. Эти диагонали могут храниться в двумерном массиве DIAG(n, Nd), где Nd – количество диагоналей. Смещения для каждой из диагональной относительно главной диагонали должны быть известны, они хранятся в массиве IOFF(Nd). Таким образом, элемент a_(i,i+ioff(j)) изначальной матрицы расположен в позиции (i, j) массива DIAG.
	Все диагонали, кроме главной, будут иметь меньшее, чем n, число элементов, поэтому в DIAG будут неиспользуемые позиции.
 
 
Ellpack-Itpack формат
	Для такого формата предполагается, что в диагоналях матрицы большинство элементов в каждой строке ненулевые. В данном случае используется два массива. Первый массив COEF похож на массив DIAG из диагонально структурированного формата и содержит ненулевые элементы матрицы A. Ненулевые элементы каждой строки матрицы могут храниться в массиве COEF(n, Nd), заполняя некоторые строки нулями, если необходимо. Целочисленный массив JCOEF(n, Nd) должен быть заполнен и содержать индексы столбцов для каждого элемента в COEF.
 
	Определенное число столбцов должно быть выбрано для каждого нулевого элемента, который должен быть добавлен, чтобы дополнить сокращенные строки в матрице A (в примере: 1, 4, 5 строки).
Случай с симметричными матрицами
	Когда обрабатываемая матрица симметрична, то есть A=A^T, достаточно хранить лишь ее верхнюю треугольную подматрицу. При этом для хранения можно использовать любую из рассмотренных схем.
	Если симметричная матрица является положительно определенной, то все ее диагональные элементы отличны от нуля, и могут храниться в отдельном массиве, а некоторым разреженным форматом представления будет представления только верхний треугольник матрицы.
Варианты представления матриц
Блочное представление
	Блочная матрица – матрица, которая разбивается на сектора, называемые блоками. Визуально это можно отобразить, проведя множество горизонтальных и вертикальных линий в изначальной матрице, которые будут ее разбивать на отдельные блоки матриц меньшего размера.
	Матрица:
 
	Может быть разделена на 4 блока с матрицами 2×2:
 
	После чего разбитая на блоки матрица может быть записана в следующем виде:
 
	Операции с блочными матрицами производятся по тем же правилам, как если бы на месте блоков были числовые элементы. Для выполнимости операций необходимо соответствующее согласование размеров блоков. Например, при умножении блочных матриц требуется, чтобы горизонтальные размеры блоков первого сомножителя совпадали с соответствующими вертикальными размерами второго сомножителя.
Ленточное представление
	Ленточная матрица – разреженная матрица, все ненулевые элементы которой находятся внутри диагональной ленты, расположенной вдоль главной диагонали матрицы.
	Формально ленточное представление матрицы описывается следующим образом:
 
	Где k_1,k_2 – нижняя и верхняя ширина ленты соответственно. Общая ширина ленты есть максимальное k_i.
	Для хранения ленточных матриц обычно сохраняются диагонали внутри ленты, оставшиеся элементы неглавной диагонали заполняются нулями.
	Например, есть ленточная матрица 6×6:
 
	Она может храниться как матрица 6×3:
 
Профильное представление
	Профилем матрицы называют сумму всех локальных ширин относительно главной диагонали. Совокупность позиций от локальной ширины до главной диагонали называется оболочкой. Формально это описывается следующим образом:
f_i (A)=min⁡{ j ┤| a_ij≠0 }  
β_i (A)=i-f_i (A)
	Число f_i (A) – столбцовый индекс первого ненулевого элемента i-й строки, β_i (A) – ширина ленты в i-й строке.
	Тогда профиль матрицы:
e(A)=∑_(i=1)^n▒〖β_i (A) 〗
	Например, для матрицы, представленной ниже, e(A)=11:
 

Базовые операции с разреженными матрицами
Умножение матрицы на вектор
	Данная операция является важной и требуется во многих итерационных методах для решения разреженных линейных систем. В результате получается вектор. Умножение происходит циклично путем умножения каждой строки матрицы на столбец вектора. 
	При использовании сжатого разреженного строчного формата для хранения матрицы каждая итерация цикла вычисляет отдельный компонент результирующего вектора и это является плюсом, поскольку данные компоненты могут вычисляться независимо друг от друга.
 
	При использовании сжатого разреженного столбцового формата в каждой итерации цикла произведение j-го столбца матрицы на вектор прибавляется к соответствующему компоненту результирующего вектора. В данном случае итерации не могут выполняться независимо, поскольку внутри нескольких итераций могут использоваться значения одной и той же составляющей результирующего вектора.
 
	Для использования параллельного вычисления при данном формате возможно независимое вычисление компонентов результирующего вектора во внутреннем цикле.
Умножение матрицы на матрицу
	Для умножения матрицы на матрицу P=AB необходимо, чтобы число столбцов матрицы A было равно числу строк матрицы B. Тогда результатом операции умножения будет матрица P. 
	Пусть A имеет размерность (M×N), B-(N×R). Тогда размерность матрицы P, которая является результатом умножения AB, будет (M×R). Элементы матрицы P вычисляются следующим образом:
P_ij=∑_(k=1)^N▒〖A_ik B_kj 〗,i=1,2,…,M,j=1,2,…,R.
	Компоненты матриц A и B не изменяются и только считываются, поэтому вычисления элементов P могут выполняться независимо друг от друга, что позволяет использовать параллельное вычисление для ускорения работы.
Решение треугольной системы
	Решение верхней или нижней треугольной системы является важным аспектом при работе с разреженными матрицами. Происходят итерации, на каждой из которых вычисляется значение компонента вектора. Вычисление происходит путем вычитания из соответствующего вектора правых частей вычисленных ранее компонентов вектора неизвестных, умноженных на соответствующие коэффициенты из разреженной матрицы, после чего происходит деление результата на коэффициент при вычисляемом компоненте вектора неизвестных.
Типовые задачи
МКР
	Пусть есть исследуемая двумерная область, состояние которой описывается уравнением:
 
	Где Ω – прямоугольная область размером l_1×l_2, а Г – ее границы. Обе стороны области могут быть дискретизированы путем нанесения n_1+2 и n_2+2 точек в соответствующих направлениях x_1,x_2:
 
	Тогда нанесенная на область сетка будет выглядеть следующим образом:
 
	Где
 
	Поскольку значения на границах известны, то считаются только внутренние точки, которые пронумерованы. Для пронумерованных узлов сетки записываются разностные уравнения по формуле центральной разности, аппроксимирующие состояния в узлах:
 
	После записи аппроксимирующих уравнений для всех узлов получается система вида Ax=b. Матрица A имеет следующий вид:
 
	Проведенные горизонтальные и вертикальные линии, разделяющие узлы матрицы, показывают, что матрица имеет блочную структуру следующего вида (для случая, когда h_1=h_2):
A=1/h^2  (■(B&-I&@-I&B&-I@&-I&B)),B=(■(■(4&-1@-1&4)&■(&&@-1&&)@■(&-1@&@&)&■(4&-1&@-1&4&-1@&-1&4)))
LU разложение
	При решении системы линейных уравнений вида Ax=b матрица A может быть разложена в произведение матриц с использованием некоторого алгоритма. Универсальным разложением квадратной матрицы является A=LU, где L- нижняя треугольная матрица, U- верхняя треугольная матрица. Система представляется в виде: LUx=b. Производится замена Ux=y, после чего находится x с помощью решения двух треугольных систем:
{█(Ly=b@Ux=y)┤
	Найти матрицы L и U можно следующим образом:
	Сперва элементы матриц инициализируются таким образом, что диагональные элементы матрицы L равны 1, все остальные – 0:
u_ij=0,l_ij=0,l_ii=1;     i=1,…,n,      j=1,…,n
	Далее циклично вычисляются элементы матриц по формулам:
if (i≤j):     u_ij=a_ij- ∑_(k=1)^(i-1)▒〖l_ik*u_kj 〗
if (i>j):     l_ij=(a_ij- ∑_(k=1)^(j-1)▒〖l_ik*u_kj 〗)/u_jj
Разложение Холецкого
Метод Холецкого
Разложение Холецкого – представление симметричной положительно-определенной матрицы в виде A=LL^T, где L – нижняя треугольная матрица.
	Метод Холецкого – это симметричный вариант гауссова исключения, основной численный алгоритм, созданный и используемый для решения симметричных положительно определенных матричных задач. Для положительно определенных матриц разложение всегда существует.
	Пусть есть система уравнений:
Ax=b
	Где A=NxN – симметричная положительно определенная матрица коэффициентов. b – вектор размера N, называемый правой частью, x – вектор-решение размера N, который нужно вычислить. Применение метода Холецкого приводит матрицу А к треугольному разложению:
A=LL^T
	Данное разложение всегда существует, когда матрица А симметрична и положительно определена.
LL^T x=b
	Замена y=L^T x показывает, что x можно получить, решая треугольные системы
{█(Ly=b@L^T x=y)┤
	При использовании метода Холецкого для разреженной матрицы A множитель Холецкого L обычно имеет ненулевые элементы в позициях, где в A стояли нули. 
Алгоритм разложения
	Симметричная матрица A – положительно определенная, если x^T Ax>0 для всех ненулевых векторов x. Диагональные элементы всегда положительны. 
	Для такой матрицы существует единственное ее треугольное разложение LL^T, где L – нижняя треугольная матрица. С положительными диагональными элементами.
Форма внешних произведений
	Один из алгоритмов разложения происходит в форме внешних произведений. Схема описывается на матричном языке рекурентно:
A=A_0=H_0=[■(d_1&v_1^T@v_1&¯(H_1 ))]=[■(√(d_1 )&0@v_1/√(d_1 )&I_(N-1) )][■(1&0@0&¯(H_1 )-(v_1 v_1^T)/d_1 )][■(√(d_1 )&v_1/√(d_1 )@0&I_(N-1) )]=L_1 [■(1&0@0&H_1 )] L_1^T=L_1 A_1 L_1^T
A_1=[■(1&0@0&H_1 )]=[■(1&0&0@0&d_2&v_2^T@0&v_2&¯(H_2 ))]=[■(1&0&0@0&√(d_2 )&0@0&v_2/√(d_2 )&I_(N-2) )][■(1&0&0@0&1&0@0&0&¯(H_2 )-(v_2 v_2^T)/d_2 )][■(1&0&0@0&√(d_2 )&v_2/√(d_2 )@0&0&I_(N-2) )]=L_2 A_2 L_2^T
…
A_(N-1)=L_N I_N L_N^T
	Для значений i от 1 до N: 
	d_i – положительный скаляр
	v_i – вектор длины N-i
	H_i – положительно определенная симметричная матрица порядка N-i
	По завершении N итераций матрица A становится представима в виде:
A=L_1 L_2…L_N I_N L_N^T…L_2^T L_1^T=LL^T
	В этой схеме вычисляются один за другим столбцы L. В то же время каждый шаг требует модификации подматрицы ¯(H_i ), посредством внешнего произведения (v_i v_i^T)/d_i . Результатом является H_i, то есть как раз та матрица, которую остается разложить. Обращение к элементам происходит следующим образом:
 
	Метод окаймления
	Процесс разложения может быть сформулирован иначе – с помощью метода окаймления.
	Пусть матрица представлена в виде:
A=[■(M&u@u^T&s)]
	Причем уже получено симметричное разложение L_M L_M^T ведущей главной подматрицы M порядка N-1. В таком случае разложением A будет:
A=[■(L_M&0@w^T&t)][■(L_M&w@0&t)]
	Где 
w=L_M^(-1) u,t=√(s-w^T w)
	В схеме поочередно вычисляются строки L. К еще неразложенной части матрицы обращение не производится до тех пор, пока не будет вычислена соответствующая часть L. Последовательность вычислений:
 
Форма скалярных произведений
	Также разложение осуществимо в форме скалярных произведений, когда элементы матрицы L вычисляются, начиная с верхнего левого угла матрицы, сверху вниз, слева направо, по формулам:
l_11=√(a_11 )
l_j1=a_j1/l_11 ,j∈[2,n]
l_ii=√(a_ii-∑_(p=1)^(i-1)▒l_ip^2 ),i∈[2,n]
l_ji=1/l_ii  (a_ji-∑_(p=1)^(i-1)▒〖l_ip l_jp 〗),i∈[2,n-1],j∈[i+1,n]
	Последовательность вычислений выглядит следующим образом:
 
	Количество операций при использовании данного разложения для матрицы размерности N:
	N- квадратных корней
	(N-1)N/2- делений
	(N^2 (N-1))/6- умножений
	(N^2 (N-1))/6- сложений/вычитаний

LDL^T разложение
	Альтернативная форма разложения Холецкого, исключающая необходимость вычисления квадратных корней, выглядит следующим образом:
A=LDL^T
	Где L- нижняя треугольная матрица, D – диагональная матрица.
	Вычисление элементов происходит последовательно слева-направо и сверху-вниз:
d_1=a_11
l_j1=a_j1/d_1 ,j∈[2,n]
d_i=a_ii-∑_(k=1)^(i-1)▒〖l_ik^2 d_k 〗,i∈[2,n]
l_ji=1/d_i  (a_ji-∑_(k=1)^(i-1)▒〖l_jk l_ik d_k 〗),i∈[2,n-1],j∈[i+1,n] 
	Количество операций при использовании данного разложений для матрицы размерности N:
	(N(N-1))/2- делений
	(N-1)N(N+1)/3- умножений
	((N-1)N(N+1))/6- сложений/вычитаний
Сравнение разложений
	При LU разложении элементы матрицы A могут быть заменены элементами матриц L и U, при этом отсутствует необходимость выделения дополнительной памяти:
 
	В случае, когда матрица A симметрична и положительно определена, то есть A_ij=A_ji (i,j=1,…,n) и P^T AP>0  (∀P≠0), целесообразно использовать разложение Холецкого вида A=LL^T, поскольку оно выгоднее LU разложения, ведь необходимо найти и хранить только одну треугольную матрицу.
	К примеру, в задачах, где применяются МКР или МКЭ, матрица жесткости A является симметричной и положительно определенной, поэтому в данных случаях стоит применять именно разложение Холецкого в целях экономии.
Перестановки
	Применение перестановок дает возможность значительно уменьшить заполнение матрицы после разложения Холецкого, снижая необходимый объем памяти для хранения.
	Важный момент состоит в том, что треугольный множитель L разрежен в точности в той мере, что и нижний треугольник исходной матрицы A. 
	Например:
  
 
	Упорядочение можно определить до начала разложения, поэтому можно также заранее определить и местоположение заполнения, которое произойдет при разложении, поэтому способ хранения L можно выбрать до реального разложения и зарезервировать место для элементов заполнения. Дальнейшие вычисления проводятся при структуре хранения, остающейся статичной (неизменной). Таким образом, три задачи могут быть разделены как самостоятельные объекты исследования и как разные модули.
	Во многих инженерно-конструкторских приложениях приходится решать большое количество различных задач с положительно определенными матрицами, имеющих одинаковую структуру. Упорядочение и формирование схемы хранения нужно выполнить лишь однажды, поэтому желательно, чтобы эти этапы были изолированы от реальных вычислений.
Профильное упорядочение
	Обратный алгоритм Катхилла-Макки широко используется для упорядочения элементов матриц, уменьшает объем профиля матрицы перед применением разложения.
	Схему Катхилла-Макки можно рассматривать как метод уменьшения ширины ленты матрицы. Упорядочение, получаемое обращением упорядочения Катхилла-Макки, часто гораздо сильнее уменьшает профиль, чем первоначальное упорядочение, хотя ширина ленты, остается неизменной. Это упорядочение называется обратным упорядочением Катхилла-Макки. 
	Алгоритм описывается следующими шагами:
	Определяется начальный узел r и выполняется присваивание r→x_1
	(Основной цикл). Для i=1,…,N найти всех ненумерованных соседей узла x_i и занумеровать их в порядке возрастания степеней.
	(Обратное упорядочение). Обратное упорядочение Катхилла-Макки есть y_1,y_2,…,y_N, где y_i=x_(N-i+1) для i=1,…,N.
	По завершении алгоритма объем профиля матрицы снижается, что снижает требуемое для разложения количество памяти.