通过普通矩阵运算复制“np.einsum”结果答案

【问题标题】：Replicating `np.einsum` result via normal matrix operations通过普通矩阵运算复制“np.einsum”结果
【发布时间】：2020-05-23 11:57:32
【问题描述】：

我已经通过 Numpy 在 Python 中实现了 TCB Spline。代码的关键部分如下所示：

np.einsum('km,km,kl,lm->m',xdiffpow_knot, h_pow_knot[:,i], hermite_matrix, lag_knot[:,i])

其中 k 和 l 始终为 4（k 是 0 到 3 的幂，l 是 TCB 样条使用的 4 个控制点），m 是数组的长度 @987654328 @我想插值。

当时我是通过np.einsum 实现的，因为如果没有np.einsum，我无法找出必要的矩阵运算来完成它。似乎我在结果中留下了一个额外的 m（请注意，前两个 km 术语只是逐元素乘法）。

现在我在没有 einsum 的情况下在 Julia 中重新实现（因此我可以利用 ForwardDiff、ReverseDiff 等中的算法微分）。如何通过矩阵运算复制上述 einsum？

我尝试了什么？

只考虑所涉及的维度并使点积起作用，感觉好像我缺少一个m-element 向量。唯一有意义的m 元素向量是 1-s 向量，我相信它可以作为求和。但在我依赖它之前，我想验证这在理论上是正确的。

完整代码。好丑……

HERMITE_MATRIX = np.array([[ 2.,-2., 1., 1.],
                           [-3., 3.,-2.,-1.],
                           [ 0., 0., 1., 0.],
                           [ 1., 0., 0., 0.]])


def hermite(x_knot, y_knot, tension=0.0, continuity=0.0, bias=0.0, weight_fwd_knot=None, hermite_matrix=HERMITE_MATRIX):
    order = 3
    h_knot = np.diff(x_knot)
    h_pow_knot = (1.0 / h_knot) ** np.arange(order, -1, -1)[:,None]
    is_first_knot = np.isclose(x_knot, x_knot[0])
    is_last_knot = np.isclose(x_knot, x_knot[-1])
    y_next_knot = np.roll(y_knot,-1)
    y_prev_knot = np.roll(y_knot, 1)
    x_next_knot = np.roll(x_knot,-1)
    x_prev_knot = np.roll(x_knot, 1)
    if weight_fwd_knot is None:
        weight_fwd_knot = np.where(is_last_knot, 0.0, np.where(is_first_knot, 1.0, (x_next_knot - x_knot)/(x_next_knot - x_prev_knot)))
    weight_bak_knot = 1.0 - weight_fwd_knot
    dydxfwd_knot = np.where(is_last_knot, (y_knot - y_prev_knot)/(x_knot - x_prev_knot), (y_next_knot - y_knot)/(x_next_knot - x_knot))
    dydxbak_knot = np.where(is_first_knot, (y_next_knot - y_knot)/(x_next_knot - x_knot), (y_knot - y_prev_knot)/(x_knot - x_prev_knot))
    dy_in_knot = (1 - tension) * ((1 + continuity) * (1 - bias) * dydxfwd_knot * weight_fwd_knot + (1 - continuity) * (1 + bias) * dydxbak_knot * weight_bak_knot)
    dy_out_knot = (1 - tension) * ((1 - continuity) * (1 - bias) * dydxfwd_knot * weight_fwd_knot + (1 + continuity) * (1 + bias) * dydxbak_knot * weight_bak_knot)
    lag_knot = np.array([y_knot[:-1], y_next_knot[:-1], dy_out_knot[:-1] * h_knot, dy_in_knot[1:] * h_knot])
    def f(x):
        i = np.maximum(np.minimum(np.searchsorted(x_knot, x, side="right") - 1, x_knot.size - 2), 0)
        xdiffpow_knot = (x - x_knot[i]) ** np.arange(order, -1, -1)[:,None]
        return np.einsum('km,km,kl,lm->m',xdiffpow_knot, h_pow_knot[:,i], hermite_matrix, lag_knot[:,i])
    return f

【问题讨论】：

我认为OMEinsum.jl 确实提供了您正在寻找的东西。但是这个函数不是不可微分的，主要是因为你调用了maximum、minimum和searchsorted吗？
通常不需要区分x 数组，它们是影响或受不可微函数影响的数组。 jacobian 几乎总是类似于dyknotdy。到目前为止，我已经在 Julia 中实现了线性、三次和指数张力样条曲线，所有 dyknotdy 似乎都按预期工作。
对不起，把那个翻过来。 dydyknot 是雅可比人。

标签： julia numpy-einsum

【解决方案1】：

考虑：

>>> a = np.array([[1,2],[3,4]])
>>> np.einsum('km,km,kl,lm->m',a,a,a,a)
array([142, 392])

这可以使用基本的线性代数运算来计算。
观察 'kl,lm' 部分是传统的矩阵乘法，可以进行子计算以产生 'km'：

>>> x = np.matmul(a,a)

现在剩下的 'km,km,km->m' 是元素乘法并在索引 'k' 上求和

>>> y = a * a * x
>>> y
array([[  7,  40],
       [135, 352]])

>>> np.sum(y, axis=0)
array([142, 392])

【讨论】：