создание симметричной треугольной (или вообще квадратной) матрицы - PullRequest
Новая
       4

создание симметричной треугольной (или вообще квадратной) матрицы

0 голосов
/ 29 сентября 2018

Я пытаюсь создать симметричную матрицу из нижней треугольной матрицы в R.

В предыдущем Q & A ( Преобразовать верхнюю треугольную часть матрицы в симметричную матрицу ) пользователя 李哲源 сказал, что для больших матриц это не следует делать в R, и предложил решение в C. Но я не понимаю C и никогда раньше не использовал, например, Rccp, поэтому не знаю, какинтерпретировать ответ.Все же очевидно, что код C там генерирует случайные числа (rnorm), которые я не хочу. Я хочу поместить квадратную матрицу и получить симметричную матрицу.

Для данной квадратной матрицы A, в которой данные находятся в нижнем треугольнике, как бы я создал симметричную матрицу?эффективно в C и использовать его в R?

Ответы [ 2 ]

0 голосов
/ 30 сентября 2018

Перед выполнением, возможно, преждевременной оптимизации с помощью C / C ++, проверьте, достаточна ли плотная матрица 

A + t(A)

(при условии, что только ненулевые элементы А находятся ниже диагонали или выше нее.

Также,если проблема с памятью, то пакет Matrix имеет упакованный симметричный класс dspMatrix, который можно создать следующим образом: 

library(Matrix)

A <- matrix(c(0, 2, 3, 0, 0, 4, 0, 0, 0), 3) # dense lower triangular test input
dspA <- as(A + t(A), "dspMatrix")

, что дает: 

> str(dspA)
Formal class 'dspMatrix' [package "Matrix"] with 5 slots
  ..@ x       : num [1:6] 0 2 0 3 4 0   <- only 6 elements stored, not 9
  ..@ Dim     : int [1:2] 3 3
  ..@ Dimnames:List of 2
  .. ..$ : NULL
  .. ..$ : NULL
  ..@ uplo    : chr "U"
  ..@ factors : list()

или его можно создать напрямуюот верхней треугольной части: 

# use upper triangular part since we already created dspA that way
tA <- t(A)
dspA2 <- new("dspMatrix", Dim = as.integer(c(3,3)), 
  x = tA[upper.tri(tA, diag = TRUE)])

identical(dspA, dspA2)
## [1] TRUE
0 голосов
/ 29 сентября 2018

Быстро адаптируется из функции dist2mat в , поскольку as.matrix на объекте расстояния очень медленный;как сделать это быстрее? . 

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), который по существу создает копию матрицы в ОЗУ;
  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.Кроме того, в режиме «перезаписи» не производится выделение памяти.

Как упоминает Дж. Гротендик , мы также можем использовать «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

Снова *Метод 1065 * является неоптимальным из-за использования 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".

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.

Нет похожих вопросов

...