UUT-разложение Холецкого
Продолжаю биться головой об стену, искать наиболее подходящий способ работы с симметричной положительно определённой матрицей 6х6 в целых числах.
Первая идея была - LDLT-разложение, за ним обращение унитреугольной матрицы и получение обратной исходной матрицы (или решение системы уравнений). Также между делом разобрались с индийским методом и обнаружили, что он давным-давно был реализован на фортране.
Потом подумалось, что для нашего процессора удобнее сделать UDUT-разложение - концептуально то же самое, но все циклы делаются "снизу вверх", от бОльших значений до нуля, что позволяет использовать "инструкции" iLoop/jLoop/kLoop (если регистр больше нуля - вычесть единичку и перейти по заданному адресу, иначе ничего не делать).
Но увы, оба метода оказались очень "капризными" - как я ни пытаюсь масштабировать неизвестные, а при разложении обязательно где-то нарываюсь на переполнение! В общем-то, задним числом я мог это предвидеть... Если исходная матрица:
)
то в процессе разложения диагональные элементы могут только УМЕНЬШАТЬСЯ. И хорошо, если не до нуля (и того ниже), поскольку это будет означать, что матрица не является положительно определённой. Возможно, она была такой, но ошибки вычислений изменили её вид.
А поскольку при нахождении элементов ВНЕ диагонали, их рано или поздно делим на диагональные, то и получаем огромный размах промежуточных значений, который упорно не лезет в наш диапазон [-1;+1].
Ясно, что можно избежать такой проблемы, если нарастить разрядность чисел, либо перейти к плавающей точке, либо как-то более удачно ввести 6 степеней свободы.
Но есть и ещё одна идея: может, надо перестать бояться и полюбить квадратный корень? Да, операция не самая приятная в плане реализации в целых числах, но зато у нас одним махом падает разброс между значениями, какой наблюдается сейчас: вместо соотношения почти 1 : 10000 (еле-еле помещается в 16 бит), у нас выйдет всего лишь 1 : 100, и теперь операции деления станут совсем не страшны.
Вместо этих LDLT и UDUT-разложений, которые были придуманы, чтобы избавиться от квадратных корней, попробуем старое доброе разложение Холецкого, но не классический вариант LLT, а вариант "задом наперёд" UUT для нашего Quat Core.
Чтобы вконец не одуреть, вбивая формулы TeX, снова ограничимся матрицей 4×4 - на ней уже видны все 3 вложенных цикла.
Нужно представить симметричную положительно определённую матрицу A в виде произведения верхней треугольной матрицы U на себя же транспонированную:
\left(\begin{array}{cccc}u_{00}&0&0&0\\u_{01}&u_{11}&0&0\\u_{02}&u_{12}&u_{22}&0\\u_{03}&u_{13}&u_{23}&u_{33}\end{array}\right))
Перемножим их:
)
Как и должно быть, результат умножения - симметричная матрица. Для наглядности, уберём повторяющиеся значения из верхнего треугольника:
=\left(\begin{array}{cccc}u_{00}^2+u_{01}^2+u_{02}^2+u_{03}^2&-&-&-\\u_{01}u_{11}+u_{02}u_{12}+u_{03}u_{13}&u_{11}^2+u_{12}^2+u_{13}^2&-&-\\u_{02}u_{22}+u_{03}u_{23}&u_{12}u_{22}+u_{13}u_{23}&u_{22}^2+u_{23}^2&-\\u_{03}u_{33}&u_{13}u_{33}&u_{23}u_{33}&u_{33}^2\end{array}\right))
Здесь соотношения даже проще, чем при LDLT или UDUT-разложении, никаких временных массивов вводить не надо, и работать можно как по строкам, так и по столбцам. Главное - начинать с правого нижнего угла и двигаться к левому верхнему.
Проще всего, видимо, начать с нижней строки, пройдя по ней справа налево, перейти строкой выше, и так далее.
Каждый раз мы начинаем с элемента на главной диагонали, он обрабатывается особо:
Далее идут "обычные" элементы:
После этого значение u33 нам больше не понадобится. Следовательно, мы можем сразу же найти обратную к нему величину, и все деления заменить на умножения. Пока этого не делаю, потому что не очень понятно, как поступать с масштабом. Если исходно значения на диагонали были в диапазоне от 1/32768 до 1 (мы работает с числами 1.15, т.е с помощью целочисленных 16 бит изображаем значения от -1 до 1-2-15), после извлечения квадратного корня они окажутся в диапазоне 0,0055..1, значит, после взятия обратной величины получим диапазон 1..181. Потом ещё подумаем, а пока что просто будем делить два целых числа одно на другое и надеяться, что переполнения не случится.
Переходим на строку выше:
/u_{22},)
/u_{22})
Всё хорошо: в каждой формуле в правой части только те величины, которые нам уже известны, причём исходное значение aij используется всего один раз, после чего вместо него мы используем uji. Значит, это преобразование элементарно делается "на месте".
И ещё на строку выше:
/u_{13})
и наконец, один элемент верхней строки:
Каждый раз, когда мы работаем с элементом (i,j), i-номер строки, j-номер столбца, мы берём его начальное значение и вычитаем скалярное произведение двух векторов - первый образован элементами ПОД текущим. Это элементы (i+1,j), (i+2,j), ... (N-1, j).
Второй вектор образован из элементов (i+1, i), (i+2,i), ..., (N-1, i), см. рисунок. Мы как бы идём от текущего элемента вправо, затем "отражаемся" от главной диагонали:
В матрице N×N нужно произвести лишь N извлечений квадратного корня, так что даже если эта операция окажется трудозатратной, вряд ли именно она окажет решающее влияние на производительность данного алгоритма.
Посмотрим, как выглядит UUT-разложение для матрицы из начала поста:
)
Такая матрица производит хорошее впечатление - на главной диагонали нет совсем маленьких, околонулевых значений, и в целом почти все значения имеют один порядок величины.
Если попытаться для той же матрицы выполнить LDLT-разложение, мы "упрёмся" почти что сразу, попытавшись заменить a40=478 (что означает 478/32768≈0,0146) на 478/38, что СИЛЬНО больше единицы, из-за чего получим переполнение. UDUT-разложение провести можно (это одна из причин, почему я начал его рассматривать), получим такую матрицу (в ней объединены нижняя унитреугольная и диагональная):
)
Как видно, на главной диагонали стоит две единицы (то есть значения 1/32768), а это означает, что ненулевые значения вектора 6х1 в этих позициях при решении системы линейных уравнений сразу же приведут к переполнению.
Элементы вне диагонали существенно возросли по амплитуде, тоже едва не приводя к переполнению, причём, если я начинаю масштабировать первые 3 параметра, то переполнение наступает.
Так что по первому впечатлению, UUT (или LLT) - разложение в данном случае лучше подойдёт. Посмотрим...
Первая идея была - LDLT-разложение, за ним обращение унитреугольной матрицы и получение обратной исходной матрицы (или решение системы уравнений). Также между делом разобрались с индийским методом и обнаружили, что он давным-давно был реализован на фортране.
Потом подумалось, что для нашего процессора удобнее сделать UDUT-разложение - концептуально то же самое, но все циклы делаются "снизу вверх", от бОльших значений до нуля, что позволяет использовать "инструкции" iLoop/jLoop/kLoop (если регистр больше нуля - вычесть единичку и перейти по заданному адресу, иначе ничего не делать).
Но увы, оба метода оказались очень "капризными" - как я ни пытаюсь масштабировать неизвестные, а при разложении обязательно где-то нарываюсь на переполнение! В общем-то, задним числом я мог это предвидеть... Если исходная матрица:
то в процессе разложения диагональные элементы могут только УМЕНЬШАТЬСЯ. И хорошо, если не до нуля (и того ниже), поскольку это будет означать, что матрица не является положительно определённой. Возможно, она была такой, но ошибки вычислений изменили её вид.
А поскольку при нахождении элементов ВНЕ диагонали, их рано или поздно делим на диагональные, то и получаем огромный размах промежуточных значений, который упорно не лезет в наш диапазон [-1;+1].
Ясно, что можно избежать такой проблемы, если нарастить разрядность чисел, либо перейти к плавающей точке, либо как-то более удачно ввести 6 степеней свободы.
Но есть и ещё одна идея: может, надо перестать бояться и полюбить квадратный корень? Да, операция не самая приятная в плане реализации в целых числах, но зато у нас одним махом падает разброс между значениями, какой наблюдается сейчас: вместо соотношения почти 1 : 10000 (еле-еле помещается в 16 бит), у нас выйдет всего лишь 1 : 100, и теперь операции деления станут совсем не страшны.
Вместо этих LDLT и UDUT-разложений, которые были придуманы, чтобы избавиться от квадратных корней, попробуем старое доброе разложение Холецкого, но не классический вариант LLT, а вариант "задом наперёд" UUT для нашего Quat Core.
Чтобы вконец не одуреть, вбивая формулы TeX, снова ограничимся матрицей 4×4 - на ней уже видны все 3 вложенных цикла.
Нужно представить симметричную положительно определённую матрицу A в виде произведения верхней треугольной матрицы U на себя же транспонированную:
Перемножим их:
Как и должно быть, результат умножения - симметричная матрица. Для наглядности, уберём повторяющиеся значения из верхнего треугольника:
Здесь соотношения даже проще, чем при LDLT или UDUT-разложении, никаких временных массивов вводить не надо, и работать можно как по строкам, так и по столбцам. Главное - начинать с правого нижнего угла и двигаться к левому верхнему.
Проще всего, видимо, начать с нижней строки, пройдя по ней справа налево, перейти строкой выше, и так далее.
Каждый раз мы начинаем с элемента на главной диагонали, он обрабатывается особо:
Далее идут "обычные" элементы:
После этого значение u33 нам больше не понадобится. Следовательно, мы можем сразу же найти обратную к нему величину, и все деления заменить на умножения. Пока этого не делаю, потому что не очень понятно, как поступать с масштабом. Если исходно значения на диагонали были в диапазоне от 1/32768 до 1 (мы работает с числами 1.15, т.е с помощью целочисленных 16 бит изображаем значения от -1 до 1-2-15), после извлечения квадратного корня они окажутся в диапазоне 0,0055..1, значит, после взятия обратной величины получим диапазон 1..181. Потом ещё подумаем, а пока что просто будем делить два целых числа одно на другое и надеяться, что переполнения не случится.
Переходим на строку выше:
Всё хорошо: в каждой формуле в правой части только те величины, которые нам уже известны, причём исходное значение aij используется всего один раз, после чего вместо него мы используем uji. Значит, это преобразование элементарно делается "на месте".
И ещё на строку выше:
и наконец, один элемент верхней строки:
Каждый раз, когда мы работаем с элементом (i,j), i-номер строки, j-номер столбца, мы берём его начальное значение и вычитаем скалярное произведение двух векторов - первый образован элементами ПОД текущим. Это элементы (i+1,j), (i+2,j), ... (N-1, j).
Второй вектор образован из элементов (i+1, i), (i+2,i), ..., (N-1, i), см. рисунок. Мы как бы идём от текущего элемента вправо, затем "отражаемся" от главной диагонали:
В матрице N×N нужно произвести лишь N извлечений квадратного корня, так что даже если эта операция окажется трудозатратной, вряд ли именно она окажет решающее влияние на производительность данного алгоритма.
Посмотрим, как выглядит UUT-разложение для матрицы из начала поста:
Такая матрица производит хорошее впечатление - на главной диагонали нет совсем маленьких, околонулевых значений, и в целом почти все значения имеют один порядок величины.
Если попытаться для той же матрицы выполнить LDLT-разложение, мы "упрёмся" почти что сразу, попытавшись заменить a40=478 (что означает 478/32768≈0,0146) на 478/38, что СИЛЬНО больше единицы, из-за чего получим переполнение. UDUT-разложение провести можно (это одна из причин, почему я начал его рассматривать), получим такую матрицу (в ней объединены нижняя унитреугольная и диагональная):
Как видно, на главной диагонали стоит две единицы (то есть значения 1/32768), а это означает, что ненулевые значения вектора 6х1 в этих позициях при решении системы линейных уравнений сразу же приведут к переполнению.
Элементы вне диагонали существенно возросли по амплитуде, тоже едва не приводя к переполнению, причём, если я начинаю масштабировать первые 3 параметра, то переполнение наступает.
Так что по первому впечатлению, UUT (или LLT) - разложение в данном случае лучше подойдёт. Посмотрим...