1

我有一个大而密集的矩阵 A，我的目标是使用迭代方法找到线性系统 Ax=b 的解决方案（在 MATLAB 中是使用其内置 GMRES 的计划）。对于超过 10,000 行，这对于我的计算机来说存储在内存中太多了，但我知道 A 中的条目是由长度为 N 的两个已知向量 x 和 y 构造的，并且条目满足： A(i,j) = .5*(x[i]-x[j])^2+([y[i]-y[j])^2 * log(x[i]-x[j])^2+([y [i]-y[j]^2)。

MATLAB 的 GMRES 命令接受一个函数调用作为输入，该函数调用可以计算矩阵向量乘积 A*x，这使我能够处理比我可以存储在内存中更大的矩阵。为了编写矩阵向量乘积函数，我首先在 matlab 中逐行尝试并使用了一些向量化，但我避免生成整个数组 A（因为它会太大）。不幸的是，这在我申请 GMRES 时相当缓慢。我的计划是用 C 语言为 MATLAB 编写一个 mex 文件，理想情况下应该比 matlab 代码快得多。我对 C 比较陌生，所以这很糟糕，而且我在 C 中编写代码的天真尝试比我在 Matlab 中的部分向量化尝试要慢。

#include <math.h>
#include "mex.h"
void Aproduct(double *x, double *ctrs_x, double *ctrs_y, double *b, mwSize n)
{
    mwSize i;
    mwSize j;
    double val;
    for (i=0; i<n; i++) {
        for (j=0; j<i; j++) {
            val = pow(ctrs_x[i]-ctrs_x[j],2)+pow(ctrs_y[i]-ctrs_y[j],2);

            b[i] = b[i] + .5* val * log(val) * x[j];
        }
        for (j=i+1; j<n; j++) {
            val = pow(ctrs_x[i]-ctrs_x[j],2)+pow(ctrs_y[i]-ctrs_y[j],2);

            b[i] = b[i] + .5* val * log(val) * x[j];
        }
    }
}

以上是matlab mex文件代码的计算部分（如果我理解正确的话，稍微修改了C）。请注意，我跳过了 i=j 的情况，因为在这种情况下，变量 val 将是 0*log(0)，对我来说应该解释为 0，所以我只是跳过它。

有没有更有效或更快的方法来写这个？当我通过 matlab 中的 mex 文件调用这个 C 函数时，它非常慢，甚至比我使用的 matlab 方法还要慢。这让我感到惊讶，因为我怀疑 C 代码应该比 matlab 快得多。

我与之比较的部分矢量化的替代 matlab 方法是

function Ax = Aprod(x,ctrs)
n = length(x);
Ax = zeros(n,1);
for j=1:(n-3)
    v = .5*((ctrs(j,1)-ctrs(:,1)).^2+(ctrs(j,2)-ctrs(:,2)).^2).*log((ctrs(j,1)-ctrs(:,1)).^2+(ctrs(j,2)-ctrs(:,2)).^2);

    v(j)=0;
    Ax(j) = dot(v,x(1:n-3);
end

（n-3是因为实际上有3个额外的组件，但它们是分开处理的，所以我排除了该代码）。这是部分矢量化的，只需要一个 for 循环，所以它更快是有道理的。但是，我希望使用 C+mex 文件可以更快。

任何建议或帮助将不胜感激！谢谢！

编辑：我应该更清楚。我对任何可以帮助我使用 GMRES 反转我感兴趣的矩阵的更快方法持开放态度，这需要一种更快的方法来执行矩阵向量乘积，而无需将数组显式加载到内存中。谢谢！

4

1 回答 1

1

如果您有Parallel Computing Toolbox和MATLAB Distributed Computing Server，您可以直接使用反斜杠求解大型密集线性系统。（如果您没有可用的集群，您可能希望使用Amazon EC2 机器）。像这样： http: //www.mathworks.co.uk/help/distcomp/examples/benchmarking-ab.html

于 2013-11-01T08:27:28.277 回答