Bri*_*ich 5 algorithm optimization idl-programming-language matrix-multiplication
我正在寻找最有效的IDL代码来替换IDL矩阵乘(#)运算符,用于特定的,对角线(非对角线或对角线)矩阵,具有3个不同的值:对角线上的单位; 统一加上对角线右边的三角形; 统一减去左边相同的delta.
IDL(已修复;不可协商;对不起); 快门式CCD成像系统上的图像模糊.
一个1024x1024的矩阵,"EMatrix",对角线上有一个整体; 对角线左边的(1-delta); (1 + delta)向右; delta = 0.044.
另一个1024x1024矩阵,Image
2014-09-16:请参阅下面的更新
更大的问题陈述(其中矩阵乘法似乎只是最慢的部分,优化整个例程不会受到影响):
第9.3.1.2节(PDF第47页;内部页面34)和同一目录中的其他文件(抱歉作为新手我只能发布两个链接)
现在(2014-09-26)比1024x1024矩阵的IDL#运算符快一个数量级.
天真的操作是O(n ^ 3)并执行大约十亿(2 ^ 30)双精度乘法和大约相同数量的加法; 维基百科进一步告诉我,斯特拉森的算法将其降低到O(n ^ 2.807),或者约为282M +,n = 1024.
对于一个简单的3x3案例进行细分,比如图像和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),移位1并且第一个元素设置为零.
即对于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基本上遵循相同的路径,但是在向上和向下翻转之前和之后使用旋转(...,7).
一旦我们有了这三个矩阵(原始图像,其后代imgminusShift01Cum和imgplusShift01Cum),所需的结果是
m * imgminusShift01Cum + 1 * image + p * imgplusShift01Cum
为了进行移位和旋转,我使用索引,移位并自行旋转并存储在公共块中以用于后续调用,从而节省另外10-20%.
2014-09-16:请参阅下面的更新
这样可以加速5+.
我期待更多,因为我认为我已经下降到3M倍增和2M添加,所以也许内存分配是昂贵的部分,我应该做一些像REPLICATE_INPLACE(我的IDL生锈 - 我自6.4以来没有做太多和早期的7.0).
或者也许IDL的矩阵乘法比理论更好？
替代方法和其他想法:
我们可以做一些事情,即EMatrix等于单位加上一个矩阵,对角线上的零点和上三角形和下三角形的+/-三角形？
通过对列进行求和,我以"错误"的方式顺序访问数据,但它实际上是否会节省时间来首先转置Image(它只有8MB)？
显然,选择另一种语言,获得GPU阵列来帮助,或者编写DLM,将是其他选择,但现在让我们将其保留给IDL.
advTHANKSance(是的,我老了;-),
Brian Carcich 2014-07-20
我想我们得到了它; 我的第一次传球太复杂了,所有的轮换.
要使用迭戈的符号,但转置,我正在寻找
[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和y以上所有值的总和.NB 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)).
所以到目前为止我基本上复制了迭戈的工作; 该分析的其余部分将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)
即[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]快一个数量级.
NB delta =(EMatrix [1,0] - 1d0)
呜啊!
我直接使用数学符号，因为我认为这比用文字解释更清楚：
[ +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))
请注意，我们看到最终每列的结果仅来自同一列的数据。
现在让我们看看[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。
让我们通过一些分析了解一下它是如何进行的！