如何将线性索引转换为支持负步长的下标

Question

是否存在一种算法，可将线性索引转换为支持负步长的子脚本列表？

背景

环境（如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的情况。我正在寻找更通用的东西，可以适用于跨步数组视图。

如果有现有实现/算法的任何指针，将不胜感激。

- kgryte

---

抱歉，我看到那一大堆文字就有些恍惚了，但我认为你正在寻找numpy.lib.stride_tricks.as_strided。也许吧。整数倍的维度步幅可以工作。我不认为负步幅会起作用，但as_strided会创建一个视图，你可以通过使用花式索引来创建该视图的视图 - view[::-1]。 - Daniel F

@DanielF 感谢您的评论，但不是我要找的。正如 OP 中所述，我对一种泛化到负步长的算法感兴趣。理想情况下，这个算法应该是与编程语言/库无关的。你提出的建议与NumPy密切相关。 - kgryte

啊，我明白了。我想你可能需要标记一些比numpy更低级的语言，因为这些内存级算法通常是用一些低到中级语言如C或FORTRAN实现的。 - Daniel F

@DanielF 是的，numpy 是SO推荐的，所以我选择了它。我可能会在明天或后天更新标签。感谢您的建议。 - kgryte

Answer 1

（编辑 - 我可能正在处理更简单的nd索引到平面索引的情况，而您则专注于反向。现在探索这个任务已经太晚了 - 我明天会重新审视这个问题。）

---

如果偏移量正确，我认为相同的公式，用于将n-d索引转换为平面索引，可以适用于负步幅和正步幅：

比较一个（3,4）数组的索引与其双重翻转：

In [32]: x = np.arange(12).reshape(3,4)
In [33]: y = x[::-1, ::-1]
In [34]: x.strides
Out[34]: (16, 4)
In [35]: y.strides
Out[35]: (-16, -4)

数据缓冲区“start”可以在__array_interface__中找到：

In [36]: x.__array_interface__['data']
Out[36]: (166934688, False)
In [37]: y.__array_interface__['data']
Out[37]: (166934732, False)
In [38]: 732-688
Out[38]: 44

缓冲区中有48个字节，但y的偏移量为44，即x [2,3]元素（在此示例中为11）的开头。

现在测试一个x元素的平面索引：

In [39]: x[1,2]
Out[39]: 6            # value
In [40]: x[1:2, 2:3].__array_interface__['data']    # a view
Out[40]: (166934712, False)
In [41]: 688+(1*16)+(2*4)    # offset + si*i + sj*j 
Out[41]: 712

现在对y做同样的操作：

In [42]: y[1:2, 2:3].__array_interface__['data']
Out[42]: (166934708, False)
In [43]: 732+(1*-16)+(2*-4)
Out[43]: 708

- hpaulj

Answer 2

我写了一段代码，我认为可以解决你的问题。除了你要求的那个函数之外，还有其他函数略微不同，所以我把它们都放在这里，以避免潜在的不一致性。

以下代码没有用任何特定的语言编写。但是，你可以在其中找到一些C语法元素。

function calcStrides(shape[]) {
    strides[]; # Result array
    currentStride = 1;
    for(i = shape.size; 0 < i;) {
        --i;
        if(0 < shape[i]) {
            strides[i] = currentStride;
            currentStride *= shape[i];
        } else {
            strides[i] = -currentStride;
            currentStride *= -shape[i];
        }
    }
    return strides;
}

function calcOffset(shape[], strides[]) {
    offset = 0;
    for(i = 0; i < shape.size; ++i) {
        if(shape[i] < 0) {
            offset += strides[i] * (shape[i] + 1);
        }
    }
    return offset;
}

function sub2ind(strides[], offset, subs[]) {
    ind = offset;
    for(i = 0; i < strides.size; ++i) {
        ind += strides[i] * subs[i];
    }
    return ind;
};

function ind2sub(shape[], strides[], ind) {
    subs[]; # Result array
    for(i = 0; i < shape.size; ++i) {
        if(0 < strides[i]) {
            subs[i] = ind / strides[i];
            ind -= subs[i] * strides[i];
        } else {
            absSub = ind / -strides[i];
            subs[i] = -shape[i] - 1 - absSub;
            ind -= absSub * -strides[i];
        }
    }
    return subs;
}

- navyblue

---

感谢您的回答。我可以确认 ind2sub 有以下注意事项：（1）在 ind2sub 中的 / 是整数除法（或截断浮点除法）。 （2）shape[i] 永远不为负数，正如在 calcStrides 和 calcOffset 中所示。 shape[] 中的每个元素对应于维度长度，其始终为非负数。（3）鉴于 shape [i] >= 0，在 ind2sub 的 else 块中，应该是 shape [i]，而不是 -shape [i]。 经过这些更改，我可以确认计算出了针对 2x2、3x3 和 2x2x2 分步数组的预期下标。 - kgryte

(1) 这是整数除法，我认为没有任何问题。 (2) & (3) 我的假设与你略有不同。在我的版本中，形状不仅存储数组范围的信息，还存储数组遍历的方向。但如果我们采用你的方法，我同意你的观点。这些差异来自于我实际上是从一个更大的解决方案中剪切出了这段代码。 - navyblue

关于我之前的评论（我无法再编辑它），特别是第（1）点，如果结果没有朝零舍入，或者您无法控制舍入模式，整数除法可能会成为问题。 - navyblue

是的，你在整数除法方面是正确的。因为我在本地实现时使用了浮点除法并截断，所以结果总是朝零四舍五入。感谢您抽出时间提供答案。我可能会发布一个基于您提供的解决方案，并与 OP 中描述的函数一致的解决方案。 - kgryte