如何找到所有反对角线的总和?
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
我有一个矩阵 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