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?
D*f
; on the contrary a C implementation could give the expected speedup. Is writing a mex file an option? $\endgroup$