Почему посимвольное умножение менее эффективно, чем холистическое умножение для вычисления дополнения Шура? - PullRequest
Новая
       9

Почему посимвольное умножение менее эффективно, чем холистическое умножение для вычисления дополнения Шура?

0 голосов
/ 23 мая 2018

Мы пытаемся вычислить дополнение Шура для большой симметричной положительно определенной матрицы H, состоящей из диагональных блоков U и V и смещенных диагональных блоков W и W '(транспонирование W), а V - диагональная положительно определенная матрица.Формула дополнения Шура: 

H <-- U - W*V^(-1)*W'                       (1)

. При реализации операции одним способом является прямое использование вышеприведенного уравнения.

Во втором способе используется тот факт, что в V все нули отключены.диагональные элементы.Уравнение 1 записано и реализовано следующим образом: 

H <-- U - \Sigma W_i * V_{ii}^(-1) * (W_i)'      (2)

, где W_i - это i-й столбец W, V_ {ii} - это i-й диагональный элемент V, а суммы \ Sigma по каждому столбцу W. Очевидно,Ожидается, что (2) будет более эффективным, чем (1), поскольку он позволяет избежать умножения целых больших матриц.

Оба подхода реализованы с помощью библиотеки Eigen на языке C ++ и протестированы на ПК и мобильном устройстве сПроцессор Intel Atom.

Удивительно, но на 4-ядерном ПК Intel Core i7 фрагмент кода (скомпилированный с gcc, -O3, режим выпуска) для (2) менее эффективен, чем для (1),Среднее время вычисления для (1) и (2) составляет 0,092 мс, 0,212 мс, соответственно, обратите внимание на соотношение 2,30.На 4-ядерном процессоре Intel Atom наблюдается подобное (фрагмент кода, скомпилированный с помощью clang -O3).Среднее время вычислений для (1) и (2) составляет 0,575 мс, 1,923 мс соответственно, обратите внимание на соотношение 3.34.

Может ли кто-нибудь помочь мне объяснить, почему (2) медленнее, чем (1)?

Для вашего сведения ниже приведен фрагмент кода для обоих подходов. 

        // preconditioner
        Eigen::VectorXd p = (H_.diagonal().array() > 1.0e-9).select(H_.diagonal().cwiseSqrt(),1.0e-3);
        Eigen::VectorXd p_inv = p.cwiseInverse();

        // scale H and b
        H_ = p_inv.asDiagonal() * H_ * p_inv.asDiagonal();
        b0_ = p_inv.asDiagonal() * b0_;

        Eigen::MatrixXd U = H_.topLeftCorner(H_.rows() - marginalizationParametersLandmarks,
                                             H_.rows() - marginalizationParametersLandmarks);
        Eigen::MatrixXd V = H_.bottomRightCorner(marginalizationParametersLandmarks,
                                                 marginalizationParametersLandmarks);
        Eigen::MatrixXd W = H_.topRightCorner(H_.rows() - marginalizationParametersLandmarks,
                                              marginalizationParametersLandmarks);
        Eigen::VectorXd b_a = b0_.head(H_.rows() - marginalizationParametersLandmarks);
        Eigen::VectorXd b_b = b0_.tail(marginalizationParametersLandmarks);

        // split preconditioner
        Eigen::VectorXd p_a = p.head(H_.rows() - marginalizationParametersLandmarks);

        // invert the marginalization block
        static const int sdim = SIZE_PARAMETERIZATION_SIZE_FEATURE;
        b0_.resize(b_a.rows());
        b0_ = b_a;
        H_.resize(U.rows(), U.cols());
        H_ = U;

// surprisingly, the below operation is more efficient than its alternative on both
// an intel atom processor and on an intel core i7 PC
#if defined(ANDROID) || defined(__ANDROID__)
        Eigen::VectorXd V_inv(V.cols());
        for (int jack = 0; jack < V.cols(); ++jack) {
            V_inv[jack] = 1.0/V(jack, jack);
        }
        H_ -= W*V_inv.asDiagonal()*W.transpose();
        b0_ -= W*V_inv.asDiagonal()*b_b;
#else
        // assume V off-diagonal elements are all zeros
        Eigen::MatrixXd delta_H(U.rows(), U.cols());
        Eigen::VectorXd delta_b(b_a.rows());
        delta_H.setZero();
        delta_b.setZero();
        for (int i = 0; i < V.cols(); ++i) {
            double V_inv = 1.0 / V(i, i);
            double V_inv_sqrt = std::sqrt(V_inv);
            Eigen::VectorXd M = W.block(0, i, W.rows(), 1) * V_inv_sqrt;
            Eigen::VectorXd M1 = W.block(0, i, W.rows(), 1) * V_inv;
            delta_H += M * M.transpose();
            delta_b += M1 * b_b[i];
        }
        b0_ -= delta_b;
        H_ -= delta_H;
#endif
        // unscale
        H_ = p_a.asDiagonal() * H_ * p_a.asDiagonal();
        b0_ = p_a.asDiagonal() * b0_;
...