Nano Hash - криптовалюты, майнинг, программирование

Eigen Linear Solver для очень маленькой квадратной матрицы

Я использую Eigen в программе C++ для решения линейного уравнения для очень маленького квадрата. матрица (4X4).

Мой тестовый код похож на

template<template <typename MatrixType> typename EigenSolver>
Vertor3d solve(){
   //Solve Ax = b and A is a real symmetric matrix and positive semidefinite

   ... // Construct 4X4 square matrix A and 4X1 vector b

   EigenSolver<Matrix4d> solver(A);
   auto x = solver.solve(b);

   ... // Compute relative error for validating
}

Я тестирую некоторые EigenSolver, которые включают:

  1. FullPixLU
  2. PartialPivLU
  3. ДомохозяинQR
  4. ColPivHouseholderQR
  5. ColPivHouseholderQR
  6. Полная ортогональная декомпозиция
  7. LDLT
  8. Прямой обратный

Прямая инверсия это:

template<typename MatrixType>
struct InverseSolve
{
private:
    MatrixType  inv;
public:
    InverseSolve(const MatrixType &matrix) :inv(matrix.inverse()) {
    }
    template<typename VectorType>
    auto solve(const VectorType & b) {
        return inv * b;
    }
};

Я обнаружил, что быстрый метод - это DirectInverse. Даже если я связал Eigen с MKL, результат не изменился.

Это результат теста

  1. FullPixLU: 477 мс
  2. PartialPivLU: 468 мс
  3. HouseholderQR: 849 мс
  4. ColPivHouseholderQR: 766 мс
  5. ColPivHouseholderQR: 857 мс
  6. CompleteOrthogonalDecomposition: 832 мс
  7. LDLT : 477 мс
  8. Прямая инверсия: 88 мс

все они используют 1000000 матриц со случайным удвоением из равномерного распределения [0,100]. Сначала я строю верхний треугольник, а затем копирую в нижний треугольник.

Единственная проблема DirectInverse заключается в том, что его относительная ошибка немного больше, чем у других решателей, но приемлема.

Есть ли более быстрое или изящное решение для моей программы? Является ли DirectInverse быстрым решением для моей программы?

DirectInverse не использует симметричную информацию, так почему же DirectInverse намного быстрее, чем LDLT?

c++ eigen linear-algebra
18.06.2018

  • Вы вычисляете несколько solve на матрицу или только один? Пожалуйста, включите результат теста в виде текста в вопрос, а не в виде ссылки на изображение. Для матриц фиксированного размера до 4x4 inverse() вычисляется напрямую с использованием кофакторов, что очень быстро и требует только деления, но в целом ожидается, что оно будет немного менее точным. solve в том виде, в котором вы его реализовали, потребует только очень быстрого продукта Matrix-Vector. Вы можете повысить точность, выполнив шаг уточнения: x=inv*b; x+=inv*(b-A*x);. MKL обычно не помогает для небольших входных данных фиксированного размера. 18.06.2018
  • Спасибо за ваш ответ. Я бы отредактировал свою проблему с текстом эталона. В моей программе каждая матрица вычисляется и решается только один раз. Ваш совет очень конкретен и звучит как хорошая идея. Меня беспокоит то, что все мои матрицы симметричный, а inverse() его не использует. Должен ли я повторно реализовать inverse() для моей программы? 18.06.2018
  • Небольшое дополнение (независимо от основной причины): предлагаю назначить Vector4d x = solver.solve(b); вместо auto x = .... В противном случае вы могли бы x оценить несколько раз. 18.06.2018

Ответы:


1

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

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

Чтобы повысить точность (этот метод можно использовать независимо от решателя, если требуется), вы можете уточнить исходное решение, снова решив систему с невязкой:

Eigen::Vector4d solveDirect(const Eigen::Matrix4d& A, const Eigen::Vector4d& b)
{
    Eigen::Matrix4d inv = A.inverse();
    Eigen::Vector4d x = inv * b;
    x += inv*(b-A*x);
    return x;
}

Я не думаю, что Eigen напрямую предоставляет способ использовать симметрию A здесь (для непосредственно вычисляемого обратного). Вы можете попытаться намекнуть на это, явно скопировав самосопряженное представление A во временное и надеясь, что компилятор достаточно умен, чтобы найти общие подвыражения:

Eigen::Matrix4d tmp = A.selfadjointView<Eigen::Upper>();
Eigen::Matrix4d inv = tmp.inverse();

Чтобы уменьшить некоторые деления, вы также можете скомпилировать с помощью -freciprocal-math (на gcc или clang), это, конечно, немного снизит точность.

Если это действительно критично для производительности, попробуйте реализовать вручную настроенный метод inverse_4x4_symmetric. Использование симметрии inv * b вряд ли принесет пользу для таких маленьких матриц.

18.06.2018
  • Фактически, обратная матрица симметрии также является матрицей симметрии. selfadjointView<Eigen::Upper> также хорошо работает для inv. Я протестирую решение и сравню его с решением, настроенным вручную. Еще раз спасибо за вашу помощь. 19.06.2018
  • inv.selfadjointView<Eigen::Upper>() * b конечно работает, но скорее всего будет работать медленнее. Использование симметрии для матричных векторных произведений в основном имеет смысл, если это экономит место в кэше. 19.06.2018
    • Новые материалы
    Кластеризация: более глубокий взгляд
    Кластеризация — это метод обучения без учителя, в котором мы пытаемся найти группы в наборе данных на основе некоторых известных или неизвестных свойств, которые могут существовать. Независимо от..
    Как написать эффективное резюме
    Предложения по дизайну и макету, чтобы представить себя профессионально Вам не позвонили на собеседование после того, как вы несколько раз подали заявку на работу своей мечты? У вас может..
    Частный метод Python: улучшение инкапсуляции и безопасности
    Введение Python — универсальный и мощный язык программирования, известный своей простотой и удобством использования. Одной из ключевых особенностей, отличающих Python от других языков, является..
    Как я автоматизирую тестирование с помощью Jest
    Шутка для победы, когда дело касается автоматизации тестирования Одной очень важной частью разработки программного обеспечения является автоматизация тестирования, поскольку она создает..
    Работа с векторными символическими архитектурами, часть 4 (искусственный интеллект)
    Hyperseed: неконтролируемое обучение с векторными символическими архитектурами (arXiv) Автор: Евгений Осипов , Сачин Кахавала , Диланта Хапутантри , Тимал Кемпития , Дасвин Де Сильва ,..
    Понимание расстояния Вассерштейна: мощная метрика в машинном обучении
    В обширной области машинного обучения часто возникает необходимость сравнивать и измерять различия между распределениями вероятностей. Традиционные метрики расстояния, такие как евклидово..
    Обеспечение масштабируемости LLM: облачный анализ с помощью AWS Fargate и Copilot
    В динамичной области искусственного интеллекта все большее распространение получают модели больших языков (LLM). Они жизненно важны для различных приложений, таких как интеллектуальные..