我正在寻找最有效的 IDL 代码来替换具有 3 个不同值的特定对角线(非对角线或对角对称)矩阵的 IDL 矩阵乘法 (#) 运算符:对角线上的统一； unity 加上对角线右侧的增量； unity 减去左侧的相同增量。
问题域
IDL(固定；不可协商；抱歉)；无快门 CCD 成像系统上的图像拖尾。
基本问题陈述
给定的
在对角线的左边； (1+delta) 向右；增量 = 0.044。
询问
2014-09-16:见下方更新
背景
更大的问题陈述(其中矩阵乘法似乎只是最慢的部分，优化整个例程不会有什么坏处):
第 9.3.1.2 节(PDF 第 47 页；内部第 34 页)，以及同一目录中的其他文档(抱歉作为新手，我只能发布两个链接)
我的工作到目前为止
现在 (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 的矩阵乘法比理论更好？
替代方法和其他想法:
加上一个矩阵，对角线上有零，上部有 +/- delta
和较低的三角形？
首先转置 Image 真的会节省时间吗(只有 8MB)？
编写 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)
呜呼!
最佳答案
第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] :所以计算它的有效方法是:
TOTAL(Image, 2, /CUMULATIVE)
类似但有点不同的是
[UT] ## [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。让我们知道它是如何进行一些分析的!
关于algorithm - 如何最好地优化使用对角线(非对角线)矩阵的 IDL 中的矩阵乘法，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/24854289/