如何有效地实现这种奇偶矩阵分解？

Question

关于数值求解偏微分方程的问题，我遇到了以下问题：如何在Matlab中有效地实现特定的矩阵矢量积算法。为了符号上的方便，我将使用Matlab符号来表达问题。

问题的背景是我有一个矩阵D，其大小为NxN（N为偶数，N的典型值为4、6和8），其性质为

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

对于1 <= i，j <=N。为了解决微分方程，我使用了一种数值方法，该方法要求将D乘以许多不同的向量，因此我想利用上述性质来使矩阵向量乘积有效尽可能。

我想到的利用D对称性的算法如下：假设我们要计算乘积D * f，其中f是Nx1列向量。此外，定义大小为N / 2x1的列向量e和o

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

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

以及大小为N / 2xN / 2的矩阵De和Do

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

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

使用矩阵向量积的定义，可以证明矩阵向量积D * f可以计算为

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

并且由于乘积De * e和Do * o各自需要乘积乘以乘积D * f四分之一，所以我希望使用这种算法可以将计算时间减少大约2倍。

但是，到目前为止，我还无法以比不使用D的特殊属性的内置矩阵向量乘积更快的方式来实现该算法。第一个障碍似乎是在构造向量e和o。如果我直接做到这一点，即

 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)))

它比直接计算D * f所需的时间长了一个数量级。将稀疏矩阵eMat和oMat形成为

 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]

垂直线表示该图案应连续，其中e和o可计算为

 e = eMat*f

 o = oMat*f

这大约需要直接计算D * f的时间的三分之一。使用稀疏矩阵eMat和oMat，我可以计算De * e和Do * o为

 Dee = De*eMat*f

 Doo = Do*oMat*f

并且使用内置矩阵向量乘积计算D * f大约需要花费时间的三分之一。现在出现了困难，因为进行计算

 Dee + Doo

与直接计算D * f所需的时间大致相同（通常甚至更多）。不用说，将D * f计算为

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

而不是内置矩阵向量乘积。

我的问题是，是否有人可以提出一种在Matlab中有效实现上述算法的方法，从而使该算法的总体运行速度比内置矩阵矢量乘积快？