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)))


# Internal Links

 Concepts Tutorials Tips Related Articles

# External links

 References 1. Gilli, M., D. Maringer and E. Schumann. (2011). Numerical Methods and Optimization in Finance. Elsevier. 2. R Development Core Team (2008). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. http://www.R-project.org. Weblinks
page revision: 42, last edited: 29 Sep 2011 11:26
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License