如何向量化双线性和二次形式的计算？

Question

给定一个任意的 n x n 的实数矩阵 A，我们可以定义一个双线性形式 b_A : Rⁿ x Rⁿ → R，它的定义是：

b_A(x, y) = x^TAy ,

还有一个二次形式 q_A : Rⁿ → R，它的定义是：

q_A(x) = b_A(x, x) = x^TAx .

（在大多数常见的二次形式应用中，矩阵 A 是对称的，或者甚至是对称正定的，所以如果这对你的答案有影响，可以假设其中之一是成立的。）

（另外，顺便提一下，b_I 和 q_I（其中 I 是 n x n 的单位矩阵）分别是标准内积和平方的 L² 范数在 Rⁿ 上，也就是 x^Ty 和 x^Tx。）

现在假设我有两个 n x m 的矩阵，X 和 Y，还有一个 n x n 的矩阵 A。我想优化计算 b_A(x_,i, y_,i) 和 q_A(x_,i)（这里 x_,i 和 y_,i 表示 X 和 Y 的第 i 列），我猜想在一些环境中，比如 numpy、R 或 Matlab，这将涉及某种形式的向量化处理。

我能想到的唯一解决方案是生成对角块矩阵 [X]、[Y] 和 [A]，它们的维度分别是 mn x m、mn x m 和 mn x mn，并且（块）对角元素分别是 x_,i、y_,i 和 A。然后，想要的计算就是矩阵乘法 [X]^T[A][Y] 和 [X]^T[A][X]。这个策略确实不太灵光，但如果有更高效的时间和空间利用方式，我很想看看。（不言而喻，任何不利用这些块矩阵稀疏性的实现都是注定失败的。）

有没有更好的方法呢？

我希望使用的系统是 numpy，但如果有其他支持高效矩阵计算的系统的答案，比如 R 或 Matlab，也可以（前提是我能搞明白如何将它们移植到 numpy）。

谢谢！

_{当然，计算 X^TAY 和 X^TAX 会得到想要的 bA(x,i, y,i) 和 qA(x,i)（作为结果 m x m 矩阵的对角元素），以及 O(m²) 的无关的 bA(x,i, y,j) 和 bA(x,i, x,j)（对于 i ≠ j），所以这不是个好主意。}

Answer 1

如果你想在MATLAB里做这件事，其实非常简单：

你只需要输入

b = x'*A*y;
q = x'*A*x;

我不太确定这样做是否值得，但如果你想加快速度，可以试试这个：

M = x'*A;
b = M*y;
q = M*x;

Answer 2

你想要实现的目标不是特别清楚，但在R语言中，你可以用 crossprod 来计算矩阵的交叉乘积。假设你有两个矩阵 X 和 Y，它们的尺寸是可以匹配的，使用 crossprod(X, Y) 就可以得到 X^TY 的结果。类似地，矩阵相乘可以用 %*% 这个符号来实现：X %*% Y 会返回 XY 的结果。因此，你可以通过 crossprod(X, A %*% Y) 来得到 X^TAY，而不需要担心矩阵相乘的具体过程、循环或者其他复杂的东西。

如果你的矩阵有特定的结构，比如对称、三角形、稀疏或者带状等，可以考虑使用 Matrix 这个包，它对这些情况有一些支持。

我没有用过Matlab，但我相信它也会有类似的功能来完成这些操作。

Answer 3

这里有一个用numpy实现的解决方案，可以满足你的需求：

((np.matrix(X).T*np.matrix(A)).A * Y.T.A).sum(1)

这个代码首先计算了X的转置（X^T）和A的矩阵乘法，然后进行逐元素的数组相乘，接着再乘以Y的转置（Y^T）。最后，得到的数组的每一行会被加起来，最终形成一个一维数组。