Sum of the all of the outer product (for Rcpp)

sumouterprod2.cpp

#include "RcppArmadillo.h"
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
double sumouterprod(const arma::field<arma::vec> & V){
  int K = V.n_rows;
  arma::vec tout = V(0);
  for(int j=1; j<K; j++){
    arma::vec V1 = V(j);
    tout = vectorise(tout * V1.t());
  }
  return sum(tout);
}

double sumouter(arma::vec & x, arma::vec & y){
  double res = 0;
  for(int i = 0; i < x.n_rows; i++) {
    res += sum(x(i) * y);
  }
  return res;
}

// [[Rcpp::export]]
double sumouterprod2(const arma::field<arma::vec> & V){
  int K = V.n_rows;
  arma::vec tout = V(0);
  for(int j=1; j<K-1; j++){
    arma::vec V1 = V(j);
    tout = vectorise(tout * V1.t());
  }
  double out;
  if(K>1){
    arma::vec VK =V(K-1);
    out = sumouter(tout, VK);
  }else{
    out = sum(tout);
  }
  return out;
}