PyTorch 的 torch.as_strided 在制作 Toeplitz 矩阵方面取得了负面进展

PyTorch's torch.as_strided with negative strides for making a Toeplitz matrix

我正在编写 scipy.linalg.toeplitz 的临时装配 PyTorch 版本,目前具有以下形式:

def toeplitz_torch(c, r=None):
    c = torch.tensor(c).ravel()
    if r is None:
        r = torch.conj(c)
    else:
        r = torch.tensor(r).ravel()

    # Flip c left to right.
    idx = [i for i in range(c.size(0)-1, -1, -1)]
    idx = torch.LongTensor(idx)
    c = c.index_select(0, idx)

    vals = torch.cat((c, r[1:]))
    out_shp = len(c), len(r)
    n = vals.stride(0)

    return torch.as_strided(vals[len(c)-1:], size=out_shp, stride=(-n, n)).copy()

但是torch.as_strided目前不支持负步幅。因此,我的函数抛出错误:

RuntimeError: as_strided: Negative strides are not supported at the moment, got strides: [-1, 1].

我对 as_strided 的(可能不正确的)理解是,它将第一个参数的值插入到一个新数组中,该数组的大小由第二个参数指定,并且它通过在原始数组并将它们放置在最后一个参数给出的下标索引步长处。

关于 as_strided 的 NumPy 和 PyTorch 文档都有关于使用该函数时“非常小心”的可怕警告,我不完全理解这个函数,所以我想问一下:

  1. 我对as_strided的理解正确吗?
  2. 是否有一种简单的方法可以重写此内容以便负步幅有效?
  3. 我可以通过 toeplitz_torch 传递渐变 w.r.t c(或 r)吗?

> 1.我对as_strided的理解是否正确?

步幅是张量访问底层连续数据缓冲区的接口。它不会 插入 值,torch.as_strided 不会复制值,步长定义了我们称之为多维数组(在 NumPy 中)的人工布局或张量(在 PyTorch 中)。

Andreas K. puts it in another :

Strides are the number of bytes to jump over in the memory in order to get from one item to the next item along each direction/dimension of the array. In other words, it's the byte-separation between consecutive items for each dimension.

如果您在步幅方面遇到问题,请随时阅读那边的答案。在这里,我们将以您的示例为例,看看它是如何使用 as_strided.

实现的

example given by Scipy for linalg.toeplitz如下:

>>> toeplitz([1,2,3], [1,4,5,6])
array([[1, 4, 5, 6],
       [2, 1, 4, 5],
       [3, 2, 1, 4]])

为此,他们首先构造值列表(我们可以称之为基础值,实际上不是基础数据):vals作为 [3 2 1 4 5 6] Toeplitz 列和行变平。

现在注意传递给 np.lib.stride_tricks.as_strided 的参数:

  • valuesvals[len(c)-1:] 注意切片:张量显示较小,但基本值仍然存在,并且它们对应于 vals 的值。继续用storage_offset比较两者:只是偏移了2,数值还在!这是如何工作的,它基本上改变了索引,使得 index=0 将引用值 1index=14,等等......

  • shape:由 column/row 输入给出,此处为 (3, 4)。这是生成的对象的形状。

  • strides:这是最重要的部分:(-n, n),在本例中为 (-1, 1)

用步幅做最直观的事情是描述多维space:(i, j) ∈ [0,3[ x [0,4[和扁平一维space:k ∈ [0, 3*4[之间的映射。由于步幅等于 (-n, n) = (-1, 1),映射为 -n*i + n*j = -1*i + 1*j = j-i。从数学上讲,您可以将矩阵描述为 M[i, j] = F[j-i],其中 F 是扁平值向量 [3 2 1 4 5 6].

例如,让我们尝试使用 i=1j=2。如果你看上面的 Topleitz 矩阵 M[1, 2] = 4。确实F[k] = F[j-i] = F[1] = 4

如果你仔细观察,你会发现负步幅背后的诀窍:它们允许你 'reference' 负指数:例如,如果你采用 j=0i=2,那么你看k=-2。记住 vals 是如何通过切片 vals[len(c)-1:] 得到偏移量 2 的。如果你查看它自己的底层数据存储,它仍然是 [3 2 1 4 5 6],但有一个偏移量。由于偏移量,vals(在本例中为 i: 1D -> k: 1D)的映射将是 M'[i] = F'[k] = F'[i+2]。这意味着 M'[-2] = F'[0] = 3.

在上面我将M'定义为vals[len(c)-1:]这基本上等同于以下张量:

>>> torch.as_strided(vals, size=(len(vals)-2,), stride=(1,), storage_offset=2)
tensor([1, 4, 5, 6])

同样,我将 F' 定义为基础值的扁平向量:[3 2 1 4 5 6].

strides的使用确实是定义Toeplitz矩阵的一种非常巧妙的方式!


> 2. 是否有一种简单的方法可以重写此代码以便负步幅有效?

问题是,PyTorch 中没有实现负步长...我不相信 torch.as_strided 可以解决这个问题,否则扩展当前的实现并提供支持该功能。

不过,还有其他方法可以解决该问题。完全可以在 PyTorch 中构建托普利茨矩阵,但 torch.as_strided.

则不行

我们自己做映射:对于(i, j)索引的M的每个元素,我们会找出对应的索引k,也就是j-i。这可以轻松完成,首先从 M:

收集所有 (i, j)
>>> i, j = torch.ones(3, 4).nonzero().T
(tensor([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2]),
 tensor([0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3]))

现在我们基本上有 k:

>>> j-i
tensor([ 0,  1,  2,  3, -1,  0,  1,  2, -2, -1,  0,  1])

我们只需要从行 r 和列 c 输入中构建所有可能值的扁平张量。负索引值(c 的内容)放在最后并翻转:

>>> values = torch.cat((r, c[1:].flip(0)))
tensor([1, 4, 5, 6, 3, 2])

最后用 k 索引 values 并重塑:

>>> values[j-i].reshape(3, 4)
tensor([[1, 4, 5, 6],
        [2, 1, 4, 5],
        [3, 2, 1, 4]])

总而言之,我建议的实施方式是:

def toeplitz(c, r):
    vals = torch.cat((r, c[1:].flip(0)))
    shape = len(c), len(r)
    i, j = torch.ones(*shape).nonzero().T
    return vals[j-i].reshape(*shape)

> 3. 我可以通过 w.r.t c(或 r)到 toeplitz_torch 传递梯度吗?

这是个有趣的问题,因为 torch.as_strided 没有实现后向函数。这意味着您将无法反向传播到 cr!然而,使用上述方法,它使用 'backward-compatible' 内置函数,反向传递是免费的。

注意输出中的 grad_fn

>>> toeplitz(torch.tensor([1.,2.,3.], requires_grad=True), 
             torch.tensor([1.,4.,5.,6.], requires_grad=True))
tensor([[1., 4., 5., 6.],
        [2., 1., 4., 5.],
        [3., 2., 1., 4.]], grad_fn=<ViewBackward>)

这是一个草稿(写下来确实花了一些时间),我会做一些修改。
如果您有任何问题或意见,请随时发表评论!
我对看到其他答案很感兴趣,因为我不是大步向前的专家,这只是我对这个问题的看法。