Эксплуатация симметрии в броненосце - PullRequest
Новая
       12

Эксплуатация симметрии в броненосце

0 голосов
/ 28 января 2019

Есть ли способ использовать тот факт, что результирующая матрица в броненосце симметрична?У меня часто есть продукты вида A * B * A.t() (где B симметрична и положительно определена), и мне интересно, есть ли какой-нибудь способ сделать умножение быстрее, используя этот факт.Я сделал несколько простых попыток, но они либо одинаково быстры, либо существенно медленнее.Было бы интересно увидеть функцию, которая обеспечивает улучшение скорости (или объясняет, почему ее невозможно получить). 

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
arma::mat m1(const arma::mat & A, const arma::mat & B) {
  arma::mat C = (A * B) * A.t();
  return C;
}

// [[Rcpp::export]]
arma::mat m2(const arma::mat & A, const arma::mat & B) {
  arma::mat C1 = (A * B);
  arma::mat C2 = arma::mat(arma::size(C1), arma::fill::zeros);
  arma::uword n = C1.n_rows;

  for (arma::uword i = 0; i < n; i++) {
    for (arma::uword j = 0; j <= i; j++) {
      C2(i, j) = arma::as_scalar(C1.row(i) * A.row(j).t());
    }
  }
  C2 = arma::symmatl(C2);
  return C2;
}

// [[Rcpp::export]]
arma::mat m3(const arma::mat & A, const arma::mat & B) {
  arma::mat C1 = (A * B);
  arma::mat C2 = arma::mat(arma::size(C1), arma::fill::zeros);
  arma::uword n = C1.n_rows;

  for (arma::uword i = 0; i < n; i++) {
    for (arma::uword j = 0; j <= i; j++) {
      C2.submat(i, j, i, j) = C1.row(i) * A.row(j).t();
    }
  }
  C2 = arma::symmatl(C2);
  return C2;
}

/*** R
n <- 10
A <- matrix(rnorm(n*n), n, n)
B <- crossprod(matrix(rnorm(n*n), n, n))


microbenchmark::microbenchmark(m1(A, B), m2(A, B), m3(A, B))

all.equal(m1(A, B), m2(A, B))
all.equal(m1(A, B), m3(A, B))
*/
...