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

我正在寻找最有效的IDL代码来替换IDL矩阵乘(#)运算符,用于特定的,对角线(非对角线或对角线)矩阵,具有3个不同的值:对角线上的单位; 统一加上对角线右边的三角形; 统一减去左边相同的delta.

IDL(已修复;不可协商;对不起); 快门式CCD成像系统上的图像模糊.

特定

一个1024x1024的矩阵,"EMatrix",对角线上有一个整体; 对角线左边的(1-delta); (1 + delta)向右; delta = 0.044.

另一个1024x1024矩阵,Image

询问

什么是最快的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)双精度乘法和大约相同数量的加法; 维基百科进一步告诉我,斯特拉森的算法将其降低到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)

呜啊!