是否有计算二次形式（x'*A*x）的内置matlab？

Question

非常简单的问题：给定一个 N x N 对称矩阵 A 和一个 N 向量 x，是否有内置的 Matlab 函数来计算 x'*A*x？即，不是y = x'*A*x，而是有一个函数quadraticform s.t。 y = quadraticform(A, x)?

显然我可以做到y = x'*A*x，但我需要性能，而且似乎应该有办法利用

1. A 是对称的
2. 左右乘法器是同一个向量

如果没有一个内置函数，有没有比x'*A*x 更快的方法？或者，Matlab 解析器是否足够智能以优化x'*A*x？如果是这样，你能指出我在文档中验证事实的地方吗？

Answer 1

我找不到这样的内置函数，我知道为什么。

y=x'*A*x 可以写成n^2 项A(i,j)*x(i)*x(j) 的总和，其中i 和j 从1 运行到n（其中A 是nxn 矩阵） . A 是对称的：A(i,j) = A(j,i) 对于所有 i 和 j。由于对称性，除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来实现，所以特殊的“二次型”函数并没有真正的优势。

=== 更新 ===

我用 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)，但前者的因数较小：

Answer 2

MATLAB 足够聪明，可以识别和优化某些类型的复合矩阵表达式，我相信（尽管我无法确定）二次形式是它所做的优化之一。

但是，这不是 MathWorks 倾向于记录的那种东西，因为 a) 它通常只会在函数内进行优化，而不是在脚本中、命令行或调试中进行优化 b) 它可能只在某些情况下有效，例如对于真正的非稀疏 A c）它可能会随着版本的变化而变化，所以他们不希望你依赖它 d) 它是使 MATLAB 如此出色的专有技术之一。

要确认，您可以尝试比较 y=x'*A*x 和 B=A*x; y=x'*B 的时间。您也可以尝试feature('accel','off')，这将关闭大部分此类优化。

最后，如果您联系 MathWorks 支持，您也许可以让其中一位开发人员确认是否正在进行优化。