是否存在一种算法,可将线性索引转换为支持负步长的子脚本列表?
背景
环境(如MATLAB、Julia等)和库(如NumPy)提供对跨步数组(也称为ndarrays)的支持。跨步数组由线性内存支持(例如单个基础缓冲区),这与嵌套数组形成对比,其中每个嵌套数组对应一个维度。例如,考虑以下2x2矩阵。
[ 1 2
3 4 ]
实现为一个数组的数组。
A = [ [ 1, 2 ], [ 3, 4 ] ]
在哪里(使用从零开始的索引)
a01 = A[0][1] = 2
我们可以将相同的2x2矩阵表示为步幅数组,如下所示(假设行优先)。
A = [ 1, 2,
3, 4 ]
在哪里
a01 = A[ 2*0 + 1*1 ] = 2
通常情况下,对于一个步长为 NxM 的矩阵,元素
(i,j)
可以通过以下方式访问。function get( i, j ) {
return buffer[ si*i + sj*j ];
}
这里的 buffer
是底层数据缓冲区,si
和 sj
分别对应于沿 i
和 j
维度的步长。假设按行主序分步的数组,对于上面的 2x2 矩阵,si = 2
,sj = 1
(省略元素字节长度)。
一般来说,步长可以根据数组形状计算如下:
function shape2strides( shape, order ) {
var out = new Array( shape.length );
var s = 1;
var i;
if ( order === 'column-major' ) {
for ( i = 0; i < shape.length; i++ ) {
out[ i ] = shape[ i ];
s *= shape[ i ];
}
return out;
} else { // row-major
for ( i = shape.length-1; i >= 0; i-- ) {
out[ i ] = shape[ i ];
s *= shape[ i ];
}
}
}
为了方便处理步进数组,环境/库通常提供便利函数,允许在线性索引和下标之间轻松转换。例如,在 MATLAB 中,要从下标转换为线性索引。
idx = sub2ind( size( A ), i, j )
同样地,要将线性索引转换为MATLAB中的下标
s = ind2sub( size( A ), idx )
Julia还具有sub2ind和ind2sub功能。在NumPy中,您可以使用ravel_multi_index和unravel_index。
除了数据局部性外,跨步数组很方便,因为它们允许通过操作跨度是否为负来创建数组“视图”。当步幅为负时,我们沿该维从右到左迭代,而不是从左到右迭代。为了支持这种迭代行为,我们需要确定底层数据缓冲区中的第一个索引元素在哪里。按照惯例,我们将称此索引为“偏移量”,其计算方法如下
function strides2offset( shape, strides ) {
var offset = 0;
var i;
for ( i = 0; i < shape.length; i++ ) {
if ( strides[ i ] < 0 ) {
offset -= strides[i] * ( shape[i]-1 ); // increments the offset
}
}
return offset;
}
一旦我们有了偏移量,我们需要修改我们的
get( i, j )
函数如下:function get( i, j ) {
return buffer[ offset + si*i + sj*j ];
}
对于步幅为
2,1
的2x2矩阵A,偏移量为0
,因此返回原始的get
函数。当步幅为2,-1
时,偏移量为1
;对于-2,1
,偏移量为2
;对于-2,-1
,偏移量为3
。因此,我们可以生成以下矩阵视图(假设按行主序)。Dims: 2x2
Strides: 2,1
Offset: 0
A = [ 1, 2,
3, 4 ]
Strides: 2,-1
Offset: 1
A = [ 2, 1,
4, 3 ]
Strides: -2,1
Offset: 2
A = [ 3, 4,
1, 2 ]
Strides: -2,-1
Offset: 3
A = [ 4, 3,
2, 1 ]
上述观点突出了步幅数组的优势之一:O(1)操作。例如,要将矩阵从左到右翻转,我们只需要翻转第二维度的步幅符号(假设按行主序)。要将其从上到下翻转,我们需要翻转第一维度的步幅符号(假设按行主序)。要同时进行左右和上下翻转,我们需要翻转两个步幅符号。所有上述操作都不涉及触碰底层数据缓冲区;我们只是更改步幅数组元数据。
sub2ind
即使考虑到负步幅(即步幅数组视图),从下标转换为线性索引也很简单。例如,对于任意维度的步幅数组,
function sub2ind( ...subscripts ) {
var sub;
var idx;
var s;
var n;
idx = offset;
for ( n = 0; n < shape.length; n++ ) {
sub = subscripts[ n ];
s = strides[ n ];
if ( s < 0 && offset === 0 ) { // assume want "view" index
idx -= sub * s; // always increments `idx`
} else { // assume want underlying data buffer index
idx += sub * s; // may increment or decrement `idx`
}
}
return idx;
}
在这里,我们允许从视图或基础数据缓冲区的角度返回线性索引。当“偏移量”为
0
时,我们假设始终返回视图中的线性索引(这可能与基础数据缓冲区中的线性索引不对应)。换句话说,对于2x2矩阵视图,(0,0) => 0,(0,1) => 1,(1,0) => 2,(1,1) => 3
,总是如此。从使用视图的角度来看,这种映射符合直觉。当我想要A(0,0)
时,我希望该元素位于“第一个”线性索引处,即使它实际上并不存储在基础数据缓冲区中。您可以自行证明,在将元素查找扩展到负步长时,
sub2ind
会像上面详细说明的那样为任何偏移量返回相同的索引。例如实现,请参见Julia, NumPy, 和 stdlib。
ind2sub
被问及的问题是如何实现
sub2ind
的反向操作,并支持负步长。对于正向步幅(因此,偏移量为
0
),我们可以使用模算术来恢复下标。例如,考虑解决NxMxL步进数组的线性索引方程。idx = offset + si*i + sj*j + sk*k
假设按行主序排列,其中
si = nj*nk, sj = nk, sk = 1
,ni, nj, nk
分别为维度大小 N, M, L
。代入数值,idx = 0 + (nj*nk)*i + nk*j + k
可以重新排列
idx = nk*(nj*i + j) + k
如果我们使用
nk
对两边取模,则有:idx % nk = k
知道 k
,让我们重新排列初始方程式
(idx - k) = nk*(nj*i + j)
(idx - k)/nk = nj*i + j
如果我们使用
nj
对两边取模,那么:((idx - k)/nk) % nj = j
了解变量 j
后,我们可以重新排列初始方程以解出变量 i
(((idx - k)/nk) - j)/nj = i
上述算法可以推广到任意维度,并且很容易实现(也可以参考Julia和NumPy)。
function ind2sub( idx, order ) {
var out = new Array( shape.length );
var s;
var i;
if ( order === 'column-major' ) {
for ( i = 0; i < shape.length; i++ ) {
s = idx % shape[ i ];
idx -= s;
idx /= shape[ i ];
out[ i ] = s;
}
} else { // row-major
for ( i = shape.length-1; i >= 0; i-- ) {
s = idx % shape[ i ];
idx -= s;
idx /= shape[ i ];
out[ i ] = s;
}
}
return out;
}
上述算法使用模算术,但不支持负步长。如果我们使用相同的过程来解决下标,则需要从以下方程开始:
idx = offset + nk*(nj*i + j) + k
可以简化为
idx-offset = nk*(nj*i + j) + k
这里的问题当然是idx-offset
可能为负数,有效地移动了可能的i,j,k
值的范围(i
应该在半开区间[0,N)上;j
在区间[0,M)上;和k
在区间[0,L)上)。
这引发了一个问题,即是否存在一种算法,用于将线性索引转换为支持负步长的下标。或者换句话说,是否有一种算法,可以在给定底层数据缓冲区中的线性索引的情况下,返回相应的视图下标?
其他语言/库中的实现(例如Julia和NumPy)似乎只支持offset = 0
的情况。我正在寻找更通用的东西,可以适用于跨步数组视图。
如果有现有实现/算法的任何指针,将不胜感激。
numpy.lib.stride_tricks.as_strided
。也许吧。整数倍的维度步幅可以工作。我不认为负步幅会起作用,但as_strided
会创建一个视图,你可以通过使用花式索引来创建该视图的视图 -view[::-1]
。 - Daniel Fnumpy
更低级的语言,因为这些内存级算法通常是用一些低到中级语言如C
或FORTRAN
实现的。 - Daniel Fnumpy
是SO推荐的,所以我选择了它。我可能会在明天或后天更新标签。感谢您的建议。 - kgryte