1
$\begingroup$

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?

$\endgroup$
1
  • 1
    $\begingroup$ 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? $\endgroup$
    – Stefano M
    Commented Aug 2, 2019 at 23:12

1 Answer 1

3
$\begingroup$

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.

run times

and the speedup

enter image description here

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.

$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$ Commented Aug 6, 2019 at 6:31
  • $\begingroup$ Indeed, nice answer! $\endgroup$
    – Bob
    Commented Nov 17, 2022 at 20:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.