Matrix Multiplication with Diagonal Matrices

## Keywords

Matrix algebra, Matlab, R

Unreviewed

# Description

There are many instances where one encounters a matrix multiplication of the type $DMD$, where $D$ is a diagonal matrix and $M$ is some full (dense) matrix.

(1)
\begin{align} \nonumber DMD = \underbrace{\left[ \begin{array}{cccc} d_{1} & 0 & \ldots & 0\\ 0 & d_{2} & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & d_{n}\\ \end{array} \right] }_{D} \underbrace{ \left[ \begin{array}{cccc} m_{11} & m_{12} & \ldots & m_{1n}\\ m_{21} & m_{22} & \ldots & m_{2n}\\ \vdots & \vdots & \ddots & \vdots\\ m_{n1} & m_{n2} & \ldots & m_{nn}\\ \end{array} \right] }_{M} \underbrace{ \left[ \begin{array}{cccc} d_{1} & 0 & \ldots & 0\\ 0 & d_{2} & \ldots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \ldots & d_{n}\\ \end{array} \right]}_{D} \end{align}

One example arises when working with variance-covariance matrices: these matrices can be decomposed such that $D$ holds the random variables' standard deviations as its diagonal elements, and $M$ is the correlation matrix.

Numerically, a matrix multiplication $DMD$ can often be sped up considerably by computing $\operatorname{diag}(D) \operatorname{diag}(D)'$ and then multiplying this rank-one matrix elementwise with $M$. (Here, $\operatorname{diag}$ is an operator that creates a column vector out of a matrix's main diagonal.) This computation requires $3n^2$ operations, while the operation count for the full matrix product is $4n^3$. Of course, for very small matrices, the speed-up may be inconsequential. Note also that if $M$ has a particular structure (e.g., is symmetric as in the case of a correlation matrix), even more efficient algorithms may be applied.

## R

# set size of matrix
N <- 500

# create matrices
D <- diag(runif(N))
M <- array(runif(N * N), dim = c(N, N))

# compute product / compare time
system.time(Z1 <- D %*% M %*% D)
system.time(Z2 <- outer(diag(D), diag(D)) * M)

# check difference between matrices
max(abs(Z1 - Z2))


## Matlab

% set size of matrix
N = 500;

% create matrices
D = diag(rand(N, 1));
M = rand(N, N);

% compute product / compare time
tic, Z1 = D * M * D; toc
tic, Z2 = diag(D) * diag(D)' .* M; toc

% check difference between matrices
max(max(abs(Z1 - Z2)))