vectoriced norm/matrix multiplication

Question

I have a norm that is described by the matrix sigma

sigma <- matrix(c(1,0.5,0,0.5,1,0,0,0,1),3,3))

To compute the norm of a vector I compute

t(x) %*% sigma %*% x

which works fine for a vector e.g. x = 1:3.

However i want to compute the norm of many vectors at the same time, that is I have

x <- t(matrix(rep(1:3, 10),3,10))

(of course filled with different entries).

Is there a way to compute the norm of each vector simultanuously? I.e. something like

lapply(1:10, function(i) t(x[i,]) %*% sigma %*% x[i,])

jogo · Accepted Answer · 2017-10-20T09:53:51.727

You can do:

sigma <- matrix(c(1,0.5,0,0.5,1,0,0,0,1),3,3)
x <- t(matrix(rep(1:3, 10),3,10))

mynorm <- function(x, sig) t(x) %*% sig %*% x
apply(x, 1, mynorm, sig=sigma)

Here is a variant with tcrossprod():

mynorm <- function(x, sig) tcrossprod(x, sig) %*% x
apply(x, 1, mynorm, sig=sigma)

And here is the benchmark (including variants of the solution from compute only diagonals of matrix multiplication in R , thanks to @Benjamin for the link):

mynorm1 <- function(x, sig) t(x) %*% sig %*% x
mynorm2 <- function(x, sig) tcrossprod(x, sig) %*% x

microbenchmark(n1=apply(x, 1, mynorm1, sig=sigma), 
               n2=apply(x, 1, mynorm2, sig=sigma), 
               n3 = colSums(t(x) * (sigma %*% t(x))),
               n4 = rowSums(x * t(sigma %*% t(x))),
               n5 = rowSums(x * (x %*% t(sigma) )),
               n6 = rowSums(x * tcrossprod(x, sigma)),
               Eugen1 = diag(x %*% sigma %*% t(x)),
               Eugen2 = diag(x %*% tcrossprod(sigma, x)),
               unit="relative")

EugenR · Answer 2 · 2017-10-20T09:32:39.097

2

What do you think about this simple matrix multiplication:

 diag(t(x) %*% sigma %*% x)

Edit: after the matrix multiplications you need the diagonal (of course).

Then it is faster then solutions with apply

edited Oct 20 '17 at 09:32

answered Oct 20 '17 at 08:49

EugenR

176
6

This gives me "Error in t(x) %*% sigma : non-conformable arguments" – Jürg W. Spaak Oct 20 '17 at 08:50
@EugenR You are calculating unneeded elements. – jogo Oct 20 '17 at 09:34
@EugenR Yes, could be even faster (see variant Eugen2 in the benchmark). The speed eventually depends from the dimensions of the matrices. – jogo Oct 20 '17 at 09:55

score 2 · Answer 3 · answered Oct 20 '17 at 09:10

2

This should do

> sigma <- matrix(c(1,0.5,0,0.5,1,0,0,0,1),3,3)
> x <- t(matrix(rep(1:30, 10),3,10))
> 
> # should give
> t(x[1, ]) %*% sigma %*% x[1, ]
     [,1]
[1,]   16
> t(x[2, ]) %*% sigma %*% x[2, ]
     [,1]
[1,]   97
> 
> # which you can get by
> rowSums((x %*% sigma) * x)
 [1]   16   97  250  475  772 1141 1582 2095 2680 3337

answered Oct 20 '17 at 09:10

Benjamin Christoffersen

4,703
15
37

1

The question is close to being a duplicate of https://stackoverflow.com/q/21708489/5861244 – Benjamin Christoffersen Oct 20 '17 at 09:11

vectoriced norm/matrix multiplication

3 Answers3