如何使用符号方法或数值方法加速Matlab积分计算？

Question

两个月前，我成为了一个Matlab的初学者。我用它来处理MRI图像的夏季项目。最近，我编写了下面展示的积分代码。但是，这两种方法都非常缓慢。运行它们需要一整天的时间。我该如何改进它们以缩短运行时间呢？

1. Symbolic method:

   syms t;
   T=t*ones(79,95,78);
   RF1=(1/2)*double(int(((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-1/2)),t,0,inf));
   RD1=(3/2)*double(int(((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-3/2)),t,0,inf));
   

2. Numeric method:

   fun1=@(T) ((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-1/2));
   RF1=(1/2)*integral(fun1,0,inf,'ArrayValued',true);
   fun2=@(T) ((T+x1).^(-1/2)).*((T+y1).^(-1/2)).*((T+z).^(-3/2));
   RD1=(3/2)*integral(fun2,0,inf,'ArrayValued',true);
   

其中x1，y1，z是大小为79x95x78的实数矩阵。

- Paul Cai

Answer 1

您正在尝试计算555,750个积分，这需要一定的时间。以下是一些建议，可能会有所帮助。

1. 首先，您可以更高效地表达方程（这可能会影响您的结果数值精度，具体取决于您的值大小和问题的稳定性）：

fun1 = @(T)1./sqrt((T+x1).*(T+y1).*(T+z));
fun2 = @(T)1./sqrt((T+x1).*(T+y1).*(T+z).*(T+z).*(T+z));

一样的转换也可以用于符号情况。在测试中，仅使用较小的（4,524元素数组）x1、y1和z数组的随机值，您的两组积分在我的计算机上计算速度都快了四倍。

2. 将匿名函数转换为常规M文件函数，也可能会使速度略微增加。在单独的文件中：

function RF1 = fun1(T,x1,y1,z)
RF1 = 1./sqrt((T+x1).*(T+y1).*(T+z));

并且

function RD1 = fun2(T,x1,y1,z)
RD1 = 1./sqrt((T+x1).*(T+y1).*(T+z).*(T+z).*(T+z));

然后调用它们，将x1，y1和z作为参数传递：

RF1 = 0.5*integral(@(T)fun1(T,x1,x2,z),0,Inf,'ArrayValued',true);
RD1 = 1.5*integral(@(T)fun2(T,x1,x2,z),0,Inf,'ArrayValued',true);

我在我的机器上看到了一个小的但不可忽略的差异：我尝试的数据大约比正常函数快10%。这可能是因为M-File被单独JIT编译并且常规函数可能会与匿名函数稍微有所不同。

3. 当然，您还应该尝试适当调整integral的绝对和相对容差。您还可以尝试使用向量化的quadv，这将提供更准确的结果并且速度更快（请注意，这个函数将在Matlab的未来版本中删除）。根据数据和函数的性质，甚至可能将问题离散化（使用较大的上界值而不是Inf）并使用trapz之类的方法来执行积分，虽然我不知道这是否会更快。

对于符号计算的情况，您可以尝试将t定义为实数（即syms t real;），并可能指定'IgnoreAnalyticConstraints'选项为true，以查看这些是否能提高速度。符号数学也可能在内存方面效率低下，因此您还可以尝试分解问题并逐行积分。除了int之外，还有其他可能更快的选项，但如果您不需要精度，纯数字方法会快得多（而且没有实际数据值很难确定）。最后，如果您不需要太高的精度，只关心速度，您可以在C mex中实现"快速求幂平方根技巧"。

- horchler