Multiply by diagonal matrix

Multiply by diagonal matrix using efficient algorithm

Usage

dmult(M, v, side = c("left", "right"))

Arguments

M

matrix

v

vector with entries forming a diagonal matrix matching the dimensions of M depending on the value of side

side

is the matrix M "left" or "right" multiplied by the diagonal matrix

Value

matrix product

Details

Naively multiplying by a diagonal matrix with p entries takes \(\mathcal{O}(p^2)\), even though must values in the diagonal matrix are zero. R has built in syntax to scale columns, so diag(v) %*% M can be evaluated with v * M. This is often difficult to remember, hard to look up, and scaling rows requires t(t(M) * v). This is hard to read and write, and requires two transpose operations.

Here, dmult() uses Rcpp to evaluate the right multiplication efficiently, and uses v * M for left multiplication. This gives good performance and readability.

In principle, the Rcpp code can be modified to use BLAS for better performance by treating a NumericMatrix as a C array. But this is not currently a bottleneck

Examples

# right multiply
# mat %*% diag(v)
n <- 30
p <- 20
mat <- matrix(n * p, n, p)
v <- rnorm(p)

A <- dmult(mat, v, side = "right")
B <- mat %*% diag(v)
range(A - B)
#> [1] 0 0

# left multiply
# diag(v) %*% mat
n <- 30
p <- 20
mat <- matrix(n * p, p, n)
v <- rnorm(p)

A <- dmult(mat, v, side = "left")
B <- diag(v) %*% mat
range(A - B)
#> [1] 0 0