Python/Scipy：查找矩阵的“有界”最小值/最大值

Question

我认为最容易说明我的问题，笼统的情况很难解释。

假设我有一个矩阵

a with dimensions NxMxT,

可以将 T 视为时间维度（使问题更容易）。令 (n,m) 为通过 NxM 的索引。我可以称 (n,m) 为状态空间标识符。然后我需要找到对应的 python/scipy

for each (n,m):
     find a*(n,m) = min(a(n,m,:) s.t. a*(n,m) > a(n,m,T)

也就是说，对于整个状态空间，找到仍然高于最后（时间维度中）观察值的最小状态空间值。

我的第一次尝试是先解决内部问题（找到高于a[...,-1]的a）：

aHigherThanLast = a[ a > a[...,-1][...,newaxis] ]

然后我想为每个 (n,m) 找到所有这些中最小的。不幸的是，aHigherThanLast 现在包含所有这些值的一维数组，所以我不再有 (n,m) 对应关系。有什么更好的方法来解决这个问题？

另外一个问题：状态空间是可变的，它也可能是 3 维或更多维（NxMxKx...），我无法对此进行硬编码。所以任何一种

for (n,m,t) in nditer(a):

不可行。

非常感谢！

/编辑：

a = array([[[[[[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]],



          [[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]]],




         [[[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]],



          [[[[ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.],
             [ 0.,  2.,  1.]]]]]]]])
# a.shape = (1L, 1L, 2L, 2L, 1L, 1L, 10L, 3L). so in this case, T = 3.
# expected output would be the sort of
# b.shape = (1L, 1L, 2L, 2L, 1L, 1L, 10L), which solves

• b[a,b,c,d,e,f,g] > a[a,b,c,d,e,f,g,-1]（b高于最新观测值）

  • a中的i中没有元素同时满足

    -- a[a,b,c,d,e,f,g,t] > a[a,b,c,d,e,f,g,-1]

    -- a[a,b,c,d,e,f,g,t]

所以，鉴于前一个数组是一个简单的堆栈，如果 [0,2,1] 沿着最后一次观察，我希望

b = ones((1,1,2,2,1,1,10))*2

然而， - 如果在一些 (a,b,c,d,e,f,g) 中，不仅有 {0,1,2} 的值，还有 {3}，那么我仍然想要 2 （因为它是满足 i > 1 的 i = {2,3} 中的较小者。 - 如果在一些 (a,b,c,d,e,f,g) 中只有值 {0,1,3}，我想要 3，因为 i = 3 将是满足 i 的最小数> 1.

希望这能澄清一点吗？

/edit2：

非常感谢答案，它有效。如果我想要相反的情况，即较小的那些中最大的，我将如何调整它？我没有尝试通过那个复杂的索引逻辑，所以我只改变前三行的（弱）尝试没有成功：

        b = sort(a[...,:-1], axis=-1)
        b = b[...,::-1]
        mask = b < a[..., -1:]
        index = argmax(mask, axis=-1)
        indices = tuple([arange(j) for j in a.shape[:-1]])
        indices = meshgrid(*indices, indexing='ij', sparse=True)
        indices.append(index)
        indices = tuple(indices)
        a[indices]

另外，我的第二次尝试 a[...,::-1][indices] 也没有结果。

Answer 1

我认为 E 先生是在正确的轨道上。您肯定会首先对没有最后一次值的数组进行排序：

b = np.sort(a[..., :-1], axis=-1)

您现在最好使用`np.searchsorted 来查找比最终值大的第一项在哪里，但不幸的是np.searchsorted 仅适用于扁平数组，因此我们必须做更多工作，例如创建布尔掩码，然后使用np.argmax找到第一个True：

mask = b > a[..., -1:]
index = np.argmax(mask, axis=-1)

您现在有了索引，要提取实际值，您需要做一些索引魔术：

indices = tuple([np.arange(j) for j in b.shape[:-1]])
indices = np.meshgrid(*indices, indexing='ij', sparse=True)
indices.append(index)
indices = tuple(indices)

你现在终于可以做到了：

>>> b[indices]
array([[[[[[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]],


          [[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]]],



         [[[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]],


          [[[ 2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.]]]]]]])
>>> b[indices].shape
(1L, 1L, 2L, 2L, 1L, 1L, 10L)

---

要在较小的那些中获得最大的，您可以执行以下操作：

mask = b >= a[..., -1:]
index = np.argmax(mask, axis=-1) - 1

即较小的项目中最大的是，在相同或更大的项目中，最小的项目之前的项目。第二种情况更清楚地表明，如果没有满足条件的项目，这种方法会给出垃圾结果。在第二种情况下，当这种情况发生时，您将获得一个索引的 -1，因此您可以通过 np.any(index == -1) 检查结果是否有效。

如果第一种情况不能满足条件，则可以将索引设置为-1

mask = b > a[..., -1:]
wrong = np.all(~mask, axis=-1)
index = np.argmax(mask, axis=-1)
index[wrong] = -1