UDUT-разложение (оно же LTDL)
Это как LDLT-разложение симметричной положительно определённой матрицы A, но только UDUT, т.е верхняя унитреугольная матрица (треугольная с единицами на главной диагонали) умножается на диагональную и на себя транспонированную. Его можно было бы также назвать LTDL-разложением, тут устоявшейся терминологии нет.
Такое разложение ничем принципиально не отличается от LDLT, но для этого разложения и всех последующих операций при инвертировании матрицы A, мы будем идти по элементам "в обратном порядке", с правого нижнего угла к левому верхнему. Для системы команд QuatCore, где циклы "вниз" организуются проще всего с помощью "команд" iLOOP, jLOOP и kLOOP (если индексная переменная равна нулю - идти дальше, иначе уменьшить её на единицу и прыгнуть в начало цикла), это предпочтительнее.
Итак, у нас опять есть матрица A - симметричная, положительно определённая:
=\left(\begin{array}{cccc}a_{00}&-&-&-\\a_{10}&a_{11}&-&-\\a_{20}&a_{21}&a_{22}&-\\a_{30}&a_{31}&a_{32}&a_{33}\end{array}\right))
В этот раз мы ведём нумерацию индексов с нуля, чтобы стать "ближе к железу". Верхний треугольник мы не записываем, чтобы напомнить: в памяти хранится только нижний треугольник. Это не обязательно означает треугольную адресацию:
}{2}+j])
когда нижний треугольник занимал бы последовательные ячейки памяти
)
тоже была идея реализовать такую адресацию "аппаратно", просто в виде таблицы:
itri[i]
00
11
23
36
410
515
(у меня матриц крупнее 6×6 здесь не ожидается)
Это весьма экономично.
[Код на верилоге]
module TriangleNum (input [2:0] D, output [3:0] Q);
assign Q = (D == 0)? 4'd0:
(D == 1)? 4'd1:
(D == 2)? 4'd3:
(D == 3)? 4'd6:
(D == 4)? 4'd10:
(D == 5)? 4'd15:
4'bxxxx;
endmodule
Как ни странно, этот код синтезируется в 3 ЛЭ, так как выходит Q[3] = D[2].
Есть в такой адресации определённая красота: в формуле, преобразующей индексы в адрес в памяти, нигде не фигурируют размеры матрицы! Но можно и самой обычной адресацией обойтись, а в верхний треугольник запихать прочие величины, как в таблице памяти данных ВИПС. Буду ещё громко думать, что лучше.
Далее мы хотим эту матрицу A представить в виде произведения:
\left(\begin{array}{cccc}d_0&0&0&0\\0&d_1&0&0\\0&0&d_2&0\\0&0&0&d_3\end{array}\right)\left(\begin{array}{cccc}1&0&0&0\\u_{01}&1&0&0\\u_{02}&u_{12}&1&0\\u_{03}&u_{13}&u_{23}&1\end{array}\right))
Две эти матрицы умещаются в памяти, выделенной под матрицу A:
)
и это разложение может быть выполнено "на месте", не требуя временного хранилища.
Чтобы понять, как это сделать, помножим 3 матрицы и посмотрим, как их коэффициенты соотносятся с коэффициентами исходной матрицы:
\left(\begin{array}{cccc}d_0&0&0&0\\d_1u_{01}&d_1&0&0\\d_2u_{02}&d_2u_{12}&d_2&0\\d_3u_{03}&d_3u_{13}&d_3u_{23}&d_3\end{array}\right)=)
)
Разумеется, и здесь мы получили симметричную матрицу, так что уберём ненужный верхний треугольник:
)
Как и ожидалось, в этот раз самые простые выражения расположены в правом нижнем углу и усложняются по мере движения к левому верхнему (в прошлый раз всё было наоборот).
Можно идти в разном порядке. В прошлый раз мы двигались по столбцам, но как оказалось, для эффективных вычислений нужно было вводить временные переменные. Будем двигаться по столбцам и в этот раз - упорство должно быть вознаграждено!
Начинаем с правого столбца.
При работе "на месте" делать это присвоение не требуется - d3 расположен в точности там же, где a33.
Берёмся за следующий (предпоследний) столбец, причём двигаемся снизу вверх:
откуда:
Берёмся за следующий столбец:
откуда:
Наконец, обрабатываем левый столбец:
откуда:
Вот и всё. Глядя на эти формулы, можно обнаружить, что введя вспомогательный массив из N-1 значений (N×N-размер матрицы, в рассмотренном случае N=4), мы обойдёмся без львиной доли умножений.
Покажем, как это делается. Запишем выражения для предпоследнего столбца:
Для столбца с индексом 1:
И наконец, для левого столбца (с индексом 0):
Выпишем все эти выкладки в таблицу. Вычисления здесь ведутся сверху вниз и только потом слева направо:
В такой записи получается, что значения dn нам нужны только в знаменателе. Если мы делаем это разложение, чтобы обратить матрицу "на месте", то мы можем прямо "на лету" заносить вместо самих значений dn обратные, как-то так:
}=1/t_3)
},)
},)
},)
}=1/t_2)
},)
},)
}=1/t_1)
}=1/t_0)
В такой записи нам нужно всего лишь N взятий обратной величины, так что даже при отсутствии аппаратного деления нас это не шибко затормозит (метод Ньютона - наше всё!). Размер вспомогательного массива возрос до N. По идее, можно обращать dj "сразу же", не запихивая их сначала в tj, но так получается более логичный код.
Все остальные выражения очень хорошо ложатся на наше "умножение с накоплением" (Fused Multiply-Add, FMA, и Fused Multiply-Subtract, FMS), тогда как исходные выражения этому не удовлетворяли, поскольку требовали перемножить 3 значения, и лишь затем результат прибавлять к результату. В результате мы не только получаем компактный и быстрый код, он ещё и обладает хорошей точностью за счёт того, что накопление производится во всех 32 битах, и лишь результат записывается в 16 бит.
Не будем подсчитывать отдельно сложения и умножения, нас интересует число операций FMA/FMS.
При обработке столбца под номером n, нам требуется n-1 умножений, чтобы временные значения tj помножить на диагональные элементы. Ещё мы вычисляем n временных значений, причём на первое не требуется арифметических операций вообще, на второе - одна операция FMS, на третье - 2 операции, и на значение n: n-1 операций. В общей сложности это получается n(n-1)/2, и ещё n-1, это получается (n-1)(n+2)/2 = n2/2+n/2-1.
Просуммировав n от 1 до N, получаем количество операций FMA/FMS: N3/6+N2/2-(2/3)N. Асимптотически совпадает с тем, что мы посчитали в LDLT-разложении (т.е коэф 1/6 при N3), но на малых значениях N чувствуется существенный выигрыш:
При N=4 в прошлый раз у нас вышло 6 делений, 16 умножений и 10 сложений. В этот раз: 4 взятия обратной величины и 16 операций FMA/FMS.
При N=6 в прошлый раз получалось 15 делений, 50 умножений и 35 сложений, а теперь: 6 взятий обратной величины и 50 операций FMA/FMS.
Причём, это ЕДИНСТВЕННЫЕ 6 взятий обратной величины, которые нужны для обращения матрицы A, больше не будет!
Понятно, что разница здесь не в том, используется ли LDLT или UDUT, а в последовательности выполнения операций и способа введения вспомогательного массива.
Такое разложение ничем принципиально не отличается от LDLT, но для этого разложения и всех последующих операций при инвертировании матрицы A, мы будем идти по элементам "в обратном порядке", с правого нижнего угла к левому верхнему. Для системы команд QuatCore, где циклы "вниз" организуются проще всего с помощью "команд" iLOOP, jLOOP и kLOOP (если индексная переменная равна нулю - идти дальше, иначе уменьшить её на единицу и прыгнуть в начало цикла), это предпочтительнее.
Итак, у нас опять есть матрица A - симметричная, положительно определённая:
В этот раз мы ведём нумерацию индексов с нуля, чтобы стать "ближе к железу". Верхний треугольник мы не записываем, чтобы напомнить: в памяти хранится только нижний треугольник. Это не обязательно означает треугольную адресацию:
когда нижний треугольник занимал бы последовательные ячейки памяти
тоже была идея реализовать такую адресацию "аппаратно", просто в виде таблицы:
itri[i]
00
11
23
36
410
515
(у меня матриц крупнее 6×6 здесь не ожидается)
Это весьма экономично.
[Код на верилоге]
module TriangleNum (input [2:0] D, output [3:0] Q);
assign Q = (D == 0)? 4'd0:
(D == 1)? 4'd1:
(D == 2)? 4'd3:
(D == 3)? 4'd6:
(D == 4)? 4'd10:
(D == 5)? 4'd15:
4'bxxxx;
endmodule
Как ни странно, этот код синтезируется в 3 ЛЭ, так как выходит Q[3] = D[2].
Есть в такой адресации определённая красота: в формуле, преобразующей индексы в адрес в памяти, нигде не фигурируют размеры матрицы! Но можно и самой обычной адресацией обойтись, а в верхний треугольник запихать прочие величины, как в таблице памяти данных ВИПС. Буду ещё громко думать, что лучше.
Далее мы хотим эту матрицу A представить в виде произведения:
Две эти матрицы умещаются в памяти, выделенной под матрицу A:
и это разложение может быть выполнено "на месте", не требуя временного хранилища.
Чтобы понять, как это сделать, помножим 3 матрицы и посмотрим, как их коэффициенты соотносятся с коэффициентами исходной матрицы:
Разумеется, и здесь мы получили симметричную матрицу, так что уберём ненужный верхний треугольник:
Как и ожидалось, в этот раз самые простые выражения расположены в правом нижнем углу и усложняются по мере движения к левому верхнему (в прошлый раз всё было наоборот).
Можно идти в разном порядке. В прошлый раз мы двигались по столбцам, но как оказалось, для эффективных вычислений нужно было вводить временные переменные. Будем двигаться по столбцам и в этот раз - упорство должно быть вознаграждено!
Начинаем с правого столбца.
При работе "на месте" делать это присвоение не требуется - d3 расположен в точности там же, где a33.
Берёмся за следующий (предпоследний) столбец, причём двигаемся снизу вверх:
откуда:
Берёмся за следующий столбец:
откуда:
Наконец, обрабатываем левый столбец:
откуда:
Вот и всё. Глядя на эти формулы, можно обнаружить, что введя вспомогательный массив из N-1 значений (N×N-размер матрицы, в рассмотренном случае N=4), мы обойдёмся без львиной доли умножений.
Покажем, как это делается. Запишем выражения для предпоследнего столбца:
Для столбца с индексом 1:
И наконец, для левого столбца (с индексом 0):
Выпишем все эти выкладки в таблицу. Вычисления здесь ведутся сверху вниз и только потом слева направо:
В такой записи получается, что значения dn нам нужны только в знаменателе. Если мы делаем это разложение, чтобы обратить матрицу "на месте", то мы можем прямо "на лету" заносить вместо самих значений dn обратные, как-то так:
В такой записи нам нужно всего лишь N взятий обратной величины, так что даже при отсутствии аппаратного деления нас это не шибко затормозит (метод Ньютона - наше всё!). Размер вспомогательного массива возрос до N. По идее, можно обращать dj "сразу же", не запихивая их сначала в tj, но так получается более логичный код.
Все остальные выражения очень хорошо ложатся на наше "умножение с накоплением" (Fused Multiply-Add, FMA, и Fused Multiply-Subtract, FMS), тогда как исходные выражения этому не удовлетворяли, поскольку требовали перемножить 3 значения, и лишь затем результат прибавлять к результату. В результате мы не только получаем компактный и быстрый код, он ещё и обладает хорошей точностью за счёт того, что накопление производится во всех 32 битах, и лишь результат записывается в 16 бит.
Не будем подсчитывать отдельно сложения и умножения, нас интересует число операций FMA/FMS.
При обработке столбца под номером n, нам требуется n-1 умножений, чтобы временные значения tj помножить на диагональные элементы. Ещё мы вычисляем n временных значений, причём на первое не требуется арифметических операций вообще, на второе - одна операция FMS, на третье - 2 операции, и на значение n: n-1 операций. В общей сложности это получается n(n-1)/2, и ещё n-1, это получается (n-1)(n+2)/2 = n2/2+n/2-1.
Просуммировав n от 1 до N, получаем количество операций FMA/FMS: N3/6+N2/2-(2/3)N. Асимптотически совпадает с тем, что мы посчитали в LDLT-разложении (т.е коэф 1/6 при N3), но на малых значениях N чувствуется существенный выигрыш:
При N=4 в прошлый раз у нас вышло 6 делений, 16 умножений и 10 сложений. В этот раз: 4 взятия обратной величины и 16 операций FMA/FMS.
При N=6 в прошлый раз получалось 15 делений, 50 умножений и 35 сложений, а теперь: 6 взятий обратной величины и 50 операций FMA/FMS.
Причём, это ЕДИНСТВЕННЫЕ 6 взятий обратной величины, которые нужны для обращения матрицы A, больше не будет!
Понятно, что разница здесь не в том, используется ли LDLT или UDUT, а в последовательности выполнения операций и способа введения вспомогательного массива.