如何找到所有反对角线的总和?

How to find the sum of all anti-diagonals?

我有一个矩阵 M:

n = 3    
x=c(0.85, 0.1, 0.05)
M <- matrix(NA, n, n); 

for(i in 1:n){
for(j in 1:n){
M[i,j] = x[i] * x[j]
}}

#       [,1]  [,2]   [,3]
# [1,] 0.7225 0.085 0.0425
# [2,] 0.0850 0.010 0.0050
# [3,] 0.0425 0.005 0.0025

我需要求所有反对角线的和,包括 M[1,1] 和 M[n, n]。我的尝试是

d <-matrix(c(0, 1, 2, 1, 2, 3, 2, 3, 4), n)
tapply(M, d, sum)

     0      1      2      3      4 
0.7225 0.1700 0.0950 0.0100 0.0025 

结果对我来说是正确的。

问题。如何定义矩阵d的元素?可以作为 col(M) 和 row(M) 的函数。

你可以这样做:

sapply(seq(3), function(x) seq(3) + x - 2)
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    2    3
#> [3,]    2    3    4

或者更一般地说,

anti_diagonal <- function(n) sapply(seq(n), function(x) seq(n) + x - 2)

例如:

anti_diagonal(6)
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    0    1    2    3    4    5
#> [2,]    1    2    3    4    5    6
#> [3,]    2    3    4    5    6    7
#> [4,]    3    4    5    6    7    8
#> [5,]    4    5    6    7    8    9
#> [6,]    5    6    7    8    9   10

首先请注意 outer 可以生成矩阵 d 而无需明确列出其元素。

matrix(c(0, 1, 2, 1, 2, 3, 2, 3, 4), 3)
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    2    3
#> [3,]    2    3    4
outer(0:2, 0:2, `+`)
#>      [,1] [,2] [,3]
#> [1,]    0    1    2
#> [2,]    1    2    3
#> [3,]    2    3    4

reprex package (v2.0.1)

于 2022-03-24 创建

并在函数中使用它。

sumAntiDiag <- function(M){
  nr <- nrow(M)
  nc <- ncol(M)
  d <- outer(seq.int(nr), seq.int(nc), `+`)
  tapply(M, d, sum)
}

n <- 3    
x <- c(0.85, 0.1, 0.05)
M <- matrix(NA, n, n); 

for(i in 1:n){
  for(j in 1:n){
    M[i,j] = x[i] * x[j]
  }}

sumAntiDiag(M)
#>      2      3      4      5      6 
#> 0.7225 0.1700 0.0950 0.0100 0.0025

reprex package (v2.0.1)

于 2022-03-24 创建

您可以使用 sequence:

function(n) matrix(sequence(rep(n, n), seq(n) - 1), nrow = n)

输出

f <- function(n) matrix(sequence(rep(n, n), seq(n) - 1), nrow = n)
f(3)
     [,1] [,2] [,3]
[1,]    0    1    2
[2,]    1    2    3
[3,]    2    3    4

f(5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    2    3    4
[2,]    1    2    3    4    5
[3,]    2    3    4    5    6
[4,]    3    4    5    6    7
[5,]    4    5    6    7    8

通过使用基数 R 中的 embed 定义函数 f 来尝试下面的代码,即

f <- function(n) embed(seq(2 * n - 1) - 1, n)[, n:1]

这样

> f(3)
     [,1] [,2] [,3]
[1,]    0    1    2
[2,]    1    2    3
[3,]    2    3    4

> f(4)
     [,1] [,2] [,3] [,4]
[1,]    0    1    2    3
[2,]    1    2    3    4
[3,]    2    3    4    5
[4,]    3    4    5    6

> f(5)
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    1    2    3    4
[2,]    1    2    3    4    5
[3,]    2    3    4    5    6
[4,]    3    4    5    6    7
[5,]    4    5    6    7    8

正如您在问题中提到的,可以使用 row(M)col(M),尽管它们 rows/columns 从 1 而不是零开始,因此您需要减去 2(1每个)给予:

tapply(M, row(M) + col(M) - 2, sum)
#     0      1      2      3      4 
#0.7225 0.1700 0.0950 0.0100 0.0025

使用索引代替 tapply 会加快速度。或者 Rcpp:

sumdiags <- function(mat, minor = TRUE) {
  m <- ncol(mat)
  
  if (minor) {
    n <- nrow(mat)
    lens <- c(pmin(1:n, m), pmin((m - 1L):1, n))
    c(mat[1], diff(cumsum(mat[sequence(lens, c(1:n, seq(2L*n, by = n, length.out = m - 1L)), n - 1L)])[cumsum(lens)]))
  } else {
    Recall(mat[,m:1])
  }
}

# compare to tapply solution
sumdiags2 <- function(mat, minor = TRUE) {
  if (minor) {
    as.numeric(tapply(mat, row(mat) + col(mat), sum))
  } else {
    Recall(mat[,ncol(mat):1])
  }
}

# or Rcpp
Rcpp::cppFunction('NumericVector sumdiagsRcpp(const NumericMatrix& mat) {
  const int n = mat.nrow();
  const int m = mat.ncol();
  NumericVector x (n + m - 1);

  for(int row = 0; row < n; row++) {
    for(int col = 0; col < m; col++) {
      x[row + col] += mat(row, col);
    }
  }
  return x;
}')

# OP data
x <- c(0.85, 0.1, 0.05)
m <- outer(x, x)
sumdiags(m)
#> [1] 0.7225 0.1700 0.0950 0.0100 0.0025
sumdiags2(m)
#> [1] 0.7225 0.1700 0.0950 0.0100 0.0025
sumdiagsRcpp(m)
#> [1] 0.7225 0.1700 0.0950 0.0100 0.0025

# bigger matrix for benchmarking
m <- matrix(runif(1e6), 1e3)
microbenchmark::microbenchmark(sumdiags = sumdiags(m),
                               sumdiags2 = sumdiags2(m),
                               sumdiagsRcpp = sumdiagsRcpp(m),
                               check = "equal")
#> Unit: milliseconds
#>         expr       min        lq      mean    median      uq        max neval
#>     sumdiags  9.985302 10.266350 13.686723 10.803401 17.5274  22.387601   100
#>    sumdiags2 55.790402 65.140051 78.763478 67.120051 70.4165 183.936801   100
#> sumdiagsRcpp  2.192201  2.378651  2.599326  2.631751  2.7050   4.038301   100