fft矩阵向量乘法

Question

我必须在MATLAB中求一个线性方程组A*x=B,其中A对称且其元素取决于指数的差异:Aij=f(i-j).

我使用迭代求解器,因为它的大小A是40000x40000.迭代求解器需要确定产品A*x在哪里x是测试解决方案.对该产品的评估结果证明是卷积,因此可以通过快速傅里叶变换(cputime~ Nlog(N)而不是N^2)来完成.我对这个问题有以下问题:

这是卷积通告吗？因为如果它是循环的,我认为我必须使用新矩阵的特定索引来获取fft.是对的吗？
我发现难以为fft编写例程,因为我无法理解我应该使用的索引.是否有任何现成的例程,我可以用fft直接评估产品A*x而不是卷积？实际上,矩阵A由3×3块构成并且是对称的.产品的现成例程A*x对我来说是最好的解决方案.
如果没有现成例程,您能否通过示例给我一个想法如何构建此例程以通过fft评估矩阵向量乘积？

先感谢您,

帕诺斯

Answer 1

非常好的且有趣的问题！:) 对于某些特殊的矩阵结构，Ax = b问题可以很快解决。

循环矩阵。

循环卷积对应的矩阵Ax = h*x（* - 为卷积符号）在傅立叶域中进行对角化，可以通过以下方式求解：

x = ifft(fft(b)./fft(h));

呈三角形和带状。

通过稀疏 LU 分解可以有效求解三角矩阵和对角主导带状矩阵：

[L,U] = lu(sparse(A)); x = U\(L\b);

泊松问题。

如果A是拉普拉斯算子的有限差分近似，则可以通过多重网格方法有效地解决该问题（例如，网络搜索“matlab multigrid”）。

Answer 2

有趣的问题！

在您的情况下，卷积不是圆形的，除非您施加额外的条件。例如，A(1,3)应该等于A(2,1)等。

您可以使用conv（仅保留带有选项的非零填充部分valid）来做到这一点，这可能也是 N*log(N)。例如，让

A = [a b c d
     e a b c
     f e a b
     g f e a];

然后A*x是一样的

conv(fliplr([g f e a b c d]),x,'valid').'

或者更一般地说，A*x与

conv(fliplr([A(end,1:end-1) A(1,:)]),x,'valid').'