R 中稀疏矩形矩阵的 LU 分解

Question

MATLAB 能够使用

[L, U, P, Q] = lu(A)

对稀疏矩形矩阵执行 LU 分解，但目前还没有 R 包可以执行此操作。 尝试在 R 中的稀疏矩阵上使用

Matrix::lu()

会返回错误，指出矩阵必须是方阵。

R 中是否存在与 MATLAB 对非满秩矩形矩阵进行 LU 分解类似的情况？澄清一下，这与内存问题无关 - 如果我将它嵌入到单位矩阵中（请参阅这个问题），那么

Matrix::lu()

会抱怨矩阵是奇异的，而 MATLAB 的

lu()

会愉快地继续。

我试图解决的问题是在 R 中实现一个需要 Moore-Penrose 伪逆的特定计量经济学估计器。

pracma

的 Moore-Penrose 伪逆因矩阵太大而失败。 显然，获得大型稀疏矩阵的 Moore-Penrose 伪逆也未解决，请参见here。 我的矩阵在任何给定的行和列中都有多个 1。

Answer 1

解决了。以下代码使用 UMFPACK，它与 MATLAB 使用的库相同。为了获得正确的元素，我们需要添加

Control[UMFPACK_SCALE] = 0;

。然后我们需要修复 R 中生成的 U 矩阵的大小。

#include <suitesparse/umfpack.h>
#include <Rcpp.h>
#include <algorithm>
#include <vector>
#include <iostream>

// [[Rcpp::export]]
Rcpp::List sparseLU(const std::vector<int> &Ap, const std::vector<int> &Ai, const std::vector<double> &Ax) {
  int n = Ap.size() - 1;
  int m = *std::max_element(Ai.begin(), Ai.end()) + 1;
  int nnz = Ax.size();
  double Control [UMFPACK_CONTROL], Info [UMFPACK_INFO];
  umfpack_di_defaults(Control);
  Control[UMFPACK_PRL] = 2;
  Control[UMFPACK_SCALE] = 0;
  Control[UMFPACK_STRATEGY] = UMFPACK_STRATEGY_SYMMETRIC;

  void *symbolic, *numeric;
  int status;
  status = umfpack_di_symbolic(m, n, Ap.data(), Ai.data(), Ax.data(), &symbolic, Control, Info);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      Rcpp::stop("umfpack_di_symbolic failed with status %d", status);
    }
  }

  status = umfpack_di_numeric(Ap.data(), Ai.data(), Ax.data(), symbolic, &numeric, Control, Info);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      Rcpp::stop("umfpack_di_numeric failed with status %d", status);
      umfpack_di_free_symbolic(&symbolic);
    }
  }

  umfpack_di_free_symbolic(&symbolic);

  int lnz, unz, nz_udiag;
  status = umfpack_di_get_lunz(&lnz, &unz, &m, &n, &nz_udiag, numeric);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      Rcpp::stop("umfpack_di_get_lunz failed with status %d", status);
      umfpack_di_free_numeric(&numeric);
    }
  }


  std::vector<int> Lp(m+1), Lj(lnz), Up(n+1), Ui(unz), P(m), Q(n);
  std::vector<double> Lx(lnz), Ux(unz), D(std::min(m, n)), Rs(m);
  int do_recip;

  status = umfpack_di_get_numeric(Lp.data(), Lj.data(), Lx.data(), Up.data(), Ui.data(), Ux.data(), P.data(), Q.data(), D.data(), &do_recip, Rs.data(), numeric);
  if (status != UMFPACK_OK) {
    umfpack_di_report_status(Control, status);
    if (status < 0) {
      umfpack_di_free_numeric(&numeric);
      Rcpp::stop("umfpack_di_get_numeric failed with status %d", status);
    }
  }

  umfpack_di_free_numeric(&numeric);

  Rcpp::List ret = Rcpp::List::create(
      Rcpp::Named("L") = Rcpp::List::create(Rcpp::Named("j") = Lj, Rcpp::Named("p") = Lp, Rcpp::Named("x") = Lx),
      Rcpp::Named("U") = Rcpp::List::create(Rcpp::Named("i") = Ui, Rcpp::Named("p") = Up, Rcpp::Named("x") = Ux),
      Rcpp::Named("P") = P,
      Rcpp::Named("Q") = Q
      );


  return ret;
}

然后我们可以像这样在R中加载它

# A is a sparseMatrix() from package Matrix
LUPQ <- sparseLU(A@p, A@i, A@x) 
n <- ncol(A)
L <- sparseMatrix(j = LUPQ$L$j, p = LUPQ$L$p, x = LUPQ$L$x, index1 = FALSE) 
U <- sparseMatrix(i = LUPQ$U$i, p = LUPQ$U$p, x = LUPQ$U$x, index1 = FALSE, dims = c(n, n))
Q <- sparseMatrix(i = LUPQ$Q + 1, j = 1:length(LUPQ$Q), x = 1)