Допустимо свободное использование материалов данного раздела, в том числе и для создания законченных коммерческих программ с указанием источника. http://drobotenko.com/ . Для написания платных книг, библиотек и плагинов требуется разрешение автора.
Метод формальных манипуляций для обращения матриц.
Линейные системы основа всех численных методов. Решение линейных систем проходят еще в школе. Предлагаемый алгоритм обращения матриц отличается просто более тщательной группировкой членов и совмещением операций. Как и LU разложение он обращает матрицу "на месте" , без дополнительной памяти, но более быстрый. Все операции выполняются в одном цикле, нет декомпозиций, прямого-обратного прохода. Число операций ровно N3 - столько же сколько и умножение матриц. Коды много лет используются в коммерческих проектах.
Примерный алгоритм
- Заводятся 2 целых вектора длинной N. Векторы расписываются номерами строк, столбцов.
- Среди всех непомеченных строк, столбцов ищется максимальный элемент Mpq
-
Вычисления
-
Для всех элементов
Mjk-= Mpk*Mjq/Mpq - Для строк и столбцов ведущего элемета
Mjq/= Mpq
Mpj/= Mpq -
сам ведущий
Mpq= -1/ Mpq
-
Для всех элементов
- строка и столбец ведущего помечаются как обработанные (значение -1) и производятся перестановки
-
Если остались необработанные - возврат к новому поиску ведущего
В результате, на месте матрицы, получается обратная матрица с обратным знаком.
Помимо отсутствия обратного хода этот алгоритм отличается еще и тем, что не делает различия между строками и столбцами.
При перестановках строка меняется с номером столбца, а столбец с номером строки, при этом используются исходные номера.
Номера приписываются строкам/столбцам изначально в двух массивах.
Строка идет при перестановке на истинный номер столбца и остается там ,получая номер -1.
Вытесняемая строка/столбец идут на ее место и пока располагаются по несвоему номеру.
Чтоб пояснить перестановки возьмем предельный случай, в каждой строке и столбце по одному ненулевому элементу, все они ведущие.
Для прямого преобразования i->aji*j .
Обратное j->1/aji*i .
Значит 1/aji в обратной матрице стоит на месте ij.
Отсюда ясно где, в конечном счете, должна оказаться строка и столбец ведущего элемента.
Вычисления в текстах эквивалентны, но несколько соптимизированы.
В этих кодах выполняется полный выбор максимального элемента.
Исключительно из за этого он несколько медленнее LU с частичным выбором ведущего, иначе был бы быстрее.
Есть примеры, где частичный выбор ведущего срабатывает плохо.
И я не вижу универсального критерия для частичного выбора, например при симметричной матрице ведущий нужно искать на диагонали, и поэтому частичный выбор ведущего - это уже специализация зависящая от задачи.
Примерно половина времени уходит на поиск ведущего при полном переборе, и при этом скорость процентов на 30 медленнее чем у алгоритмов с LU с частичным выбором.
Обоснование методов – элементарные преобразования + приписывание параллельной единичной матрицы.
При элементарных преобразованиях зануляется элемент в исходной матрице и образуется новый в единичной.
Это позволяет использовать только одну матрицу.
До сих пор известный вариант совмещения этого LU разложение.
Но имеются 2 стратегии – либо занулять только ниже диагонали как в LU, либо занулять сразу весь столбец.
Как оказывается, мой вариант более производителен в случае обращения матриц.
Для решения линейных уравнений рациональней двухпроходный метод.
Но LU-разложение проигрывает обычному методу Гаусса с обратным ходом из-за абстрактности метода.
Производятся операции над коэффициентами, вместо конкретных чисел.
Соотношение 1/3N³ к 2/3N³ на LU-разложение (только LU, без обращения матрицы).
Здесь классическое LU-обращение матриц ( netlib.org ) - помимо использования нерационального LU, эти коды весьма громоздки и обращаются к огромному количеству низкоуровневых процедур написанных на ассемблере.
Но современные компиляторы обычно не нуждаются в этом.
Мат библиотеки от Intel и AMD отличаются большей скоростью.
Они неплохо работают с большими матрицами, которые не помещаются в кэш.
Основная операция это вычитание из строки другой строки, помноженной на число, которое в регистре с обратной запись результата.
Выпадение строк из кэша, а это несколько мегабайт в зависимости от процессора приводит к замедлению на порядок.
Это на одно ядро, тоесть для четырехядерных можно сказать, что производительность падает в сорок раз.
Библиотеки от Intel и AMD неплохо справляются с этой проблемой.
Чтоб уменьшить замедление нужно за один проход по матрице делать как можно больше операций над одной строкой.
Можно искать следующий разрешающий элемент параллельно с преобразованиями и вычитать сразу несколько строк.
Правда, при вычитание нескольких строк за один проход невозможен поиск ведущего по всей матрице.
Возможно в каких-то случаях использовать сравнение кодов вместо сравнения плавающих чисел при поиске ведущего.
Если удалить знаковый бит, то поскольку порядок в старшей части результат плавающего сравнения такой же как и целого, но вряд ли это актуально.
Вот оптимизированное обращение матрицы корреляций, которая симметрична а ведущие элементы положительны и на главной диагонали. Используется верхняя часть, потом переписывается в нижнюю CovarianceMtrix
В этих кодах производится полный разбор всех случаев.
И обращение, и решение СЛАУ возвращают дефект линейного преобразования, детерминант и частичное решение в том случае если детерминант равен 0.
Частичное решение - это решение на подпространстве, где детерминант еще не равен 0.
Это частичное решение, например, подходит для метода наименьших квадратов, даже в том случае если ищется приближение функциями, которые линейно зависимы.
Здесь две версии кодов.
Шаблонная версия кодов отличается тем, что может работать с родными массивами c++ и любыми контейнерами.
При этом допустимо их смешение в одном вызове.
Для подключения любого контейнера достаточно дописать пару строчек.
Реализован оператор определения типа в шаблонах – typeof.
Это решение соответствует последним стандартам и проверено в VC6 – VC8.
Версия double содержит удобные контейнеры, также совместима с простыми массивами.
#include <math.h>
int gaus_ostatok;
// дефект линейного преобразования
// 0 == ОК
long double gaus_deter;
// в случае успешного обращения детерминант
double gaus_minved;
// минимальный ведущий элемент
// можно использовать для оценки точности
#define _NaN() (double&)*"Implementation dependent";
// заполнитель для неразрешенных измерений
// при невозможности обращения
template<int N>
inline static void invers(double m[N][N])
{
// (c) Drobotenko http://drobotenko.com
int rn[N],cn[N];
int j,k;
for(j=N;j--;)
rn[j]=cn[j]=j;
gaus_minved=1e99;
gaus_deter=1;
for(gaus_ostatok=N;gaus_ostatok;gaus_ostatok--)
{ int jved,kved;
double vved=-1,t;
// поиск ведущего
for(j=N;j--;)
{ if(~rn[j])
for(k=N;k--;)
if(~cn[k])
if(vved<fabs(m[j][k]))
vved=fabs(m[j][k]),jved=j,kved=k;
}
if(gaus_minved>vved)
gaus_minved=vved;
gaus_deter*=m[jved][kved];
if(vved<1e-99)
{ for(j=N;j--;)
{ if(~rn[j]) for(k=N;k--;)
m[j][k]=_NaN();
if(~cn[j]) for(k=N;k--;)
m[k][j]=_NaN();
}
return;
}
int jt=rn[jved],kt=cn[kved];
if(jt!=kved)
gaus_deter=-gaus_deter;
// меняем знак детерминанта если перестановка
// перестановки
for(j=N;j--;)
t=m[kt][j],m[kt][j]=m[jved][j],m[jved][j]=t;
for(j=N;j--;)
t=m[j][jt],m[j][jt]=m[j][kved],m[j][kved]=t;
rn[jved]=rn[kt];
cn[kved]=cn[jt];
rn[kt]=cn[jt]=-1;
vved=m[kt][jt]; m[kt][jt]=1;
for(j=N;j--;)
{ if(j==kt)
continue;
double mul=m[j][jt]/vved;
m[j][jt]=0;
for(k=N;k--;)
m[j][k]-=m[kt][k]*mul;
// самый внутренний цикл ровно N^3 операций
}
for(k=N;k--;)
m[kt][k]/=vved;
}
}
const int r=5;
double m[r][r],e[r][r];
double mval(int i,int j)
{ // значения матрицы
return int(1000./(.1+fabs(i-.95*j)))/1000.;
}
int main(int, char **)
{ for(int i=r;i--;)
for(int j=r;j--;)
m[i][j]=mval(i,j);
invers(m);
for(int i=r;i--;)
for(int j=r;j--;)
{ // проверка в Е получается единичная матрица
e[i][j]=0;
for(int k=r;k--;)
e[i][j]+=m[i][k]*mval(k,j);
}
}
Если компилятор не может справится с этим текстом, попробуйте заменить параметр на ссылку void invers(double (&m)[N][N]) , или явно укажите параметр шаблона при вызове invers<r>(m); можно наконец в параметрах шаблона задать и тип матрицы
template<int N, class C>inline static void invers(C &m)
{
}
....................
invers<r,double [r][r]>(m);
возможно это с любым компилятором работать будет, хоть каждый вызов указывать параметры шаблона неудобно - поддержка стандарта не у всех компиляторов на высоте :(
т.к. резмерность массива можно вывести через sizeof удобно написать макрос
#define INVERS(M) invers<sizeof((M)[0])/sizeof((M)[0][0])>(M)
КОДЫ - какаято старая версия