2

我想做一个扭曲的矩阵乘法。

我有这个矩阵：

A <- matrix(c(1,-1,-1,0,-1,0,1,0,0,1,0,0,0,1,-1,1,-1,0,0,-1,1,0,1,0,1,-1,-1,1,-1,1), nrow = 6, ncol = 5)

A
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    0    0    1
[2,]   -1    0    1   -1   -1
[3,]   -1    0   -1    1   -1
[4,]    0    1    1    0    1
[5,]   -1    0   -1    1   -1
[6,]    0    0    0    0    1

我想得到两个不同的矩阵。第一个矩阵是这样的：

C
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    0    0    2    0    1
[2,]    0    0    2    1    2    0
[3,]    0    2    0    0    4    0
[4,]    2    1    0    0    0    1
[5,]    0    2    4    0    0    0
[6,]    1    0    0    1    0    0

这个“收敛矩阵”类似于 A 对其转置的乘法（在 R 中是这样的A%*%t(A)），但有一点扭曲，在求和期间获得每个单元格，我只想要正值的总和。例如，对于单元格 C23，常规总和将是：

(-1)(-1) + (0)(0) + (1)(-1) + (-1)(1) + (-1)(-1) = 0

，但我只想要正积的总和，在这个例子中，第一个 [(-1)(-1)] 和最后一个 [(-1)(-1)] 得到 2。

第二个矩阵是这样的：

D
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    2    2    0    2    0
[2,]    2    0    2    1    2    1
[3,]    2    2    0    2    0    1
[4,]    0    1    2    0    2    0
[5,]    2    2    0    2    0    1
[6,]    0    1    1    0    1    0

这个“散度矩阵”与前一个类似，不同之处在于我只想对负值的绝对值求和。例如，对于单元格 D23，常规总和为：

(-1)(-1) + (0)(0) + (1)(-1) + (-1)(1) + (-1)(-1) = 0

，但我只想要负积的绝对值之和，在这个例子中，第三个 abs[(1)(-1)] 和第四个 abs[(-1)(-1)] 得到 2。

我一直在尝试应用、扫描和循环，但我做不到。感谢您的回复。

4

3 回答 3

1

它的效率会大大降低，但您可以将矩阵分解为行向量列表，这样更容易计算。使用purrr，这对于列表很方便，

library(purrr)

A <- matrix(c(1,-1,-1,0,-1,0,1,0,0,1,0,0,0,1,-1,1,-1,0,0,-1,1,0,1,0,1,-1,-1,1,-1,1), 
            nrow = 6, ncol = 5)

C <- seq(nrow(A)) %>%    # generate a sequence of row indices
    map(~A[.x, ]) %>%   # subset matrix into a list of rows
    cross2(., .) %>%    # do a Cartesian join to get pairs of rows
    # calculate products, then subset before summing. Simplify to vector
    map_dbl(~{ij <- .x[[1]] * .x[[2]]; sum(ij[ij >= 0])}) %>% 
    matrix(nrow(A))    # reassemble to matrix

C
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    3    0    0    2    0    1
#> [2,]    0    4    2    1    2    0
#> [3,]    0    2    4    0    4    0
#> [4,]    2    1    0    3    0    1
#> [5,]    0    2    4    0    4    0
#> [6,]    1    0    0    1    0    1

# same except subsetting and `-` to make negatives positive
D <- seq(nrow(A)) %>% 
    map(~A[.x, ]) %>%
    cross2(., .) %>% 
    map_dbl(~{ij <- .x[[1]] * .x[[2]]; sum(-ij[ij <= 0])}) %>% 
    matrix(nrow(A))

D
#>      [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,]    0    2    2    0    2    0
#> [2,]    2    0    2    1    2    1
#> [3,]    2    2    0    2    0    1
#> [4,]    0    1    2    0    2    0
#> [5,]    2    2    0    2    0    1
#> [6,]    0    1    1    0    1    0
于 2017-04-14T21:07:50.710 回答
1

这是在基础 R 中的一次尝试。所以基本上你遵循矩阵交叉乘积方法，但你尝试sum手动管理该步骤：

f <- function(A, convergence=TRUE){
    sapply(seq_len(nrow(A)), function(i) {
        r <- t(matrix(A[i,],ncol(A),nrow(A)))*A
        if(convergence)
            r[r<0] <- 0
        else
            r[r>0] <- 0
        rowSums(abs(r))
    })
}

> f(A, convergence = TRUE)

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    3    0    0    2    0    1
[2,]    0    4    2    1    2    0
[3,]    0    2    4    0    4    0
[4,]    2    1    0    3    0    1
[5,]    0    2    4    0    4    0
[6,]    1    0    0    1    0    1

> f(A, convergence = FALSE)

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    2    2    0    2    0
[2,]    2    0    2    1    2    1
[3,]    2    2    0    2    0    1
[4,]    0    1    2    0    2    0
[5,]    2    2    0    2    0    1
[6,]    0    1    1    0    1    0
于 2017-04-14T21:08:13.557 回答
1

另一种做法：

D <- A
D[D<0] = -1i*D[D<0]
D <- Im(tcrossprod(D))

C <- tcrossprod(A) + D

A在问题中定义。输出：

> D
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    2    2    0    2    0
[2,]    2    0    2    1    2    1
[3,]    2    2    0    2    0    1
[4,]    0    1    2    0    2    0
[5,]    2    2    0    2    0    1
[6,]    0    1    1    0    1    0
> C
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    3    0    0    2    0    1
[2,]    0    4    2    1    2    0
[3,]    0    2    4    0    4    0
[4,]    2    1    0    3    0    1
[5,]    0    2    4    0    4    0
[6,]    1    0    0    1    0    1
于 2017-04-14T21:46:47.627 回答