How to implement this even-odd matrix decomposition efficiently?

Question

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?

Answer 1

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.