用对数来考虑它怎么样：

import numpy

a = numpy.random.rand(10, 3)
b = numpy.random.rand(3, 5)

out = np.exp(np.matmul(np.log(a), b))

因为c_ij = prod(a_ik ** b_kj, k=1..K)，所以log(c_ij) = sum(log(a_ik) * b_ik, k=1..K)。
注意：在a中有零 * 可能 * 会搞乱结果（也是负数，但结果无论如何都不会被很好地定义）。我不知道NumPy是否保证了这种行为，但为了安全起见，你可以在最后添加一些东西，比如：

out[np.logical_or.reduce(a < eps, axis=1)] = 0

您可以在将这些数组扩展到3D版本后使用broadcasting-

(a[:,:,None]**b[None,:,:]).prod(axis=1)

简单地说-

(a[...,None]**b[None]).prod(1)

基本上，我们保持两个数组的最后一个轴和第一个轴对齐，同时在两个输入的第一个轴和最后一个轴之间执行元素幂。

10   x   3   x   1
   1   x   3   x   5

还有两种解决方案：
内联

np.array([
    np.prod([a[:, i]**bb[i] for i in range(len(bb))], axis=0)
    for bb in b.T
]).T

使用power.outer：

np.prod([numpy.power.outer(a[:, k], b[k]) for k in range(len(b))], axis=0

两者都比广播解决方案慢一点。

用于重现绘图的代码：

import numpy as np
import perfplot

def loop(a, b):
    m = a.shape[0]
    n = b.shape[1]
    out = np.ones((m, n))
    for i in range(m):
        for j in range(n):
            for k in range(3):
                out[i, j] *= a[i, k] ** b[k, j]
    return out

def broadcasting(a, b):
    return (a[..., None] ** b[None]).prod(1)

def log_exp(a, b):
    neg_a = np.zeros(a.shape, dtype=int)
    neg_a[a < 0.0] = 1
    odd_b = np.zeros(b.shape, dtype=int)
    odd_b[b % 2 == 1] = 1
    negative_count = np.dot(neg_a, odd_b)

    out = (-1) ** negative_count * np.exp(
        np.matmul(np.log(abs(a), where=abs(a) > 0.0), b)
    )

    zero_a = np.zeros(a.shape, dtype=int)
    zero_a[a == 0.0] = 1
    pos_b = np.zeros(b.shape, dtype=int)
    pos_b[b > 0] = 1
    zero_count = np.dot(zero_a, pos_b)
    out[zero_count > 0] = 0.0
    return out

def inline(a, b):
    return np.array(
        [np.prod([a[:, i] ** bb[i] for i in range(len(bb))], axis=0) for bb in b.T]
    ).T

def outer_power(a, b):
    return np.prod([np.power.outer(a[:, k], b[k]) for k in range(len(b))], axis=0)

b = perfplot.bench(
    setup=lambda n: (
        np.random.rand(n, 3) - 0.5,
        np.random.randint(0, 10, (3, n)),
    ),
    n_range=[2**k for k in range(13)],
    kernels=[loop, broadcasting, inline, log_exp, outer_power],
    xlabel="len(a)",
)
b.save("out.png")
b.show()

import numpy

a = numpy.random.rand(10, 3)
b = numpy.random.rand(3, 5)

out = [[numpy.prod([a[i, k]**b[k, j] for k in range(3)]) for j in range(5)] for i in range(10)]