Есть ли встроенный Matlab, который вычисляет квадратичную форму (x'*A*x)?

Я не смог найти такую ​​встроенную функцию, и я знаю, почему.

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, не было бы никакой реальной пользы от специальной функции "квадратичной формы".

=== ОБНОВЛЕНИЕ ===

Я написал функцию «quadraticform» на C как расширение 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), но у первого коэффициент меньше:

MATLAB достаточно умен, чтобы распознавать и оптимизировать некоторые виды составных матричных выражений, и я полагаю (хотя не могу точно подтвердить), что квадратичная форма является одной из оптимизаций, которые он делает.

Однако это не то, что MathWorks склонен документировать, потому что а) это обычно оптимизируется только внутри функций, а не в сценариях, в командной строке или при отладке б) это может работать только в некоторых обстоятельствах, например, для реального неразреженный A c) он может меняться от релиза к релизу, поэтому они не хотят, чтобы вы на него полагались d) это одна из проприетарных вещей, которые делают MATLAB таким хорошим.

Чтобы подтвердить, вы можете попробовать сравнить тайминги для y=x'*A*x и B=A*x; y=x'*B. Вы также можете попробовать feature('accel','off'), который отключит большинство подобных оптимизаций.

Наконец, если вы обратитесь в службу поддержки MathWorks, вы сможете попросить одного из разработчиков подтвердить, выполняется ли оптимизация.