我正在尝试从 R 中的下三角矩阵创建对称矩阵。
在之前的问答中 ( Convert upper triangular part of a matrix to symmetric matrix ) 用户 李哲源说对于大型矩阵，这不应该在 R 中完成，并在 C 中提出了一个解决方案。但我不了解 C 并且之前从未使用过例如 Rccp 所以不知道如何解释答案。然而很明显，那里的 C 代码生成了我不想要的随机数 (rnorm)。 我想输入一个方阵，得到一个对称矩阵。
对于给定的方阵 A，其下三角中有数据，我如何在 C 中有效地创建对称矩阵并在 R 中使用它？
最佳答案
快速改编自 as.matrix on a distance object is extremely slow; how to make it faster? 中的 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 方法不是最优的，原因有以下三个:
- 它需要将
x槽作为数字 vector 传递，因此我们必须执行base::c(A)，这实际上是在 RAM 中创建矩阵的副本； - 不能原地修改，所以创建一个新的矩阵副本作为输出矩阵；
- 转置复制时不阻塞缓存。
这是一个快速比较:
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 有多快。在“覆盖”模式下也没有内存分配。
As G. Grothendieck mentions ，我们也可以使用“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”的效率。
关于c - 使三角形(或通常为正方形)矩阵对称，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/52571101/