How to implement this even-odd matrix decomposition efficiently?

Note: This question has also been asked on stackoverflow - see https://stackoverflow.com/questions/57197910/how-to-implement-this-even-odd-matrix-decomposition-efficiently?noredirect=1#comment100975572_57197910

In connection with solving a partial differential equation numerically, I have come across the following question about how to implement a specific matrix-vector product algorithm efficiently in Matlab. For notational convenience I will formulate the question using Matlab notation.

The background of the question is that I have a matrix, D, which has size NxN (N is even and typical values for N are 4, 6 and 8) and which has the property that

D(i,j) = - D(N+1-i, N+1-j)


for 1 <= i,j <= N. To solve the differential equation, I use a numerical method that requires D to be multiplied onto many different vectors and I therefore want to exploit the above property to make the matrix-vector product as efficient as possible.

The algorithm which utilizes the symmetry of D that I have in mind is the following: Assume we want to compute the product D*f where f is a Nx1 column vector. Moreover, define the column vectors e and o of size N/2x1 by

e(n) = 0.5*(f(n) + f(N+1-n))

o(n) = 0.5*(f(n) - f(N+1-n))


as well as the matrices De and Do of size N/2xN/2 by

De(n,m) = D(n,m) + D(n,N+1-m)

Do(n,m) = D(n,m) - D(n,N+1-m).


Using the definition of a matrix-vector product it can be shown that the matrix vector product D * f can be computed as

Df = [De*e + Do*o; -De*e + Do*o]


and since the products De * e and Do * o each require one quarter as many multiplications as the product D * f I hope that it is possible to reduce the computational time by roughly a factor of 2 using this algorithm.

So far I have, however, not been able to implement the algorithm in a way such that it is faster than the built in matrix-vector product which does not utilize the special property of D. The first hirdle seems to be constructing the vectors e and o. If I do it the straight forward way, i.e.

e = 0.5*(f(1:N/2) + flipud(f(N/2+1:N)))

o = 0.5*(f(1:N/2) + flipud(f(N/2+1:N)))


it already takes more than an order of magnitude longer than to compute D * f directly. Forming the sparse matrices eMat and oMat as

eMat = 0.5*[1, 0, ..., 0, 0, ..., 0, 1;

0, 1, ..., 0, 0, ..., 1, 0;

|  |   |   |  |   |   |  |

0, 0, ..., 1, 1, ..., 0, 0]

oMat = 0.5*[1, 0, ..., 0,  0, ...,  0, -1;

0, 1, ..., 0,  0, ..., -1,  0;

|  |    |   |  |    |   |   |

0, 0, ..., 1, -1, ...,  0, 0]


where the vertical lines mean that the pattern should be continued, e and o can be computed as

e = eMat*f

o = oMat*f


and this takes about a third of the time it takes computing D * f directly. Using the sparse matrices eMat and oMat I can compute De * e and Do * o as

Dee = De*eMat*f

Doo = Do*oMat*f


and this also takes roughly a third of the time it takes computing D * f with the built in matrix-vector product. Now the difficulty arises, because carrying out the computation

Dee + Doo


requires about the same amount of time (and typically even a little more) than computing D * f directly. Needless to say, it also takes more time to compute D * f as

Df = [Dee + Doo; -Dee + Doo]


than with the built in matrix-vector product.

My question is if anyone can suggest a way of implementing the above algorithm efficiently in Matlab, such that the algorithm in total runs faster than the built in matrix-vector product?

• For small $N$, ($N \leq 8$) I do not think that a pure MATLAB implementation can be much faster than D*f; on the contrary a C implementation could give the expected speedup. Is writing a mex file an option? Commented Aug 2, 2019 at 23:12

statement of problem

The OP question is how to compute efficiently the dense matrix-vector product $$Df$$ where $$D\in\mathbb{R}^{N\times N}$$, with the property

$$D = -J\,D\,J$$

and $$J$$ is the $$N\times N$$ exchange matrix. Assumptions are that matrix $$D$$ is constant and $$N$$ even and small.

While reading the OP it occurred to me that the rows of Emat and Omat are the eigenvectors of the exchange matrix associated to the eigenvalues $$\lambda = 1$$ and $$\lambda = -1$$, respectively, and that after renormalization the eigenvectors are orthonormal. Therefore the OP algorithm can be rewritten as

$$V \, (\frac14 V^T\,D\,V) \, V^T f$$

with

$$\newcommand\iddots{\mathinner{ \kern1mu\raise1pt{.} \kern2mu\raise4pt{.} \kern2mu\raise7pt{\Rule{0pt}{7pt}{0pt}.} \kern1mu }}$$

$$V = \begin{bmatrix} 1 & 0 & \cdots & 0 & 1 & 0 & \cdots & 0 \\ 0 & 1 & \cdots & 0 & 0 & 1 & \cdots & 0 \\ 0 & 0 & \ddots & & & & \ddots & 0 \\ 0 & 0 & \cdots & 1 & 0 & 0 & \cdots & 1 \\ 0 & 0 & \cdots & 1 & 0 & 0 & \cdots & -1 \\ 0 & 0 & \iddots & & & & \iddots & 0 \\ 0 & 1 & \cdots & 0 & 0 & -1 & \cdots & 0 \\ 1 & 0 & \cdots & 0 & -1 & 0 & \cdots & 0 \\ \end{bmatrix}$$

and noting that $$V^T V = V V^T = 2 I$$

It is easy to show that $$V^T\,D\,V$$ is block northeast-to-southwest diagonal, i.e.

$$V^T\,D\,V = \begin{bmatrix} 0 & D_1 \\ D_2 & 0 \end{bmatrix}$$

where $$D_i$$ have dimensions $$N/2 \times N/2$$.

It is therefore confirmed that it is possible to rewrite $$Df$$ so that a $$N\times N$$ matrix-vector product is substituted by two $$N/2 \times N/2$$ products, with a theoretical speedup of 2.

MATLAB implementation

The idea here is to use a function closure to construct a function handle mfh that computes D*f. In MATLAB this is obtained with nested functions: the outer scope computes the constant matrices $$V, D_1, D_2$$, while the inner function mf computes the $$V^T f$$ decomposition, the two dense matrix-vector products, and the $$V\cdot$$ composition.

function mfh = makemf(D)

[N,M] = size(D);
assert(N==M);

HN = round(N/2);
HN1 = HN + 1;

V = sparse(N,N);
for i = 1:HN
V(i,i) = 1;
V(N+1-i,i) = 1;
V(i,HN+i) = 1;
V(N+1-i,HN+i) = -1;
end

D1 = V(:,1:HN)'*D*V(:,HN1:end) / 4;
D2 = V(:,HN1:end)'*D*V(:,1:HN) / 4;

function y = mf(x)
x = V'*x;
y = zeros(size(x));
y(1:HN,:) = D1*x(HN1:end);
y(HN1:end,:) = D2*x(1:HN);
y = V*y;
end

mfh = @mf;

end


This function returns a function handle mf so that mf(f) == D*f. E.g.

>> N = 8;
>> D = randn(N);
>> D = flipud(D) - fliplr(D);
>> f = randn(N,1);
>> mf = makemf(D);
>> norm(mf(f) - D*f) / norm(D*f)

ans =

1.9575e-16


timings

Now some timings of mf(f) (labelled 'odd even decom.') against D*f (labelled 'direct') for N from 8 to 8000.

and the speedup

Timings where obtained with timeit

t = @() D*f;
timeit(t);
t = @() mf(f);
timeit(t);


The bottom line is that this implementation attains the expected asymptotics for $$N > 1000$$ and that for $$N=8000$$ the measured speedup is very close to the theoretical one.

closing comments

This is a late answer that I wrote mainly for fun. In fact it is expected that the theoretical asymptotics can be reached only for medium sized matrices, and that in pure MATLAB for very small $$N$$ there is no way to speed up D*f.

Nevertheless I hope that some useful programming techniques are demonstrated in this answer:

• function closures,
• preallocation of output arrays (avoid constructs like [D_1*y_1; D_2*y_2] but first create an empty $$N$$ array and then assign to its $$N/2$$ portions),
• sparse matrices.

I think that for $$N=10$$ the only way to obtain a significative speedup would be a C function inside a MEX file.

• That is a very impressive answer that uses quite a few tricks I was not aware of! Thank you for your time and your help! Commented Aug 6, 2019 at 6:31
• Indeed, nice answer!
– Bob
Commented Nov 17, 2022 at 20:43