самая быстрая линейная система для небольших квадратных матриц (10x10)

9

Мне очень интересно оптимизировать решение линейных систем для маленьких матриц (10x10), которые иногда называют крошечными матрицами. Есть ли готовое решение для этого? Матрицу можно считать неособой.

Этот решатель должен выполняться более 1 000 000 раз в микросекундах на процессоре Intel. Я говорю об уровне оптимизации, используемом в компьютерных играх. Неважно, кодирую ли я его в сборке и для конкретной архитектуры, или изучаю снижение компромиссов между точностью и надежностью и использую хаки с плавающей запятой (я использую флаг компиляции -ffast-math, нет проблем). Решение может даже потерпеть неудачу примерно в 20% случаев!

Eigen'spartalPivLu является самым быстрым в моем текущем тесте, превосходя LAPACK, когда оптимизирован с -O3 и хорошим компилятором. Но сейчас я нахожусь на стадии ручного изготовления специального линейного решателя. Любой совет будет принята с благодарностью. Я сделаю свое решение открытым исходным кодом и узнаю ключевые идеи публикаций и т. Д.

Связанный: Скорость решения линейной системы с блочной диагональной матрицей Какой самый быстрый способ инвертировать миллионы матриц? https://stackoverflow.com/q/50909385/1489510

rfabbri
источник
7
Это похоже на растянутую цель. Давайте предположим, что мы используем самый быстрый Skylake-X Xeon Platinum 8180 с теоретической пиковой пропускной способностью 4 TFLOP с одинарной точностью, и что для одной системы 10x10 требуется около 700 (примерно 2n ** 3/3) операций с плавающей запятой для решения. Тогда партия 1M таких систем теоретически может быть решена за 175 микросекунд. Это число не может превышать скорость света. Можете ли вы поделиться тем, какую производительность вы достигли в настоящее время с вашим самым быстрым существующим кодом? Кстати, данные одинарной или двойной точности?
njuffa
@njuffa да, я стремился достичь почти 1 мс, но микро это другая история. Для микро я рассмотрел использование инкрементной обратной структуры в пакете, обнаруживая подобные матрицы, которые встречаются часто. Perftrely - это диапазон 10-500 мс, в зависимости от процессора. Точность двойная или даже сложная двойная. Одиночная точность делает медленнее.
rfabbri
@njuffa Я могу уменьшить или повысить точность для скорости
rfabbri
2
Кажется, что точность / точность не ваш приоритет. Для вашей цели, может быть, полезен итерационный метод, усеченный при сравнительно небольшом числе оценок? Особенно, если у вас есть разумное начальное предположение.
Спенсер Брингельсон
1
Вы разворачиваетесь? Не могли бы вы сделать QR-факторизацию вместо исключения Гаусса? Вы чередуете свои системы, так что вы можете использовать инструкции SIMD и работать с несколькими системами одновременно? Вы пишете прямые программы без циклов и косвенной адресации? Какую точность вы хотите, и как я буду определять вашу систему? Есть ли у них какая-либо структура, которая может быть использована.
Карл Кристиан

Ответы:

7

Использование собственного матричного типа, в котором число строк и столбцов кодируется в тип во время компиляции, дает преимущество над LAPACK, где размер матрицы известен только во время выполнения. Эта дополнительная информация позволяет компилятору выполнять полное или частичное развертывание цикла, устраняя множество инструкций ветвления. Если вы предпочитаете использовать существующую библиотеку, а не писать собственные ядра, возможно, важно иметь тип данных, в котором размер матрицы может быть включен в качестве параметров шаблона C ++. Единственная из известных мне библиотек, которая делает это, - пламя , так что, возможно, стоит сравнить с Эйгеном.

Если вы решите развернуть собственную реализацию, вы можете найти то, что PETSc делает для своего блочного CSR-формата, как полезный пример, хотя сам PETSc, вероятно, не будет подходящим инструментом для ваших целей. Вместо того, чтобы писать цикл, они записывают каждую отдельную операцию для небольших умножений матрицы на вектор (см. Этот файл в их хранилище). Это гарантирует, что нет никаких ветвящихся инструкций, которые вы могли бы получить с помощью цикла. Версии кода с инструкциями AVX являются хорошим примером того, как на самом деле использовать векторные расширения. Например, эта функция использует__m256dТип данных для одновременной работы на четырех дублях. Вы можете получить заметное повышение производительности, явно записав все операции с использованием векторных расширений, только для факторизации LU вместо умножения матрицы на вектор. Вместо того, чтобы на самом деле писать код на C вручную, вам лучше использовать скрипт для его генерации. Также было бы интересно посмотреть, есть ли заметная разница в производительности, когда вы переупорядочиваете некоторые операции, чтобы лучше использовать конвейерную обработку команд.

Вы также можете получить некоторое преимущество от инструмента STOKE , который будет случайным образом исследовать пространство возможных преобразований программы, чтобы найти более быструю версию.

Даниэль Шаперо
источник
ТХ. Я уже использую Eigen как Map <const Matrix <complex, 10, 10>> AA (A) успешно. проверим в другие вещи.
rfabbri
Eigen также имеет AVX и даже заголовок complex.h для него. Почему PETSc для этого? В этом случае сложно конкурировать с Эйгеном. Я еще больше специализировал Eigen для своей задачи и применил приблизительную стратегию разворота, которая вместо того, чтобы брать максимум по столбцу, меняет своп немедленно, когда он находит другой, который на 3 порядка больше.
rfabbri
1
@ rfabbri Я не предлагал вам использовать PETSc для этого, только то, что они делают в данном конкретном случае, может быть поучительным. Я отредактировал ответ, чтобы прояснить ситуацию.
Даниэль Шаперо
4

Другая идея может заключаться в использовании генеративного подхода (программа, пишущая программу). Автор (мета) программы, которая выплевывает последовательность инструкций C / C ++ для выполнения неотключенного ** LU в системе 10x10 .. в основном беря гнездо цикла k / i / j и выравнивая его в O (1000) или около того строк скалярная арифметика. Затем введите эту сгенерированную программу в любой оптимизирующий компилятор. Здесь я думаю, что это довольно интересно, так как удаление циклов раскрывает каждую зависимость от данных и избыточное подвыражение и дает компилятору максимальную возможность переупорядочить инструкции так, чтобы они хорошо отображались в реальном оборудовании (например, число исполнительных блоков, опасностей / остановок, поэтому на).

Если вы знаете все матрицы (или даже несколько из них), вы можете улучшить пропускную способность, вызывая встроенные функции / функции SIMD (SSE / AVX) вместо скалярного кода. Здесь вы будете использовать смущающий параллелизм между экземплярами, вместо того, чтобы гоняться за каким-либо параллелизмом внутри одного экземпляра. Например, вы можете выполнить 4 LU двойной точности одновременно, используя встроенные функции AVX256, упаковав 4 матрицы «поперек» регистра и выполнив одинаковые операции ** на всех из них.

** Отсюда акцент на непивотных LU. Поворот портит этот подход двумя способами. Во-первых, он вводит ответвления из-за выбора центра, что означает, что ваши зависимости от данных не так хорошо известны. Во-вторых, это означает, что разные «слоты» SIMD должны будут выполнять разные действия, потому что экземпляр A может вращаться не так, как экземпляр B. Поэтому, если вы решите что-то из этого, я бы предложил статически повернуть ваши матрицы перед вычислением (переставить наибольшую запись каждого столбца по диагонали).

rchilton1980
источник
Поскольку матрицы настолько малы, возможно, поворот можно убрать, если они предварительно масштабированы. Даже без предварительного поворота матриц. Все, что нам нужно, это то, что записи находятся в пределах 2-3 порядков друг от друга.
rfabbri
2

Ваш вопрос приводит к двум различным соображениям.

Во-первых, вам нужно выбрать правильный алгоритм. Следовательно, вопрос о том, имеют ли матрицы какую-либо структуру, должен быть рассмотрен. Например, когда матрицы симметричны, разложение Холецкого более эффективно, чем LU. Когда вам нужна только ограниченная точность, итерационный метод может быть быстрее.

Во-вторых, вам необходимо эффективно реализовать алгоритм. Для этого вам нужно знать узкое место вашего алгоритма. Ваша реализация ограничена скоростью передачи памяти или скоростью вычислений. Поскольку вы считаете только10×10матрицы, ваша матрица должна полностью вписаться в кэш процессора. Таким образом, вы должны использовать блоки SIMD (SSE, AVX и т. Д.) И ядра вашего процессора, чтобы выполнять как можно больше вычислений за цикл.

В целом, ответ на ваш вопрос сильно зависит от оборудования и матриц, которые вы рассматриваете. Вероятно, нет однозначного ответа, и вам нужно попробовать несколько вещей, чтобы найти оптимальный метод.

Х. Риттих
источник
Пока что Eigen сильно оптимизирует, использует SEE, AVX и т. Д., И я пробовал итерационные методы в предварительном тесте, но они не помогли. Я попробовал Intel MKL, но не лучше Eigen с оптимизированными флагами GCC. В настоящее время я пытаюсь сделать что-то лучше и проще, чем Эйген, и провести более подробные тесты с итерационными методами.
rfabbri
1

Я бы попробовал блочную инверсию.

https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

Eigen использует оптимизированную процедуру для вычисления обратной матрицы 4x4, что, вероятно, является лучшим, что вы получите. Попробуйте использовать это как можно больше.

http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html

Слева вверху: 8х8. Вверху справа: 8х2. Внизу слева: 2x8. Справа внизу: 2x2. Инвертируйте 8x8, используя оптимизированный код инверсии 4x4. Остальное - это матричные продукты.

РЕДАКТИРОВАТЬ: Использование блоков 6x6, 6x4, 4x6 и 4x4 оказалось немного быстрее, чем я описал выше.

using namespace Eigen;

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> blockwise_inversion(const Matrix<Scalar, tl_size, tl_size>& A, const Matrix<Scalar, tl_size, br_size>& B, const Matrix<Scalar, br_size, tl_size>& C, const Matrix<Scalar, br_size, br_size>& D)
{
    Matrix<Scalar, tl_size + br_size, tl_size + br_size> result;

    Matrix<Scalar, tl_size, tl_size> A_inv = A.inverse().eval();
    Matrix<Scalar, br_size, br_size> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<tl_size, tl_size>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<tl_size, br_size>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<br_size, tl_size>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<br_size, br_size>() = DCAB_inv;

    return result;
}

template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> my_inverse(const Matrix<Scalar, tl_size + br_size, tl_size + br_size>& mat)
{
    const Matrix<Scalar, tl_size, tl_size>& A = mat.topLeftCorner<tl_size, tl_size>();
    const Matrix<Scalar, tl_size, br_size>& B = mat.topRightCorner<tl_size, br_size>();
    const Matrix<Scalar, br_size, tl_size>& C = mat.bottomLeftCorner<br_size, tl_size>();
    const Matrix<Scalar, br_size, br_size>& D = mat.bottomRightCorner<br_size, br_size>();

    return blockwise_inversion<Scalar,tl_size,br_size>(A, B, C, D);
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_8_2(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 8, 8>& A = input.topLeftCorner<8, 8>();
    const Matrix<Scalar, 8, 2>& B = input.topRightCorner<8, 2>();
    const Matrix<Scalar, 2, 8>& C = input.bottomLeftCorner<2, 8>();
    const Matrix<Scalar, 2, 2>& D = input.bottomRightCorner<2, 2>();

    Matrix<Scalar, 8, 8> A_inv = my_inverse<Scalar, 4, 4>(A);
    Matrix<Scalar, 2, 2> DCAB_inv = (D - C * A_inv * B).inverse();

    result.topLeftCorner<8, 8>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<8, 2>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<2, 8>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<2, 2>() = DCAB_inv;

    return result;
}

template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_6_4(const Matrix<Scalar, 10, 10>& input)
{
    Matrix<Scalar, 10, 10> result;

    const Matrix<Scalar, 6, 6>& A = input.topLeftCorner<6, 6>();
    const Matrix<Scalar, 6, 4>& B = input.topRightCorner<6, 4>();
    const Matrix<Scalar, 4, 6>& C = input.bottomLeftCorner<4, 6>();
    const Matrix<Scalar, 4, 4>& D = input.bottomRightCorner<4, 4>();

    Matrix<Scalar, 6, 6> A_inv = my_inverse<Scalar, 4, 2>(A);
    Matrix<Scalar, 4, 4> DCAB_inv = (D - C * A_inv * B).inverse().eval();

    result.topLeftCorner<6, 6>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
    result.topRightCorner<6, 4>() = -A_inv * B * DCAB_inv;
    result.bottomLeftCorner<4, 6>() = -DCAB_inv * C * A_inv;
    result.bottomRightCorner<4, 4>() = DCAB_inv;

    return result;
}

Вот результаты одного теста с использованием миллиона Eigen::Matrix<double,10,10>::Random()матриц и Eigen::Matrix<double,10,1>::Random()векторов. Во всех моих тестах обратное всегда быстрее. Моя процедура решения включает в себя вычисление обратного и затем умножение его на вектор. Иногда это быстрее, чем Эйген, иногда нет. Мой метод разметки может быть ошибочным (не отключал турбо-буст и т. Д.). Кроме того, случайные функции Эйгена могут не представлять реальные данные.

  • Собственный частичный обратный поворот: 3036 миллисекунд
  • Мой обратный с верхним блоком 8x8: 1638 миллисекунд
  • Мой инверс с верхним блоком 6x6: 1234 миллисекунды
  • Собственное частичное разворотное решение: 1791 миллисекунда
  • Мое решение с верхним блоком 8x8: 1739 миллисекунд
  • Мое решение с верхним блоком 6x6: 1286 миллисекунд

Мне очень интересно посмотреть, сможет ли кто-нибудь еще оптимизировать это, так как у меня есть приложение с конечными элементами, которое инвертирует матрицу gazillion 10x10 (и да, мне нужны индивидуальные коэффициенты обратного преобразования, так что непосредственное решение линейной системы не всегда возможно) ,

Чарли С
источник