使三角形(或通常是正方形)矩阵对称

Question

我试图从R中的下三角矩阵创建一个对称矩阵.

在之前的Q&A(将矩阵的上三角部分转换为对称矩阵)中,用户李哲源表示,对于大型矩阵,这不应该在R中完成并且在C中提出了解决方案.但是我不理解C并且从来没有例如Rccp之前用过,所以不知道如何解释答案.然而很明显,那里的C代码会生成rnorm我不想要的随机数().我想放入一个方阵并得出一个对称矩阵.

对于A在其下三角形中具有数据的给定方阵,如何在C中有效地创建对称矩阵并在R中使用它？

Answer 1

从as.matrix中的dist2mat函数快速改编为距离对象非常慢; 如何让它更快？.

library(Rcpp)

cppFunction('NumericMatrix Mat2Sym(NumericMatrix A, bool up2lo, int bf) {

  IntegerVector dim = A.attr("dim");
  size_t n = (size_t)dim[0], m = (size_t)dim[1];
  if (n != m) stop("A is not a square matrix!");

  /* use pointers */
  size_t j, i, jj, ni, nj;
  double *A_jj, *A_ij, *A_ji, *col, *row, *end;

  /* cache blocking factor */
  size_t b = (size_t)bf;

  /* copy lower triangular to upper triangular; cache blocking applied */
  for (j = 0; j < n; j += b) {
    nj = n - j; if (nj > b) nj = b;
    /* diagonal block has size nj x nj */
    A_jj = &A(j, j);
    for (jj = nj - 1; jj > 0; jj--, A_jj += n + 1) {
      /* copy a column segment to a row segment (or vise versa) */
      col = A_jj + 1; row = A_jj + n;
      for (end = col + jj; col < end; col++, row += n) {
        if (up2lo) *col = *row; else *row = *col;
        }
      }
    /* off-diagonal blocks */
    for (i = j + nj; i < n; i += b) {
      ni = n - i; if (ni > b) ni = b;
      /* off-diagonal block has size ni x nj */
      A_ij = &A(i, j); A_ji = &A(j, i);
      for (jj = 0; jj < nj; jj++) {
        /* copy a column segment to a row segment (or vise versa) */
        col = A_ij + jj * n; row = A_ji + jj;
        for (end = col + ni; col < end; col++, row += n) {
          if (up2lo) *col = *row; else *row = *col;
          }
        }
      }
    }

  return A;
  }')

对于方形矩阵A,此函数将Mat2Sym其下三角形部分(具有换位)复制到其上三角形部分以使其对称,如果up2lo = FALSE,则反之亦然up2lo = TRUE.

请注意,该函数会在不使用额外内存的情况下覆盖 A.要保留输入矩阵并创建新的输出矩阵,请A + 0不要A将其传递给函数.

## an arbitrary asymmetric square matrix
set.seed(0)
A <- matrix(runif(25), 5)
#          [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.8966972 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.2655087 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.3721239 0.9446753 0.17655675 0.7176185 0.2121425
#[4,] 0.5728534 0.6607978 0.68702285 0.9919061 0.6516738
#[5,] 0.9082078 0.6291140 0.38410372 0.3800352 0.1255551

## lower triangular to upper triangular; don't overwrite
B <- Mat2Sym(A + 0, up2lo = FALSE, 128)
#          [,1]      [,2]      [,3]      [,4]      [,5]
#[1,] 0.8966972 0.2655087 0.3721239 0.5728534 0.9082078
#[2,] 0.2655087 0.8983897 0.9446753 0.6607978 0.6291140
#[3,] 0.3721239 0.9446753 0.1765568 0.6870228 0.3841037
#[4,] 0.5728534 0.6607978 0.6870228 0.9919061 0.3800352
#[5,] 0.9082078 0.6291140 0.3841037 0.3800352 0.1255551

## A is unchanged
A
#          [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.8966972 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.2655087 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.3721239 0.9446753 0.17655675 0.7176185 0.2121425
#[4,] 0.5728534 0.6607978 0.68702285 0.9919061 0.6516738
#[5,] 0.9082078 0.6291140 0.38410372 0.3800352 0.1255551

## upper triangular to lower triangular; overwrite
D <- Mat2Sym(A, up2lo = TRUE, 128)
#           [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.89669720 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.20168193 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.06178627 0.2059746 0.17655675 0.7176185 0.2121425
#[4,] 0.76984142 0.4976992 0.71761851 0.9919061 0.6516738
#[5,] 0.77744522 0.9347052 0.21214252 0.6516738 0.1255551

## A has been changed; D and A are aliased in memory
A
#           [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.89669720 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.20168193 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.06178627 0.2059746 0.17655675 0.7176185 0.2121425
#[4,] 0.76984142 0.4976992 0.71761851 0.9919061 0.6516738
#[5,] 0.77744522 0.9347052 0.21214252 0.6516738 0.1255551

关于使用Matrix包

Matrix对于稀疏矩阵特别有用.为了兼容性,它还为密集矩阵定义了一些类,如"dgeMatrix","dtrMatrix","dtpMatrix","dsyMatrix","dspMatrix".

给定方形矩阵A,Matrix使其对称的方式如下.

set.seed(0)
A <- matrix(runif(25), 5)
#          [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.8966972 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.2655087 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.3721239 0.9446753 0.17655675 0.7176185 0.2121425
#[4,] 0.5728534 0.6607978 0.68702285 0.9919061 0.6516738
#[5,] 0.9082078 0.6291140 0.38410372 0.3800352 0.1255551

## equivalent to: Mat2Sym(A + 0, TRUE, 128)
new("dsyMatrix", x = base::c(A), Dim = dim(A), uplo = "U")
#5 x 5 Matrix of class "dsyMatrix"
#           [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.89669720 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.20168193 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.06178627 0.2059746 0.17655675 0.7176185 0.2121425
#[4,] 0.76984142 0.4976992 0.71761851 0.9919061 0.6516738
#[5,] 0.77744522 0.9347052 0.21214252 0.6516738 0.1255551

## equivalent to: Mat2Sym(A + 0, FALSE, 128)
new("dsyMatrix", x = base::c(A), Dim = dim(A), uplo = "L")
#5 x 5 Matrix of class "dsyMatrix"
#          [,1]      [,2]      [,3]      [,4]      [,5]
#[1,] 0.8966972 0.2655087 0.3721239 0.5728534 0.9082078
#[2,] 0.2655087 0.8983897 0.9446753 0.6607978 0.6291140
#[3,] 0.3721239 0.9446753 0.1765568 0.6870228 0.3841037
#[4,] 0.5728534 0.6607978 0.6870228 0.9919061 0.3800352
#[5,] 0.9082078 0.6291140 0.3841037 0.3800352 0.1255551

Matrix 方法是次优的,原因有三:

1. 它需要传递x槽作为数字向量,所以我们必须这样做base::c(A)基本上在RAM中创建矩阵的副本;
2. 它不能进行就地修改,因此创建一个新的矩阵副本作为输出矩阵;
3. 在进行转置复制时,它不会执行缓存阻止.

这是一个快速比较:

library(bench)
A <- matrix(runif(5000 * 5000), 5000)
bench::mark("Mat2Sym" = Mat2Sym(A, FALSE, 128),
            "Matrix" = new("dsyMatrix", x = base::c(A), Dim = dim(A), uplo = "L"),
            check = FALSE)
#  expression      min     mean   median     max `itr/sec` mem_alloc  n_gc n_itr
#  <chr>      <bch:tm> <bch:tm> <bch:tm> <bch:t>     <dbl> <bch:byt> <dbl> <int>
#1 Mat2Sym      56.8ms   57.7ms   57.4ms  59.4ms     17.3     2.48KB     0     9
#2 Matrix      334.3ms  337.4ms  337.4ms 340.6ms      2.96  190.74MB     2     2

注意有多快Mat2Sym.在"覆盖"模式下也没有进行内存分配.

正如G.格洛腾迪克提到的那样,我们也可以使用"dspMatrix".

set.seed(0)
A <- matrix(runif(25), 5)
#          [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.8966972 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.2655087 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.3721239 0.9446753 0.17655675 0.7176185 0.2121425
#[4,] 0.5728534 0.6607978 0.68702285 0.9919061 0.6516738
#[5,] 0.9082078 0.6291140 0.38410372 0.3800352 0.1255551

## equivalent to: Mat2Sym(A + 0, TRUE, 128)
new("dspMatrix", x = A[upper.tri(A, TRUE)], Dim = dim(A), uplo = "U")
#5 x 5 Matrix of class "dspMatrix"
#           [,1]      [,2]       [,3]      [,4]      [,5]
#[1,] 0.89669720 0.2016819 0.06178627 0.7698414 0.7774452
#[2,] 0.20168193 0.8983897 0.20597457 0.4976992 0.9347052
#[3,] 0.06178627 0.2059746 0.17655675 0.7176185 0.2121425
#[4,] 0.76984142 0.4976992 0.71761851 0.9919061 0.6516738
#[5,] 0.77744522 0.9347052 0.21214252 0.6516738 0.1255551

## equivalent to: Mat2Sym(A + 0, FALSE, 128)
new("dspMatrix", x = A[lower.tri(A, TRUE)], Dim = dim(A), uplo = "L")
#5 x 5 Matrix of class "dspMatrix"
#          [,1]      [,2]      [,3]      [,4]      [,5]
#[1,] 0.8966972 0.2655087 0.3721239 0.5728534 0.9082078
#[2,] 0.2655087 0.8983897 0.9446753 0.6607978 0.6291140
#[3,] 0.3721239 0.9446753 0.1765568 0.6870228 0.3841037
#[4,] 0.5728534 0.6607978 0.6870228 0.9919061 0.3800352
#[5,] 0.9082078 0.6291140 0.3841037 0.3800352 0.1255551

同样,Matrix由于使用upper.tri或,方法是次优的lower.tri.

library(bench)
A <- matrix(runif(5000 * 5000), 5000)
bench::mark("Mat2Sym" = Mat2Sym(A, FALSE, 128),
            "Matrix" = new("dspMatrix", x = A[lower.tri(A, TRUE)], Dim = dim(A),
                           uplo = "L"),
            check = FALSE)
#  expression      min     mean   median     max `itr/sec` mem_alloc  n_gc n_itr
#  <chr>      <bch:tm> <bch:tm> <bch:tm> <bch:t>     <dbl> <bch:byt> <dbl> <int>
#1 Mat2Sym      56.5ms   57.9ms     58ms  58.7ms     17.3     2.48KB     0     9
#2 Matrix      934.7ms  934.7ms    935ms 934.7ms      1.07  524.55MB     2     1

特别是,我们看到使用"dspMatrix"的效率甚至低于使用"dsyMatrix".