algorithm - 如何最好地优化使用对角线(非对角线)矩阵的 IDL 中的矩阵乘法

Question

我正在寻找最有效的 IDL 代码来替换具有 3 个不同值的特定对角线(非对角线或对角对称)矩阵的 IDL 矩阵乘法 (#) 运算符:对角线上的统一； unity 加上对角线右侧的增量； unity 减去左侧的相同增量。

问题域

IDL(固定；不可协商；抱歉)；无快门 CCD 成像系统上的图像拖尾。

基本问题陈述

给定的

一个 1024x1024 矩阵，“EMatrix”，对角线上具有统一性； (1-delta)
在对角线的左边； (1+delta) 向右；增量 = 0.044。

另一个 1024x1024 矩阵，图像

询问

计算速度最快的 IDL 代码是什么(Image # EMatrix)？

2014-09-16:见下方更新

背景

更大的问题陈述(其中矩阵乘法似乎只是最慢的部分，优化整个例程不会有什么坏处):

http://pdssbn.astro.umd.edu/holdings/nh-j-lorri-3-jupiter-v1.1/document/soc_inst_icd.pdf

第 9.3.1.2 节(PDF 第 47 页；内部第 34 页)，以及同一目录中的其他文档(抱歉作为新手，我只能发布两个链接)

我的工作到目前为止

https://github.com/drbitboy/OptimizeDiagishIDLMatrixMultiply

现在 (2014-09-26) 比 1024x1024 矩阵的 IDL # 运算符快一个数量级。

细节

朴素的操作是 O(n^3) 并执行大约十亿 (2^30) 个 double 乘法和大约相同数量的加法；维基百科进一步告诉我，Strassen 的算法将其降低到 O(n^2.807)，或者对于 n=1024 的乘法约为 282M+。

将其分解为一个简单的 3x3 案例，例如 image 和 EMatrix 是

   image        EMatrix

[ 0  1  2 ]   [ 1  p  p ]
[ 3  4  5 ] # [ m  1  p ]
[ 6  7  8 ]   [ m  m  1 ]

其中 p 代表 1+delta (1.044)，m 代表 1-delta (0.956)。

由于m和p的重复，应该有一个可用的简化:查看图像的中间列，三行的结果应该是

[1,4,7] . [1,p,p] = m*(0)   + 1*1 + p*(4+7)
[1,4,7] . [m,1,p] = m*(1)   + 1*4 + p*(7)
[1,4,7] . [m,m,1] = m*(1+4) + 1*7 + p*(0)

在哪里 。表示点(内部？)乘积，即 [1,4,7].[a,b,c] = (1a + 4b + 7c)

基于此，这是我到目前为止所做的:

中间项只是列本身，乘以 m 和 p 的总和看起来很像列的连续部分的累积总和，可能反转(对于 m)，移位一个并且第一个元素设置为零。

即对于 m:

; shift image down one:
imgminusShift01 = shift(image,0,1)

; zero the top row:
imgminusShift01[*,0] = 0

; Make a cumulative sum:
imgminusCum = total( imageShift01, 2, /cumulative)

对于 p， imgplusShift01Cum 遵循基本相同的路径，但在前后翻转事物的 rotate(...,7) 。

一旦我们有了这三个矩阵(原始图像，及其后代 imgminusShift01Cum 和 imgplusShift01Cum)，所需的结果是

m * imgminusShift01Cum    +   1 * image   +   p * imgplusShift01Cum

为了进行移位和旋转，我使用索引、移位和旋转自己并存储在一个公共(public)块中以供后续调用，这又节省了 10-20%。

总结，到此为止

2014-09-16:见下方更新

这提供了 5+ 的加速。

我期待更多一点，因为我认为我只有 3M 乘法和 2M 加法，所以也许内存分配是昂贵的部分，我应该做类似 REPLICATE_INPLACE 的事情(我的 IDL 生锈了 - 自 6.4 以来我没有做太多和早期的 7.0)。

或者也许 IDL 的矩阵乘法比理论更好？

替代方法和其他想法:

我们可以对 EMatrix 等于 unity 的事实做些什么
加上一个矩阵，对角线上有零，上部有 +/- delta
和较低的三角形？

通过对列求和，我以“错误”的方式顺序访问数据，但是
首先转置 Image 真的会节省时间吗(只有 8MB)？

显然选择另一种语言，获得 GPU 阵列来帮助，或者
编写 DLM，将是其他选择，但让我们将其保留给 IDL
现在。

advTHANKSance(是的，我那么老；-)，

布赖恩·卡西奇 2014-07-20

更新 2014-09-16

Diego 的方法给我带来了一个更简单的解决方案:

我想我们明白了；我的第一遍太复杂了，所有的轮换。

使用迭戈的符号，但转置，我正在寻找

[K] = [IMG] # [E]

由于[IMG]的列乘以[E]的行，因此[IMG]的列之间没有交互作用，所以为了分析我们只需要查看[IMG]的一列，其中的点(内)积，与[E]的行，成为结果[K]的一列。将这个想法扩展到包含元素 x、y 和 z 的 3x3 矩阵的一列:

[Kx]   [x]   [1    1+d  1+d]
[Ky] = [y] # [1-d  1    1+d]
[Kz] = [z]   [1-d  1-d  1  ]

特别查看上面的元素 Ky( [K] 的一个 元素，对应于 [IMG] 中的 y，分解为仅使用标量的公式):

Ky = x * (1-d)  +  y * 1  + z * (1+d)

对任意长度列中的任意 y 进行泛化:

Ky = sumx * (1-d)  +  y * 1  + sumz * (1+d)

其中标量 sumx 和 sumz 分别是 [IMG] 的该列中所有高于和低于 y 的值的总和。注意 sumx 和 sumz 是特定于元素 y 的值。

重新排列:

Ky = (sumx + y + sumz) - d * (sumx - sumz)

Ky = tot - d * (sumx - sumz)

在哪里

tot = (sumx + y + sumz) 

即，tot 是列中所有值的总和(例如，在 IDL 中:tot = total(IMG,2))。

所以到目前为止，我基本上复制了 Diego 的工作；此分析的其余部分将 Ky 的最后一个方程转换为适合在 IDL 中进行快速评估的形式。

求解 sumz 的 tot 方程:

sumz = tot - (y + sumx)

代入 Ky:

Ky = tot - (sumx - (tot - (y + sumx)))

Ky = tot - ((2 * sumx) + y - tot)

Ky = tot + (tot - ((2 * sumx) + y)

使用 sumxy 表示从上到下的列中所有值的总和，包括 y (IDL: [SUMXY] = total([IMG],2,/CUMULATIVE))

sumxy = sumx + y

和

sumx = sumxy - y

代入 Ky:

Ky = tot + (tot - ((2 * (sumxy - y)) + y)

Ky = tot + (tot + y - (2 * sumxy))

因此，如果我们可以为 [IMG] 的每个元素计算 tot 和 sumxy，即如果我们可以计算矩阵 [TOT] 和 [SUMXY]，并且我们已经有了 [IMG] 作为 y 的矩阵版本，那么它是一个简单的这些矩阵的线性组合。

在 IDL 中，这些很简单:

[SUMXY] = TOTAL([IMG],2,/CUMULATIVE)

[TOT] = [SUMXY][*,N-1] # REPLICATE(1D0,1,N)

IE。 [TOT] 是 [SUMXY] 的最后一行，复制形成 N 行矩阵。

最终代码如下所示:

function lorri_ematrix_multiply,NxN,EMatrix

  NROWS = (size(NxN,/DIM))[1]

  SUMXY = TOTAL(NxN,2,/CUMULATIVE)

  TOT   = SUMXY[*,NROWS-1] # REPLICATE(1,NROWS,1d0)

  RETURN, TOT + ((EMatrix[1,0] - 1d0) * (TOT + NxN - (2d0 * SUMXY)))

end

在我们的系统上，它比 [IMG] # [E] 快一个数量级。

注意delta = (EMatrix[1,0] - 1d0)

呜呼!

Answer 1

第1步:

我直接使用数学符号，因为我认为这比用文字解释更清楚:

      [ +1 +d +d +d +d ]   [ 1 0 0 0 0 ]       [ 0 1 1 1 1 ]        [ 0 0 0 0 0 ]
      [ -d +1 +d +d +d ]   [ 0 1 0 0 0 ]       [ 0 0 1 1 1 ]        [ 1 0 0 0 0 ]
[E] = [ -d -d +1 +d +d ] = [ 0 0 1 0 0 ] + d * [ 0 0 0 1 1 ] - d *  [ 1 1 0 0 0 ]
      [ -d -d -d +1 +d ] | [ 0 0 0 1 0 ]       [ 0 0 0 0 1 ]        [ 1 1 1 0 0 ]
      [ -d -d -d -d +1 ] | [ 0 0 0 0 1 ]       [ 0 0 0 0 0 ]        [ 1 1 1 1 0 ]
                         |
                         | [ 1 0 0 0 0 ]       [ 1 1 1 1 1 ]        [ 1 0 0 0 0 ]
                         | [ 0 1 0 0 0 ]       [ 0 1 1 1 1 ]        [ 1 1 0 0 0 ]
                         = [ 0 0 1 0 0 ] + d * [ 0 0 1 1 1 ] - d *  [ 1 1 1 0 0 ]
                         | [ 0 0 0 1 0 ]       [ 0 0 0 1 1 ]        [ 1 1 1 1 0 ]
                         | [ 0 0 0 0 1 ]       [ 0 0 0 0 1 ]        [ 1 1 1 1 1 ]
                         |     [ID]                [UT]                 [LT]
                         |
                         = [ID] + d * [UT] - d * [LT]

==>

[Img] # [E] = [E]##[Img] = [Img] + d * [UT] ## [Img] - d * [LT] ## [Img]

现在让我们观察一下什么是 [LT] ## [Img] :

第 1 行与 [Img] 的第一行相同

第 2 行是 [Img]

的第 1 行和第 2 行的(按列)总和

第 i 行是 [Img]

的所有前 i 行的(按列)总和

第 n 行是 [Img]

的所有行的(按列)总和

所以计算它的有效方法是:

TOTAL(Image, 2, /CUMULATIVE)

类似但有点不同的是 [UT] ## [Img] 的结果:

第 1 行是 [Img]

的所有行的(按列)总和

第 i 行是 [Img]

的所有最后 i 行的(按列)总和

第 n-1 行是 [Img]

的第 1 行和第 2 行的(按列)总和

第 n 行与 [Img]

的第一行相同

所以 [UT] ## [Img] = REVERSE(TOTAL(REVERSE(Image,2), 2, /CUMULATIVE),2)
然后:

[Img] # [E] = [E]##[Img] = Image + d * (REVERSE(TOTAL(REVERSE(Image,2), 2, /CUMULATIVE),2) - TOTAL(Image, 2, /CUMULATIVE))

请注意，我们看到最终每列的结果仅来自同一列的数据。

第2步:

现在让我们看看 [K] = [UT] ## [Img] - [LT] ## [Img]看看它的样子。如果对于每个通用列，我们将列元素命名为 r(1), r(2), r(3), .... ,r(i), ..., r(n) 我们可以看到相应的 [ K] 列元素 R(i)
看起来像这样:

当 n 为偶数时 (1024就是这种情况)

row 1     => R(1) = +r(1)                 -r(1) -r(2) .... -r(n-1) -r(n) = -SUM(j=2, n, r(n))
row 2     => R(1) = +r(1) +r(2)           -r(1) -r(2) .... -r(n-1)       = -SUM(j=3, n-1, r(n))
row 3     => R(1) = +r(1) +r(2) +r(3)     -r(1) -r(2) .... -r(n-2)       = -SUM(j=4, n-2, r(n))
    :        :       :     :     :    :    :     :     :    :               :
row i (i < n/2) 
          => R(1) = +r(1) ... +r(i)       -r(1) -r(2) .... -r(n-i+1)     = -SUM(j=i+1, n-i+1, r(n))
    :        :       :     :     :    :    :     :     :    :               :
row n/2   => R(1) = +r(1) ... +r(n/2)     -r(1) -r(2) .... -r(n/2+1)     = -r(n/2+1) 
row n/2+1 => R(1) = +r(1) ... +r(n/2+1)   -r(1) -r(2) .... -r(n/2)       = +r(n/2+1) 
    :        :       :     :     :    :    :     :     :    :               :
row i (i > n/2)
          => R(1) = +r(1) ... + r(i)      -r(1) -r(2) .... -r(n-i+1)     = +SUM(j=n-i+2, i, r(n))
                                                                         = -R(n-i+1)
    :        :       :     :     :    :    :     :     :    :               :
row n     => R(1) = +r(1) ... + r(n)      -r(1)                          = +SUM(j=2, n, r(n))
                                                                         = -R(1)

当 n 为奇数时
类似但 R((n+1)/2)将全部为 0。我不会详细讨论这个。

重要的是矩阵 [K] = [UT] ## [Img] - [LT] ## [Img]相对于它的水平二等分线是反对称的。

这意味着我们只能计算半矩阵中的值(假设顶部)，然后通过镜像和更改符号来填充下部。
请注意，顶部的有效计算可以从 R(n/2) = r(n/2+1) 开始，然后向上减小索引 (R(n/2 -1), R(n/2 -2), R(n/2 -3)...)每次使用 R(i) = R(i+1) - r(i+1) - r(n-i+1)可以很好地重写为 R(i-1) = R(i) - r(i) - r(n-i+2) .

作为计算的问题，这大约完成了一半的计算，但作为实际速度的问题，需要对其进行测试，以查看使用显式操作的实现是否与内置函数的内部实现一样快，如 TOTAL(/CUMULATIVE)和类似的。很有可能它更快，因为我们在这里也可以避免 TRANSPOSE 和/或 REVERSE。

让我们知道它是如何进行一些分析的!