# Fast vector - "diagonal" matrix multiplication

Let $\mathbf{1}\in\mathbb{R}^d$ be a vector with all elements equal to $1$. Define: $$\mathbf{D} = \mathrm{diag}(\mathbf{1}^\top,\mathbf{1}^\top,\ldots,\mathbf{1}^\top) = \begin{bmatrix} 1 \cdots 1 & & \\ & 1\cdots 1 & \\ & & 1\cdots 1 \end{bmatrix} \in\mathbb{R}^{d\times d^2}$$

I would like to compute $\mathbf{D}\cdot \mathbf{x}$ for any vector $\mathbf{x}\in\mathbb{R}^{d^2}$. Could anybody suggest me a way to efficiently compute this product? (Since $\mathbf{D}$ has a very special structure, I guess there should be such a way).

This can be interpreted as summing over an index of a tensor when the vector $x$ is reshaped into a box of numbers instead of a list. In particular, if $X$ is the $d\text{-by-}d$ folded version of $x$, then the operation you are doing is, \begin{align} Dx &= \mathrm{vec}\left((I \otimes \mathbf{1})\mathrm{vec}(X)\right) \\ &= \mathrm{vec}(\mathbf{1}^T X I) \\ &= \mathrm{vec}(\mathbf{1}^T X). \end{align}

The matlab code to do this is surprisingly simple. It is,

sum(reshape(x,d,d))'


Here's an example of it in action - you can see that it far outperforms the standard dense multiply, sparse matrix multiply, and for loop versions:

>> onesmatrixquestion
dense matrix multiply
Elapsed time is 0.000873 seconds.
sparse matrix multiply
Elapsed time is 0.000115 seconds.
for loop version
Elapsed time is 0.000154 seconds.
tensorized version
Elapsed time is 0.000018 seconds.


Here's the code that generated those timing results:

%onesmatrixquestion.m
d = 100;
onevec = ones(1,d);
D = kron(eye(d,d),onevec);
Dsparse = sparse(D);

x = randn(d^2,1);

disp('dense matrix multiply')
tic
aa = D*x;
toc

disp('sparse matrix multiply')
tic
bb = Dsparse*x;
toc

disp('for loop version')
tic
cc = zeros(d,1);
ind = 1;
for kk=1:d
for jj=1:d
cc(kk) = cc(kk) + x(ind);
ind = ind + 1;
end
end
toc

disp('tensorized version')
tic
dd = sum(reshape(x,d,d))';
toc

if (norm(aa - bb) > 1e-9 || norm(bb - cc) > 1e-9 ...
|| norm(cc - dd) > 1e-9)
disp('error: different methods give different results')
end

• Hi @Nick. I've realized that it can be done even better: ee = (onevec*reshape(x,d,d))';
– f10w
Commented May 28, 2014 at 13:13
• Does it go faster that way as compared to using the sum command? Commented May 28, 2014 at 14:43
• Yes, it does. That's why I posted that comment ;)
– f10w
Commented May 28, 2014 at 15:08
• Oh, interesting. I was thinking that the 'sum' command would be specially optimized and therefore faster than a generic matvec - I guess not! Commented May 28, 2014 at 15:31
• Hmm you're right. Right multiplying a matrix by a vector of ones seems to be about twice as fast as using the sum command in Matlab. Commented May 28, 2014 at 16:11

This is equivalent to computing sums of consecutive contiguous subvectors of $\mathbf{x}$. You won't do much better than simple hand-coded nested loops if you have an automatically vectorizing compiler since you will probably be memory bandwidth limited for large $d$.

• Hi Victor. I've got an answer. Thanks anyway :D
– f10w
Commented May 28, 2014 at 9:48

Edit: removed confusing code

Yes there is a shortcut, here's some Python code:



import numpy as np

d = 3

D = np.repeat(np.identity(d), d, axis=1)
print('D:')
print D
print

print('x:')
x = np.arange(d*d, dtype=float)
print x
print

print('D * x:')
print D.dot(x)
print

print('shortcut:')
print x.reshape((d, d)).sum(axis=1)
print



output:


D:
[[ 1.  1.  1.  0.  0.  0.  0.  0.  0.]
[ 0.  0.  0.  1.  1.  1.  0.  0.  0.]
[ 0.  0.  0.  0.  0.  0.  1.  1.  1.]]

x:
[ 0.  1.  2.  3.  4.  5.  6.  7.  8.]

D * x:
[  3.  12.  21.]

shortcut:
[  3.  12.  21.]

• How is that a shortcut? It doesn't take advantage of sparsity at all and just directly computes the dot product (which is the same as matrix vector multiplication here). Commented May 27, 2014 at 18:20
• @doug Sorry my answer was unnecessarily long and confusing. First I compute the product directly, but if you keep reading you will find that I then compute the product in a second way using a trick.
– k20
Commented May 27, 2014 at 18:40
• Oh, I see now. I suggest adding at least a minimal explanation or comments in your code. Commented May 27, 2014 at 19:06