在MATLAB中加速exp(A*x)的分析方法

ber*_*ers 8 performance matlab simplify exp

我需要f(x)=exp(A*x)反复计算一个微小的,可变的列向量x和一个巨大的,恒定的矩阵A(许多行,几列).换句话说,x很少,但A*x很多.我的问题维度是这样的,它A*x占用了与exp()部分一样多的运行时间.

除了泰勒展开,并预先计算值的范围exp(y)(假设已知的范围内y的值A*x),这是我没有设法大大加快(同时保持精度)就什么MATLAB是对自己做的,我考虑分析性地重述问题,以便能够预先计算某些值.

例如,我发现了 exp(A*x)_i = exp(\sum_j A_ij x_j) = \prod_j exp(A_ij x_j) = \prod_j exp(A_ij)^x_j

这将允许我预先计算exp(A)一次,但是循环中所需的取幂与原始exp()函数调用一样昂贵,并且必须另外执行乘法(\ prod).

我还有其他想法可以遵循,或者我可能错过了MATLAB中的解决方案吗?

编辑:更多细节

A尺寸为81的是26873856(是的,那是巨大的),所以x是81乘1. nnz(A) / numel(A)0.0012,nnz(A*x) / numel(A*x)0.0075.我已经使用稀疏矩阵来表示A,但稀疏矩阵的exp()不再稀疏.所以事实上,我存储x非稀疏和我计算的exp(full(A*x))结果是快/慢full(exp(A*x))(我认为无论A*x是非稀疏的,因为x是非稀疏的.)exp(full(A*sparse(x)))是一种稀疏的方法A*x,但速度较慢.甚至更慢的变体exp(A*sparse(x))(对于稀疏类型的非稀疏矩阵具有双倍的内存影响)和full(exp(A*sparse(x))(其再次产生非稀疏结果).

sx = sparse(x);
tic, for i = 1 : 10, exp(full(A*x)); end, toc
tic, for i = 1 : 10, full(exp(A*x)); end, toc
tic, for i = 1 : 10, exp(full(A*sx)); end, toc
tic, for i = 1 : 10, exp(A*sx); end, toc
tic, for i = 1 : 10, full(exp(A*sx)); end, toc

Elapsed time is 1.485935 seconds.
Elapsed time is 1.511304 seconds.
Elapsed time is 2.060104 seconds.
Elapsed time is 3.194711 seconds.
Elapsed time is 4.534749 seconds.
Run Code Online (Sandbox Code Playgroud)

是的,我确实计算了元素明确的exp,我更新上面的等式来反映这一点.

还有一个编辑:我试图变得聪明,但收效甚微:

tic, for i = 1 : 10, B = exp(A*x); end, toc
tic, for i = 1 : 10, C = 1 + full(spfun(@(x) exp(x) - 1, A * sx)); end, toc
tic, for i = 1 : 10, D = 1 + full(spfun(@(x) exp(x) - 1, A * x)); end, toc
tic, for i = 1 : 10, E = 1 + full(spfun(@(x) exp(x) - 1, sparse(A * x))); end, toc
tic, for i = 1 : 10, F = 1 + spfun(@(x) exp(x) - 1, A * sx); end, toc
tic, for i = 1 : 10, G = 1 + spfun(@(x) exp(x) - 1, A * x); end, toc
tic, for i = 1 : 10, H = 1 + spfun(@(x) exp(x) - 1, sparse(A * x)); end, toc

Elapsed time is 1.490776 seconds.
Elapsed time is 2.031305 seconds.
Elapsed time is 2.743365 seconds.
Elapsed time is 2.818630 seconds.
Elapsed time is 2.176082 seconds.
Elapsed time is 2.779800 seconds.
Elapsed time is 2.900107 seconds.
Run Code Online (Sandbox Code Playgroud)

Eng*_*ica 2

计算机实际上并不计算指数。您可能会认为它们会这样做,但它们所做的是高精度多项式近似。

参考:

最后一个参考看起来相当不错。也许它应该是第一位的。

由于您正在处理图像,因此您可能有离散数量的强度级别(通常为 255)。这可以减少采样或查找,具体取决于“A”的性质。检查这一点的一种方法是对一组具有足够代表性的“x”值执行类似以下操作:

y=Ax
cdfplot(y(:))
Run Code Online (Sandbox Code Playgroud)

如果您能够将图像预先分割为“更有趣”和“不太有趣”——就像您正在看 X 射线一样,能够修剪掉所有“人体外部”位置并将它们夹在零来预先稀疏您的数据,这可以减少唯一值的数量。您可以考虑数据中每个独特“模式”的前一个。

我的方法包括:

  • 查看 exp(x) 的替代公式,其精度较低但速度较高
  • 如果“x”的级别足够少,请考虑表查找
  • 如果您有“稍微太多”的级别来进行表查找,请考虑插值和表查找的组合
  • 考虑基于分段模式的单个查找(或替代公式)。如果您知道它是骨头并且正在寻找静脉,那么也许应该应用较少的高成本数据处理。

现在我必须问自己为什么你会生活在 exp(A*x)*x 的这么多迭代中,我认为你可能会在频率/波数域和时间/空间域之间来回切换。您还可能会使用 exp(x) 作为基础来处理概率,并进行一些贝叶斯乐趣。我不知道 exp(x) 是一个很好的共轭先验,所以我将使用傅立叶材料。

其他选项: - 考虑使用 fft、fft2 或 fftn 给定您的矩阵 - 它们速度很快,并且可能可以满足您正在寻找的部分功能。

我确信以下内容存在前域变化:

您也许可以使用 woodbury 矩阵将查找与计算混合在一起。不过,我必须考虑一下这一点才能确定。(链接)一度我知道所有重要的事情(CFD、FEA、FFT)都与矩阵求逆有关,但后来我忘记了具体的细节。

现在,如果您生活在 MatLab 中,那么您可能会考虑使用“编码器”将 MatLab 代码转换为 C 代码。不管解释器有多么有趣,一个好的 C 编译器可以更快。我使用的助记符(希望不要太雄心勃勃)如下所示:链接从 13:49 左右开始。它非常简单,但它显示了规范解释语言(python)和相同语言的编译版本(cython/c)之间的区别。

我确信,如果我有一些更具体的信息,并且被要求这样做,那么我可以更积极地参与更具体的相关答案。

您可能没有在传统硬件上执行此操作的好方法,您可能会考虑购买诸如 GPGPU 之类的东西。CUDA 及其同类产品具有大规模并行操作,可以以几块显卡的成本大幅提高速度。你可以有数千个“核心”(夸张的管道)来完成几个 ALU 的工作,如果这项工作可以适当地并行化(就像这样),那么它可以更快地完成。

编辑:

我在考虑Eureqa。如果我有一些“大铁”用于开发而不是生产,我会考虑的一个选择是使用他们的 Eureqa 产品来得出足够快、足够准确的近似值。

如果您对“A”矩阵执行“快速”奇异值分解,您会发现主要性能由 81 个特征向量控制。我会查看特征值,看看是否只有这 81 个特征向量中的少数几个提供了大部分信息。如果是这种情况,那么您可以将其他值限制为零,并构造一个简单的转换。

现在,如果是我,我会想从指数中得到“A”。我想知道您是否可以查看 81x81 特征向量矩阵和“x”,并思考一下线性代数,以及您将向量投影到什么空间。有什么方法可以创建如下所示的函数:

f(x) = B2 * exp( B1 * x )

使得

B1*x

比你现在的等级低得多

斧头。