如何向量化双线性和二次形式的评估？答案

【问题标题】：How to vectorize the evaluation of bilinear & quadratic forms?如何向量化双线性和二次形式的评估？
【发布时间】：2012-01-17 10:21:09
【问题描述】：

给定任意n x n实系数矩阵A，我们可以定义一个双线性形式bA : Rⁿ x Rⁿ → R by

b_A(x, y) = x ^T是的,

和二次形式q_A：Rⁿ → R by

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

(对于二次形式q_A的最常见应用，矩阵A是对称的，甚至是对称正定的，所以随意如果这对您的答案很重要，则假设其中任何一种情况都是如此。）

（另外，FWIW、b_I 和 q_I（其中 I 是 n x n 单位矩阵）分别是标准内积和平方 L² -norm on 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 将计算出所需的 b_A(x_,i, y_,i) 和 q_A(x_,i) (如结果 m x m 矩阵的对角线元素），以及 O(m²) 不相关 b_A(x_,i, y_,j) 和 b_A(x_,i, x_,j), (对于 i ≠ j), 所以这是一个非首发。

【问题讨论】：

标签： python r matlab matrix numpy

【解决方案1】：

如果你想在 MATLAB 中做这件事，这真的很简单：

你可以输入

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

我怀疑这是否值得，但如果你想加快速度，你可以试试这个：

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

【讨论】：

你没有回答这个问题。如果 x 和 y 是列向量矩阵，而 F 是方阵，那么正确答案是值向量，而不是结果中的矩阵。我认为正确的论坛是sum((x'*F).*y', 2)。
或者更简单：sum(x.*(F*y))
@peci1 我不明白为什么结果应该是向量而不是矩阵。以下将使我相信矩阵应该是结果：“[X]、[Y] 和 [A]，尺寸分别为 mn xm、mn xm 和 mn x mn，以及（块）对角线元素 x， i、y、i 和 A。那么所需的计算将是矩阵乘法 [X]T[A][Y] 和 [X]T[A][X]"
只需阅读问题的第一段。 b_A : R^n x R^n → R 是单个输入向量的双线性形式。结果只是一个实数。所以如果你把m输入向量，结果应该只是m实数。

【解决方案2】：

这是 numpy 中的一个解决方案，应该可以满足您的需求：

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

这对 X^T * A 进行矩阵乘法，然后进行逐元素数组乘法以乘以 Y^T。然后将结果数组的行求和以产生一维数组。

【讨论】：

谢谢！您是说 ((np.matrix(X).T*np.matrix(A)).A * np.matrix(Y).T.A).sum(1) 吗？
如果将 X、A 和 Y 存储为矩阵，它可以更短：((X.T*A).A *Y.T.A).sum(1)。我包括 np.matrix 调用主要是为了阐明矩阵乘法与逐元素乘法。
或者只是坚持使用数组：(np.dot(X.T, A) * Y.T).sum(1)

【解决方案3】：

尚不完全清楚您要实现的目标，但在 R 中，您使用 crossprod 形成叉积：给定矩阵 X 和 Y 具有兼容的尺寸，crossprod(X, Y) 返回 X^TY。同样，矩阵乘法是通过%*% 运算符实现的：X %*% Y 返回乘积 XY。因此，您可以将 X^TAY 设为 crossprod(X, A %*% Y)，而无需担心矩阵乘法、循环等机制。

如果您的矩阵具有允许优化计算的特定结构（对称、三角形、稀疏、带状等），您可以查看Matrix 包，它对此提供了一些支持。

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

【讨论】：