log2元素的两个numpy矩阵之间的矩阵乘法方法

Question

我有两个numpy 矩阵a_log2 和b_log2 log2 元素，我想在它们之间执行矩阵乘法。

a = np.array([[0.4, 0.4, 0.2],
          [0.1, 0.5, 0.4]])

b = np.array([[0.3, 0.7],
          [0.5, 0.5],
          [0.2, 0.8]])

a_log2 = np.log2(a)
b_log2 = np.log2(b)

我曾经使用scipy.special.logsumexp 执行e-based logarithms 的矩阵乘法。这是我使用的代码（感谢 Erik Parkinson 在此线程 Handling matrix multiplication in log space in Python 中的回答）：

def log_space_product(A,B):
   Astack = np.stack([A]*A.shape[0]).transpose((2,1,0))
   Bstack = np.stack([B]*B.shape[1]).transpose((1,0,2))
   log_sum_exp = logsumexp(Astack+Bstack, axis=0)
   return log_sum_exp

现在，我需要您的帮助来执行a_log2 和b_log2 之间的矩阵乘法，因为scipy.special.logsumexp 没有为base-2 logarithms 定义。

注意： 我最初打算使用a_loge[i, j] = np.log(2**a_log2[i, j])和b_loge[i, j] = np.log(2**b_log2[i, j])将矩阵元素转换为natural logarithms，然后使用前面提到的log_space_product()方法进行矩阵乘法。

但我拒绝自己这样做，因为我最终将使用的矩阵有> 1000 行和> 20 列。 （不要对这里的行数和列数感到困惑。我确实确保矩阵乘法属性保持不变。）

Answer 1

记住 log₂(x) = ln(x)/ln(2)。我们可以用 c = ln(2) 将其重写为 ln(x) = c log₂(x)。因此，与具有不同基数的日志空间之间的转换实际上是一个乘法。

from scipy.special import logsumexp

def log_space_product(A, B, base=np.e):
    c = np.log(base)
    Astack = np.stack([A]*A.shape[0]).transpose(2,1,0)
    Bstack = np.stack([B]*B.shape[1]).transpose(1,0,2)
    return (1/c) * logsumexp(c*(Astack+Bstack), axis=0)