Блог

Главная - Вопросы по программированию - Существует ли оптимизированный способ вычисления `x ^ T A X` для симметричного `A` в собственном файле ?

Существует ли оптимизированный способ вычисления `x ^ T A X` для симметричного `A` в собственном файле ?

#matrix #vector #linear-algebra #numerical-methods #eigen3

#матрица #вектор #линейная алгебра #численные методы #eigen3

Вопрос:

Учитывая симметричную матрицу A и вектор x , мне часто нужно вычислять x^T * A * x . Я могу сделать это с помощью Eigen 3.x x.transpose() * (A * x) , но это не использует информацию, которая x одинакова с обеих сторон и A симметрична. Есть ли более эффективный способ вычислить это?

Ответ №1:

Как часто вы это вычисляете? Если очень часто для разных x , то это может немного ускорить вычисление разложения Холецкого или LDLT A и использовать, что произведение треугольной матрицы с вектором требует только половины умножений.

Или, возможно, еще проще, если вы разложите A=L D L.T , где L строго нижний треугольный и D диагональный, тогда

 x.T*A*x = x.T*D*x   2*x.T*(L*x)

где первые члены — это сумма по d[k]*x[k]**2 . Второй член, при тщательном использовании треугольной структуры, использует половину операций исходного выражения.

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

Комментарии:

1. Мой вариант использования заключается в том, что A и x меняются вместе.

2. Было бы это полезно? eigen.tuxfamily.org/dox/classEigen_1_1TriangularView.html

3. Я думаю, что способ сообщить Eigen «эй, эта матрица симметрична» заключается в использовании этого: eigen.tuxfamily.org/dox/classEigen_1_1SelfAdjointView.html

Ответ №2:

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

 double xAx_symmetric(const Eigen::MatrixXdamp; A, const Eigen::VectorXdamp; x)
{
    const auto dim = A.rows();
    if (dim < 15) {
        double sum = 0;
        for (Eigen::Index i = 0; i < dim;   i) {
            const auto x_i = x[i];
            sum  = A(i, i) * x_i * x_i;
            for (Eigen::Index j = 0; j < i;   j) {
                sum  = 2 * A(j, i) * x_i * x[j];
            }
        }
        return sum;
    }
    else {
        return x.transpose() * A.selfadjointView<Eigen::Upper>() * x;
    }
}