二等分方法的网格化应用 [英] Gridwise application of the bisection method

问题描述

我需要找到广义状态空间的根.也就是说，我有一个尺寸为grid=AxBx(...)xX的离散网格，我事先不知道它有多少尺寸(该解决方案应适用于任何grid.size).

I need to find roots for a generalized state space. That is, I have a discrete grid of dimensions grid=AxBx(...)xX, of which I do not know ex ante how many dimensions it has (the solution should be applicable to any grid.size) .

我想使用

I want to find the roots (f(z) = 0) for every state z inside grid using the bisection method. Say remainder contains f(z), and I know f'(z) < 0. Then I need to

• 如果remainder> 0
，则增加z
• 如果remainder<减小z 0

• increase z if remainder > 0
• decrease z if remainder < 0

Wlog，例如形状为(grid.shape, T)的矩阵history包含网格中每个点的z早期值的历史记录，我需要增加z(因为remainder> 0).然后，我需要在history[z, :]中选择zAlternative，这是其中最小的，大于z".用伪代码，即:

Wlog, say the matrix historyof shape (grid.shape, T) contains the history of earlier values of z for every point in the grid and I need to increase z (since remainder > 0). I will then need to select zAlternative inside history[z, :] that is the "smallest of those, that are larger than z". In pseudo-code, that is:

zAlternative =  hist[z,:][hist[z,:] > z].min()

我之前已经问过了.给出的解决方案是

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

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

newZ = history[..., -1]
criterion = 0.05
increase = remainder > 0 + criterion
decrease = remainder < 0 - criterion
newZ[increase] = 0.5*(newZ[increase] + higherZ[increase])
newZ[decrease] = 0.5*(newZ[decrease] + lowerZ[decrease])

但是，此代码对我而言不再起作用.我对接受它感到非常难过，但是我从未理解索引正在发生的魔力，因此很遗憾，我需要帮助.

However, this code ceases to work for me. I feel extremely bad about admitting it, but I never understood the magic that is happening with the indices, therefore I unfortunately need help.

代码实际上是做什么的，它给了我最低的和最高的.也就是说，如果我修复两个特定的z值:

What the code actually does, it to give me the lowest respectively the highest. That is, if I fix on two specific z values:

history[z1] = array([0.3, 0.2, 0.1])
history[z2] = array([0.1, 0.2, 0.3])

我将得到higherZ[z1] = 0.3和lowerZ[z2] = 0.1，即极值.两种情况的正确值都是0.2.这是怎么了

I will get higherZ[z1] = 0.3 and lowerZ[z2] = 0.1, that is, the extrema. The correct value for both cases would have been 0.2. What's going wrong here?

如果需要，为了生成测试数据，可以使用

If needed, in order to generate testing data, you can use something along the lines of

history = tile(array([0.1, 0.3, 0.2, 0.15, 0.13])[newaxis,newaxis,:], (10, 20, 1))
remainder = -1*ones((10, 20))

测试第二种情况.

预期结果

我在上面调整了history变量，以给出向上和向下的测试用例.预期结果将是

I adjusted the history variable above, to give test cases for both upwards and downwards. Expected outcome would be

lowerZ = 0.1 * ones((10,20))
higherZ = 0.15 * ones((10,20))

对于历史记录[z，:]中的每个点z，哪个是第二高的前一个值(higherZ)和第二个最小的前一个值(lowerZ).由于所有点z具有完全相同的历史记录([0.1, 0.3, 0.2, 0.15, 0.13])，因此对于lowerZ和higherZ它们都将具有相同的值.当然，通常，每个z的历史记录都将不同，因此，这两个矩阵在每个网格点上都可能包含不同的值.

Which is, for every point z in history[z, :], the next highest previous value (higherZ) and the next smallest previous value (lowerZ). Since all points z have exactly the same history ([0.1, 0.3, 0.2, 0.15, 0.13]), they will all have the same values for lowerZ and higherZ. Of course, in general, the histories for each z will be different and hence the two matrices will contain potentially different values on every grid point.

推荐答案

我将您在此处发布的内容与上一篇文章的解决方案进行了比较并注意到一些差异.

I compared what you posted here to the solution for your previous post and noticed some differences.

对于较小的 z，您说

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

他们说:

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

对于更大 z，您说

mask = b <= history[..., -1:]
index = argmax(mask, axis=-1)

他们说:

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

使用您上一篇文章的解决方案，我得到:

import numpy as np
history = np.tile(np.array([0.1, 0.3, 0.2, 0.15, 0.13])[np.newaxis,np.newaxis,:], (10, 20, 1))
remainder = -1*np.ones((10, 20))

a = history

# b is a sorted ndarray excluding the most recent observation
# it is sorted along the observation axis
b = np.sort(a[..., :-1], axis=-1)

# mask is a boolean array, comparing the (sorted)
# previous observations to the current observation - [..., -1:]
mask = b > a[..., -1:]

# The next 5 statements build an indexing array.
# True evaluates to one and False evaluates to zero.
# argmax() will return the index of the first True,
# in this case along the last (observations) axis.
# index is an array with the shape of z (2-d for this test data).
# It represents the index of the next greater
# observation for every 'element' of z.
index = np.argmax(mask, axis=-1)

# The next two statements construct arrays of indices
# for every element of z - the first n-1 dimensions of history.
indices = tuple([np.arange(j) for j in b.shape[:-1]])
indices = np.meshgrid(*indices, indexing='ij', sparse=True)

# Adding index to the end of indices (the last dimension of history)
# produces a 'group' of indices that will 'select' a single observation
# for every 'element' of z
indices.append(index)
indices = tuple(indices)
higherZ = b[indices]

mask = b >= a[..., -1:]
# Since b excludes the current observation, we want the
# index just before the next highest observation for lowerZ,
# hence the minus one.
index = np.argmax(mask, axis=-1) - 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)
lowerZ = b[indices]
assert np.all(lowerZ == .1)
assert np.all(higherZ == .15)

似乎可行

这篇关于二等分方法的网格化应用的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！