Быстрое решение нижней треугольной системы с использованием Rcpp Eigen и Cholesky - PullRequest
0 голосов
/ 08 июня 2019

У меня есть следующая проблема, я пытаюсь найти быстрый способ кодирования.Это в основном эта проблема (https://math.stackexchange.com/questions/2304623/how-to-easily-obtain-the-quadratic-form-of-an-inverse), где я хочу вычислить большое количество решений разреженной обратной треугольной системы, используя Rcpp и Eigen.

Моя первая попытка некоторого кода приведена ниже.более быстрый способ, чем стратегия solveInPlace, на которую я сейчас полагаюсь?

Eigen::VectorXd fast_solution(
  const Eigen::MappedSparseMatrix<double> cholL, // Lower Triangular Cholesky.
  const Eigen::MappedSparseMatrix<double> Z
){
  int size_Z = Z.cols();

  Eigen::VectorXd output(size_Z);
  for (int i = 0; i < size_Z; i++){
    Eigen::SparseVector<double> v_i = Z.col(i);
    cholL.triangularView<Eigen::Lower>().solveInPlace(v_i);
    output[i] = v_i.squaredNorm();
  }
  return output;
}
...