问题很简单：给定一个N x N对称矩阵A和一个N向量x，是否有内置的Matlab函数来计算x'*A*x？也就是说，是否有一个函数quadraticform使得y = quadraticform(A, x)，而不是y = x'*A*x？
显然，我可以直接做y = x'*A*x，但我需要更好的性能，看起来应该有一种方法利用以下两点：
A是对称的如果没有单个内置函数，是否有比x'*A*x更快的方法？或者，Matlab解析器是否足够聪明以优化x'*A*x？如果是这样，请指出文档中验证此事实的位置。
y=x'*A*x 可以写成 A(i,j)*x(i)*x(j) 的 n^2 项和，其中 i 和 j 运行从 1 到 n（其中 A 是一个 nxn 矩阵）。A 是对称的：对于所有的 i 和 j，有 A(i,j) = A(j,i)。由于对称性，在和中每个项都出现了两次，除了那些 i 等于 j 的项。因此我们有 n*(n+1)/2 种不同的项。每个项有两个浮点乘法，所以天真的方法需要总共 n*(n+1) 个乘法。很容易看出，天真的计算 x'*A*x，即计算 z=A*x，然后 y=x'*z，也需要 n*(n+1) 个乘法。然而，有一种更快的方法来求和我们的 n*(n+1)/2 不同的项：对于每个 i，我们可以分解出 x(i)，这意味着只需要 n*(n-1)/2+3*n 个乘法就足够了。但这并不能真正帮助：计算 y=x'*A*x 的运行时间仍然是 O(n^2)。O(n^2) 更快，由于这也可以通过公式 y=x'*A*x 实现，所以特殊的“quadraticform”函数没有真正的优势。
=== 更新 ===
我已经用 C 写了一个名为“quadraticform”的函数，作为 Matlab 扩展。// y = quadraticform(A, x)
#include "mex.h"
/* Input Arguments */
#define A_in prhs[0]
#define x_in prhs[1]
/* Output Arguments */
#define y_out plhs[0]
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize mA, nA, n, mx, nx;
double *A, *x;
double z, y;
int i, j, k;
if (nrhs != 2) {
mexErrMsgTxt("Two input arguments required.");
} else if (nlhs > 1) {
mexErrMsgTxt("Too many output arguments.");
}
mA = mxGetM(A_in);
nA = mxGetN(A_in);
if (mA != nA)
mexErrMsgTxt("The first input argument must be a quadratic matrix.");
n = mA;
mx = mxGetM(x_in);
nx = mxGetN(x_in);
if (mx != n || nx != 1)
mexErrMsgTxt("The second input argument must be a column vector of proper size.");
A = mxGetPr(A_in);
x = mxGetPr(x_in);
y = 0.0;
k = 0;
for (i = 0; i < n; ++i)
{
z = 0.0;
for (j = 0; j < i; ++j)
z += A[k + j] * x[j];
z *= x[i];
y += A[k + i] * x[i] * x[i] + z + z;
k += n;
}
y_out = mxCreateDoubleScalar(y);
}
我将这段代码保存为“quadraticform.c”，并使用Matlab进行编译：
mex -O quadraticform.c
我写了一个简单的性能测试，以比较这个函数与x'Ax的表现：
clear all; close all; clc;
sizes = int32(logspace(2, 3, 25));
nsizes = length(sizes);
etimes = zeros(nsizes, 2); % Matlab vs. C
nrepeats = 100;
h = waitbar(0, 'Please wait...');
for i = 1 : nrepeats
for j = 1 : nsizes
n = sizes(j);
A = randn(n);
A = (A + A') / 2;
x = randn(n, 1);
if randn > 0
start = tic;
y1 = x' * A * x;
etimes(j, 1) = etimes(j, 1) + toc(start);
start = tic;
y2 = quadraticform(A, x);
etimes(j, 2) = etimes(j, 2) + toc(start);
else
start = tic;
y2 = quadraticform(A, x);
etimes(j, 2) = etimes(j, 2) + toc(start);
start = tic;
y1 = x' * A * x;
etimes(j, 1) = etimes(j, 1) + toc(start);
end;
if abs((y1 - y2) / y2) > 1e-10
error('"x'' * A * x" is not equal to "quadraticform(A, x)"');
end;
waitbar(((i - 1) * nsizes + j) / (nrepeats * nsizes), h);
end;
end;
close(h);
clear A x y;
etimes = etimes / nrepeats;
n = double(sizes);
n2 = n .^ 2.0;
i = nsizes - 2 : nsizes;
n2_1 = mean(etimes(i, 1)) * n2 / mean(n2(i));
n2_2 = mean(etimes(i, 2)) * n2 / mean(n2(i));
figure;
loglog(n, etimes(:, 1), 'r.-', 'LineSmoothing', 'on');
hold on;
loglog(n, etimes(:, 2), 'g.-', 'LineSmoothing', 'on');
loglog(n, n2_1, 'k-', 'LineSmoothing', 'on');
loglog(n, n2_2, 'k-', 'LineSmoothing', 'on');
axis([n(1) n(end) 1e-4 1e-2]);
xlabel('Matrix size, n');
ylabel('Running time (a.u.)');
legend('x'' * A * x', 'quadraticform(A, x)', 'O(n^2)', 'Location', 'NorthWest');
W = 16 / 2.54; H = 12 / 2.54; dpi = 100;
set(gcf, 'PaperPosition', [0, 0, W, H]);
set(gcf, 'PaperSize', [W, H]);
print(gcf, sprintf('-r%d',dpi), '-dpng', 'quadraticformtest.png');
x'*A*x 还是 quadraticform(A,x) 的运行时间都收敛于 O(n^2)，但前者的因子更小:
n*(n-1)/2+3*n。 - kolMATLAB非常聪明，能够识别和优化某些类型的复合矩阵表达式，我相信（虽然我不能确定）二次形式是它所做的优化之一。
然而，这不是MathWorks通常会记录的内容，因为： a）它通常只在函数内部进行优化，而不是在脚本、命令行或调试中； b）它可能只适用于某些情况，例如对于实数非稀疏A； c）它可能会随着版本的发布而改变，因此他们不希望您依赖它； d）这是使MATLAB如此出色的专有技术之一。
要确认这一点，您可以尝试比较y=x'*A*x和B=A*x; y=x'*B的时间。您还可以尝试使用feature('accel','off')，这将关闭大多数这类优化。
最后，如果您联系MathWorks支持，您可能会得到其中一位开发人员的确认是否正在进行优化。