Sum of products which indice is not i an j

Question

I am implementing the following function in the R code:

So far I used:

sig.TOM <- function(adj, sig.adj) {
out <- matrix(nrow = nrow(adj), ncol = ncol(adj))
  for (i in 1:nrow(adj)) {
    for (j in 1:ncol(adj)) {
      out[i,j] <- abs(adj[i, j] + sum(sig.adj[i, -c(i, j)]*sig.adj[-c(i, j), i]))/(
        min(sum(sig.adj[-i, i]), sum(sig.adj[-j, j])) + 1 - abs(adj[i,j]))
    }
  }
  return(out)
}

where ~a is the following mock matrix:

sig.adj <- structure(c(1, -0.418913311940584, 1, 0.947013383275973, -1, 
-0.418913311940584, 1, -0.207962861914701, 0.584386281408348, 
-0.687223049826016, 1, -0.207962861914701, 1, 0.763551721347657, 
-0.0327147711077901, 0.947013383275973, 0.584386281408348, 0.763551721347657, 
1, 0.284466543760789, -1, -0.687223049826016, -0.0327147711077901, 
0.284466543760789, 1), .Dim = c(5L, 5L))

and adj <- abs(sig.adj), where adj in the formula is described as a and sig.adj as ~a.

But the result is not symmetric as expected so I must have implemented it wrong, I have doubts on the summation part.

How can that sum of products of values when the indices are not i or j be implemented?

The proposed solutions:

spec.mult1 <- function(A,B) {
  rA <- nrow(A); cB <- ncol(B)
  C <- A %*% B
  for (i in 1:rA) for (j in 1:cB) 
    C[i,j] <- C[i,j] - A[i,i]*B[i,j] - A[i,j]*B[j,j] + ifelse(i==j, A[i,i]*B[j,j], 0) 
  C
}

spec.mult2 <- function(A) {
  dA.A <- diag(A)*A
  crossprod(A) - dA.A - t(dA.A) + diag(diag(A)^2)
}

spec.mult3 <- function(A,B) {
  rA <- nrow(A); cB <- ncol(B)
  C <- A %*% B
  for (i in 1:rA) for (j in 1:cB) 
    C[i,j] <- C[i,j] - A[i,i]*B[i,j] - A[i,j]*B[j,j] 
  C
}

spec.mult4 <- function(A) {
  dA.A <- diag(A)*A
  crossprod(A) - dA.A - t(dA.A)
}

spec.mult5 <- function(sig.adj) {
  nr <- nrow(sig.adj); nc <- ncol(sig.adj)
  C <- matrix(NA, nr, nc)
  for (i in 1:nr) for (j in 1:nc) 
    C[i,j] <- sum(sig.adj[i, -c(i, j)]*sig.adj[-c(i, j), j])
  C
}

Comparing the results of each function :

all(res1 == res2)
[1] TRUE
> all(res1 == res3)
[1] FALSE
> all(res1 == res4)
[1] FALSE
> all(res1 == res5)
[1] FALSE
> all(res2 == res3)
[1] FALSE
> all(res2 == res4)
[1] FALSE
> all(res2 == res5)
[1] FALSE
> all(res3 == res4)
[1] TRUE
> all(res3 == res5)
[1] FALSE
> all(res4 == res5)
[1] FALSE

Results that, spec.mult1 == spec.mult2 and spec.mult3 == spec.mult4 but spec.mult5 (the one I understand, and hope it is correct) doesn't appear

Why are you calculating `... - abs(adj[i,j])` when `adj <- abs(sig.mat)`? — jogo, Sep 01 '16 at 11:03
Thanks @jogo, one answer tells me what is wrong in my formula, and the other calculates the formula without explaining what is made. (My formula corrected per Vandenman and your formula differ) And neither of them explains how to write the "sum of products of values when the indices are not i or j". Thus I am not accepting neither of them even if they were helpful. — llrs, Sep 14 '16 at 06:47
Thanks for your proposed solutions. I think now I learned how to do it, thanks for your patience. — llrs, Sep 15 '16 at 14:01
Thanks for comparing the results of the different functions. The difference `spec.mult5(sig.adj) - spec.mult4(sig.adj)`is surprizing for me. all other is clear. BTW: is is not a good idea to compare floating point numbers per `==` , see http://stackoverflow.com/questions/9508518/why-are-these-numbers-not-equal — jogo, Sep 15 '16 at 14:13
I was wrong in thinking about omitted indicies (how they work in R). I edited my answer. Now all functions give the same result. — jogo, Sep 15 '16 at 18:11

score 2 · Answer 1 · answered Sep 01 '16 at 11:03

2

I think you indexed incorrectly in the sum over u!=i, j. The part

sum(sig.adj[i, -c(i, j)]*sig.adj[-c(i, j), i])

should be

sum(sig.adj[i, -c(i, j)]*sig.adj[-c(i, j), j])

With your example the output is then a symmetric matrix for me.

answered Sep 01 '16 at 11:03

Vandenman

3,046
20
33

jogo · Accepted Answer · 2016-09-15T18:21:18.307

$C_{ij} = \sum_{u} a_{iu} b_{uj})$ is a normal matrix multiplication. So you can get $\sum_{u \ne i,j} a_{iu} b_{uj}$ by correcting the result of a matrix multiplication (i.e. subtracting the unwanted parts of the sum).
Take care of the fact that in the case of i==j only one part of $\sum_{u} a_{iu} b_{uj})$ has to be neglected.

something about omitted indicies:

A <- matrix(c(1:4, 2,5:7, 3,6,8:9, 4,7,9,10), 4,4)
A[1, -c(1,1)]

The element is omitted only once.

The part of the formula beside the special matrix multiplication is clear:

sig.TOM <- function(sig.adj) {
  adj <- abs(sig.adj)
  k <- colSums(adj) - abs(diag(adj))
  abs(adj + spec.mult(sig.adj, sig.adj)) / (outer(k,k, pmin) +1 - abs(adj))
}
sig.TOM(sig.adj)

(or spec.mult(sig.adj) for a one-argument variant of the special matrix multiplication, see below) p.s.: I copied the part ... +1 - abs(adj) from your question because I don't know if you want ... +1 - adj or ... +1 - sig.adj

Here are five variants of the special matrix multiplication:

A <- matrix(c(1:4, 2,5:7, 3,6,8:9, 4,7,9,10), 4,4)
A[1, -c(1,1)]

spec.mult1 <- function(A,B) {
  rA <- nrow(A); cB <- ncol(B)
  C <- A %*% B
  for (i in 1:rA) for (j in 1:cB) 
    C[i,j] <- C[i,j] - A[i,i]*B[i,j] - A[i,j]*B[j,j] + ifelse(i==j, A[i,i]*B[j,j], 0) 
  C
}

spec.mult2 <- function(A) {
  dA.A <- diag(A)*A
  crossprod(A) - dA.A - t(dA.A) + diag(diag(A)^2)
}

spec.mult3 <- function(A,B) {
  rA <- nrow(A); cB <- ncol(B)
  C <- A %*% B
  for (i in 1:rA) for (j in 1:cB) 
    C[i,j] <- C[i,j] - A[i,i]*B[i,j] - A[i,j]*B[j,j] 
  for (i in 1:rA) C[i,i] <- C[i,i] + A[i,i]*B[i,i]
  C
}

spec.mult4 <- function(A) {
  dA <- diag(A)
  dA.A <- dA*A
  crossprod(A) - dA.A - t(dA.A) + diag(dA^2)
}

spec.mult5 <- function(sig.adj) {
  nr <- nrow(sig.adj); nc <- ncol(sig.adj)
  C <- matrix(NA, nr, nc)
  for (i in 1:nr) for (j in 1:nc) 
    C[i,j] <- sum(sig.adj[i, -c(i, j)]*sig.adj[-c(i, j), j])
  C
}

spec.mult1(A, A) - spec.mult5(A) 
spec.mult2(A) - spec.mult5(A)
spec.mult3(A, A) - spec.mult5(A) 
spec.mult4(A) - spec.mult5(A)

Sorry for bothering you again, but I don't understand the spec.mult function, how are they equivalent? How could I calculate the sum of the absolute product of a_iu a_uj with it? — llrs, Sep 01 '16 at 14:18
yes, they are equivalent, but I don't understand why. There is another formula with a slightly different approach, see the equation 4 of the [report](https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/TechnicalReports/signedTOM.pdf) I am trying to use, which I would like to implement too — llrs, Sep 01 '16 at 14:39
I would prefer to understand how to do the $\sum_{u \ne i,j} a_{iu} a_{uj}$ or $\sum_{u \ne i,j} \abs(a_{iu} a_{uj})$. And Whatever it is spec.mult I don't understand it, or how to tweak it. — llrs, Sep 14 '16 at 12:38

Sum of products which indice is not i an j

2 Answers2