Быстрое точечное изделие для особого случая
Для вектора X размера L, где каждый скалярный элемент из X является бинарным множеством {0,1}, он должен найти точечное произведение z = dot (X, Y), если вектор Y размера L состоит из целочисленных элементов. Я предлагаю, должен существовать очень быстрый способ сделать это.
Скажем, мы имеем L=4; X[L]={1, 0, 0, 1}; Y[L]={-4, 2, 1, 0}
, и нам нужно найти z=X[0]*Y[0] + X[1]*Y[1] + X[2]*Y[2] + X[3]*Y[3]
(что в этом случае даст нам -4
).
Очевидно, что X может быть представлено с использованием двоичных цифр, например. целочисленный тип int32 для L = 32. Тогда все, что нам нужно сделать, это найти точечный продукт этого целого числа с массивом из 32 целых чисел. Есть ли у вас какие-либо идеи или предложения, как это сделать очень быстро?
Ответы
Ответ 1
Это действительно потребует профилирования, но альтернативы, которые вы можете рассмотреть:
int result=0;
int mask=1;
for ( int i = 0; i < L; i++ ){
if ( X & mask ){
result+=Y[i];
}
mask <<= 1;
}
Обычно бит-сдвиг и побитовые операции быстрее, чем умножение, однако оператор if может быть медленнее, чем умножение, хотя с предсказанием ветвления и большим L, я предполагаю, что это может быть быстрее. Вы действительно должны были бы профилировать его, однако, чтобы определить, привело ли это к ускорению.
Как указано в комментариях ниже, разворачивание цикла вручную или через флаг компилятора (например, "-unun-loop-loop" на GCC) также может ускорить это (вернувшись к условию цикла).
Изменить
В комментариях ниже предложена следующая хорошая настройка:
int result=0;
for ( int i = 0; i < L; i++ ){
if ( X & 1 ){
result+=Y[i];
}
X >>= 1;
}
Ответ 2
Помогло ли вам взглянуть на SSE2? Он уже имеет операции типа "точка-продукт", плюс вы можете тривиально сделать 4 (или, возможно, 8, я забыл размер регистра) простую итерацию вашей наивной петли параллельно.
SSE также имеет некоторые простые операции с логическим типом, поэтому он может выполнять дополнения, а не умножения, не используя никаких условных операций... снова вам нужно будет посмотреть, какие операционные системы доступны.
Ответ 3
Попробуйте следующее:
int result=0;
for ( int i = 0; i < L; i++ ){
result+=Y[i] & (~(((X>>i)&1)-1));
}
Это позволяет избежать условного оператора и использует побитовые операторы для маскировки скалярного значения с нулями или единицами.
Ответ 4
Поскольку размер явно не имеет значения, я думаю, что следующий, вероятно, наиболее эффективный код общего назначения:
int result = 0;
for (size_t i = 0; i < 32; ++i)
result += Y[i] & -X[i];
Бит-кодирование X
просто ничего не приносит в таблицу (даже если цикл может завершиться раньше, чем @Mathieu правильно отметил). Но опускание if
внутри цикла делает.
Конечно, разворот цикла может резко ускорить это, как отмечали другие.
Ответ 5
Это решение идентично, но немного быстрее (по моему тесту), чем у Micheal Aaron's:
long Lev=1;
long Result=0
for (int i=0;i<L;i++) {
if (X & Lev)
Result+=Y[i];
Lev*=2;
}
Я думал, что существует численный способ быстро установить следующий бит в слове, который должен повысить производительность, если ваши X-данные очень разрежены, но в настоящее время не могут найти указанную численную формулировку в настоящее время.
Ответ 6
Я видел несколько ответов с битовым обманом (чтобы избежать ветвления), но ни один не получил петлю вправо imho:/
Оптимизация @Goz
ответ:
int result=0;
for (int i = 0, x = X; x > 0; ++i, x>>= 1 )
{
result += Y[i] & -(int)(x & 1);
}
Преимущества:
- не нужно выполнять операции бит-сдвига
i
каждый раз (X>>i
)
- цикл останавливается раньше, если
X
содержит 0 в более высоких битах
Теперь я задаюсь вопросом, работает ли он быстрее, тем более, что преждевременная остановка цикла for может быть не такой простой для разворачивания цикла (по сравнению с константой времени компиляции).
Ответ 7
Как насчет объединения цикла переключения с небольшой таблицей поиска?
int result=0;
for ( int x=X; x!=0; x>>=4 ){
switch (x&15) {
case 0: break;
case 1: result+=Y[0]; break;
case 2: result+=Y[1]; break;
case 3: result+=Y[0]+Y[1]; break;
case 4: result+=Y[2]; break;
case 5: result+=Y[0]+Y[2]; break;
case 6: result+=Y[1]+Y[2]; break;
case 7: result+=Y[0]+Y[1]+Y[2]; break;
case 8: result+=Y[3]; break;
case 9: result+=Y[0]+Y[3]; break;
case 10: result+=Y[1]+Y[3]; break;
case 11: result+=Y[0]+Y[1]+Y[3]; break;
case 12: result+=Y[2]+Y[3]; break;
case 13: result+=Y[0]+Y[2]+Y[3]; break;
case 14: result+=Y[1]+Y[2]+Y[3]; break;
case 15: result+=Y[0]+Y[1]+Y[2]+Y[3]; break;
}
Y+=4;
}
Производительность этого будет зависеть от того, насколько хорош компилятор для оптимизации оператора switch, но по моему опыту они довольно хороши в этом в наши дни....
Ответ 8
Вероятно, нет общего ответа на этот вопрос. Вы должны прокомментировать свой код под разными целями. Производительность будет зависеть от оптимизации компилятора, такой как разматывание цикла и инструкции SIMD, которые доступны на большинстве современных процессоров (x86, PPC, ARM имеют свои собственные реализации).
Ответ 9
result = 0;
for(int i = 0; i < L ; i++)
if(X[i]!=0)
result += Y[i];
Ответ 10
Для small L вы можете использовать оператор switch вместо цикла. Например, если L = 8, вы могли бы:
int dot8(unsigned int X, const int Y[])
{
switch (X)
{
case 0: return 0;
case 1: return Y[0];
case 2: return Y[1];
case 3: return Y[0]+Y[1];
// ...
case 255: return Y[0]+Y[1]+Y[2]+Y[3]+Y[4]+Y[5]+Y[6]+Y[7];
}
assert(0 && "X too big");
}
И если L = 32, вы можете написать функцию dot32(), которая вызывает dot8() четыре раза, если возможно, по необходимости. (Если ваш компилятор отказывается от встроенного dot8(), вы можете переписать dot8() в качестве макроса для принудительной вставки.) Добавлено:
int dot32(unsigned int X, const int Y[])
{
return dot8(X >> 0 & 255, Y + 0) +
dot8(X >> 8 & 255, Y + 8) +
dot8(X >> 16 & 255, Y + 16) +
dot8(X >> 24 & 255, Y + 24);
}
Это решение, как указывает Mikera, может иметь стоимость кэша команд; если это так, использование функции dot4() может помочь.
Дальнейшее обновление. Это можно комбинировать с решением mikera:
static int dot4(unsigned int X, const int Y[])
{
switch (X)
{
case 0: return 0;
case 1: return Y[0];
case 2: return Y[1];
case 3: return Y[0]+Y[1];
//...
case 15: return Y[0]+Y[1]+Y[2]+Y[3];
}
}
Посмотрев на полученный код ассемблера с параметрами -S -O3 с gcc 4.3.4 на CYGWIN, я немного удивлен, увидев, что это автоматически встроено в dot32(), с восемь 16-позиционных таблиц перехода.
Но добавление __attribute __ ((__ noinline__)), похоже, создает привлекательный ассемблер.
Другим вариантом является использование провалов в инструкции switch, но gcc добавляет инструкции jmp, и он не выглядит быстрее.
Изменить - Полностью новый ответ: После размышления о 100-часовом штрафе, упомянутом Ants Aasma, и других ответах, вышеупомянутое, вероятно, не оптимально. Вместо этого вы можете вручную развернуть цикл, как в:
int dot(unsigned int X, const int Y[])
{
return (Y[0] & -!!(X & 1<<0)) +
(Y[1] & -!!(X & 1<<1)) +
(Y[2] & -!!(X & 1<<2)) +
(Y[3] & -!!(X & 1<<3)) +
//...
(Y[31] & -!!(X & 1<<31));
}
Это, на моей машине, генерирует 32 x 5 = 160 быстрых инструкций. Разумный компилятор мог бы, возможно, развернуть другие предложенные ответы, чтобы дать тот же результат.
Но я все еще проверяю.
Ответ 11
Весьма вероятно, что время, затраченное на загрузку X
и Y
из основной памяти, будет доминировать. Если это имеет место для вашей архитектуры процессора, алгоритм работает быстрее при загрузке меньше. Это означает, что сохранение X
в виде битовой маски и расширение его в кеш L1 ускорит алгоритм в целом.
Другой важный вопрос заключается в том, будет ли ваш компилятор создавать оптимальные нагрузки для Y
. Это зависит от процессора и компилятора. Но в целом это помогает, если компилятор может видеть, какие значения необходимы, когда. Вы можете вручную развернуть цикл. Однако, если L является contant, оставьте его компилятору:
template<int I> inline void calcZ(int (&X)[L], int(&Y)[L], int &Z) {
Z += X[I] * Y[I]; // Essentially free, as it operates in parallel with loads.
calcZ<I-1>(X,Y,Z);
}
template< > inline void calcZ<0>(int (&X)[L], int(&Y)[L], int &Z) {
Z += X[0] * Y[0];
}
inline int calcZ(int (&X)[L], int(&Y)[L]) {
int Z = 0;
calcZ<L-1>(X,Y,Z);
return Z;
}
(Конрад Рудольф поставил под сомнение это в комментарии, размышляя об использовании памяти.Это не настоящее узкое место в современных компьютерных архитектурах, пропускная способность между памятью и процессором. Этот ответ почти не имеет значения, если Y как-то уже находится в кеше.)
Ответ 12
Вы можете сохранить свой бит-вектор как последовательность int, где каждый int упаковывает пару коэффициентов в виде битов. Тогда компонентное умножение эквивалентно бит-и. При этом вам просто нужно подсчитать количество заданных битов, которое можно сделать следующим образом:
inline int count(uint32_t x) {
// see link
}
int dot(uint32_t a, uint32_t b) {
return count(a & b);
}
Немного взломать отсчет установленных бит см. http://graphics.stanford.edu/~seander/bithacks.html#CountBitsSetParallel
Изменить: Извините, я только понял, что только один из векторов содержит элементы {0,1}, а другой - нет. Этот ответ применим только к случаю, когда оба вектора ограничены коэффициентами из набора {0,1}.
Ответ 13
Ну, вы хотите, чтобы все биты прошли, если они равны 1 и никому, если это 0. Итак, вы хотите как-то превратить 1 в -1 (т.е. 0xffffffff), и 0 останется прежним. Это просто -X.... так что вы...
Y & (-X)
для каждого элемента... работа выполнена?
Edit2: Чтобы дать пример кода, вы можете сделать что-то вроде этого и избежать ветки:
int result=0;
for ( int i = 0; i < L; i++ )
{
result+=Y[i] & -(int)((X >> i) & 1);
}
Конечно, вам лучше всего оставить 1s и 0s в массиве int и, следовательно, избегать сдвигов.
Изменить: также стоит отметить, что если значения в Y имеют размер 16 бит, вы можете сделать 2 из них и операции для каждой операции (4, если у вас есть 64-разрядные регистры). Это означает отрицание значений X 1 на 1 в большее целое число.
т.е. YVals = -4, 3 в 16-бит = 0xFFFC, 0x3... помещается в 1 32-бит, и вы получаете 0xFFFC0003. Если у вас есть 1, 0 как X vals, тогда вы формируете битную маску 0xFFFF0000 и 2 вместе, и у вас есть 2 результата в 1 побитовом и op.
Другое редактирование:
ЕСЛИ вы хотите, чтобы код о том, как сделать второй метод, должен работать (хотя он использует неуказанное поведение, поэтому он не может работать на каждом компиляторе.. работает на каждом компиляторе, с которым я столкнулся).
union int1632
{
int32_t i32;
int16_t i16[2];
};
int result=0;
for ( int i = 0; i < (L & ~0x1); i += 2 )
{
int3264 y3264;
y3264.i16[0] = Y[i + 0];
y3264.i16[1] = Y[i + 1];
int3264 x3264;
x3264.i16[0] = -(int16_t)((X >> (i + 0)) & 1);
x3264.i16[1] = -(int16_t)((X >> (i + 1)) & 1);
int3264 res3264;
res3264.i32 = y3264.i32 & x3264.i32;
result += res3264.i16[0] + res3264.i16[1];
}
if ( i < L )
result+=Y[i] & -(int)((X >> i) & 1);
Надеюсь, что компилятор оптимизирует назначение (не в верхней части моей головы, я не уверен, но идея может быть переработана так, чтобы они определенно были) и дать вам небольшую скорость в том, что вам сейчас нужно только сделать 1 побитовое - и вместо 2. Ускорение будет незначительным, хотя...
Ответ 14
Представляем X
, используя связанный список мест, где x[i] = 1
.
Чтобы найти требуемую сумму, вам нужны операции O(N)
, где N
- размер вашего списка.