<p>这本质上是一个<a href="https://en.wikipedia.org/wiki/Block_matrix" rel="nofollow noreferrer">matrix block</a>运算,可以从简单的矩阵运算中概括出来。您的示例很好,因为您这样定义了块矩阵(我将它们重命名为<code>A_p</code>和<code>u_p</code>,以避免与矩阵<code>u</code>冲突):</p>
<pre><code>A_p = np.array([[Xx, Xy],
[Yx, Yy]])
u_p = np.array([[u], [v]])
</code></pre>
<p>如果你看一个标准的矩阵运算,就形状而言,你有<code>(r, s)</code>遇到<code>(s, t)</code>,这最终导致了一个形状为<code>(r, t)</code>的矩阵。“<code>s</code>”维度在该过程中减少。如果您还不熟悉<a href="https://stackoverflow.com/questions/26089893/understanding-numpys-einsum">Einstein Summation</a>和<a href="https://numpy.org/doc/stable/reference/generated/numpy.einsum.html" rel="nofollow noreferrer">^{<cd8>}</a>,那么应该快速查看一下。给定矩阵<code>M</code>和<code>N</code>,矩阵操作如下所示:</p>
<pre><code>out = np.zeros((r, t))
for i in range(r):
for j in range(s):
for k in range(t):
out[i, k] += M[i, j]*N[j, k]
</code></pre>
<p>当你意识到操作是如何进行的时,这变得非常直观:覆盖整个多维空间,我们在输出矩阵上累积产品</p>
<p>对于<code>np.einsum</code>,这将非常类似(这里是<code>M=A</code>,这里是<code>N=u</code>):</p>
<pre><code>>>> np.einsum('ij,jk->ik', A, u)
array([[3],
[9]])
</code></pre>
<p>使用坐标:<code>i</code>、<code>j</code>和<code>k</code>而不是维度大小来思考是很有用的:<code>r</code>、<code>s</code>和<code>t</code>,就像上面使用<code>np.einsum</code>时一样。所以我会在剩下的答案中保留前一个符号</p>
<hr/>
<p>现在来看块矩阵,输入矩阵本身包含矩阵块。矩阵<code>M</code>是形状<code>(i, j, (m, n))</code>(<em>即</em><code>(i, j, m, n)</code>),其中<code>(k, l)</code>表示块的形状。同样,对于<code>N</code>,我们有<code>(j, k, m, n)</code>。我们基本上用2x2块{<cd28>}替换了标量形状的{<cd27>},仅此而已</p>
<p>因此,您可以直接推断<code>A_p*u_p</code>操作为:</p>
<pre><code>>>> np.einsum('ijmn,jkmn->ikmn', A_p, u_p)
array([[[[ 0, 50],
[200, 450]]],
[[[ 0, 110],
[440, 990]]]])
</code></pre>
<p>您可以将其与“手工制作”结果进行比较:</p>
<pre><code>>>> np.stack((Xx*u + Xy*v, Yx*u + Yy*v))[:, None]
</code></pre>
<p>注意,这可以用于任何<code>(m, n)</code>块和任何矩阵大小</p>
<p>您现在可以尝试使用类似于<code>A*u</code>的循环来实现<code>A_p*u_p</code></p>
<hr/>
<p><em>编辑:</em>在<a href="https://stackoverflow.com/users/15752941/loopy-walt">@loopy walt</a>的<a href="https://stackoverflow.com/a/68825058/6331369">answer</a>上展开,他向您展示了可以通过以下方式获得相同的结果:</p>
<pre><code>>>> (u_p.T @ A_p.T).T
</code></pre>
<p>实际上<a href="https://numpy.org/doc/stable/reference/generated/numpy.matmul.html" rel="nofollow noreferrer">^{<cd33>}</a>可以处理多维数组,并将以以下方式处理(对于高于2的维度数):两个输入都将<em>视为驻留在最后两个索引</em>中的矩阵堆栈。输入形状为<code>(i, j, m, n)</code>和<code>(j, k, m, n)</code>。为了使用<code>__matmul__</code>,我们需要将块维度(即<code>(m, n)</code>)拉到第一个位置,并将矩阵维度(分别为<code>(i, j)</code>和<code>(j, k)</code>)放在最后。需要采取以下步骤</p>
<ul>
<li><p>转置将产生颠倒维度顺序的效果:<code>(n, m, j, i)</code>和<code>(n, m, k, j)</code></p>
</li>
<li><p>然后我们交换两个操作数,我们设法以<code>(*, k, j) x (*, j, i) = (*, k, i)</code>的形式计算一些东西,其中<code>* = (n, m)</code></p>
</li>
<li><p>转换输出会导致<code>(i, k, m, n)</code>的形状</p>
</li>
</ul>