3

非常に簡単な質問: N x N 対称行列 A と N ベクトル x が与えられた場合、計算する組み込みの Matlab 関数はありますx'*A*xか? つまり、 の代わりにy = x'*A*x、関数quadraticformstがありy = quadraticform(A, x)ますか?

明らかに私はそれを行うことができますy = x'*A*xが、私はパフォーマンスが必要であり、それを利用する方法があるべきだと思われます

  1. A対称です
  2. 左と右の乗数は同じベクトルです

単一の組み込み関数がない場合、より高速な方法はありx'*A*xますか? または、Matlab パーサーは最適化するのに十分スマートx'*A*xですか? もしそうなら、その事実を証明しているドキュメントの場所を教えていただけますか?

4

2 に答える 2

6

そのような組み込み関数が見つかりませんでした。その理由はわかりました。

y=x'*A*xn^2は項の和として書くことができA(i,j)*x(i)*x(j)、ここでiとjは から1まで実行されますn(ここでAはnxn行列です)。Aは対称です:A(i,j) = A(j,i)すべてのiとj. i対称性により、等しい項を除いて、すべての項が合計で 2 回表示されjます。したがって、n*(n+1)/2用語が異なります。それぞれに 2 つの浮動小数点乗算があるため、単純なメソッドn*(n+1)では合計で乗算が必要になります。x'*A*xの単純な計算、つまり を計算z=A*xしてから を計算するのにも乗算y=x'*zが必要であることは容易にわかります。ただし、さまざまな項を合計するより高速なn*(n+1)方法があります。n*(n+1)/2ix(i)、つまりn*(n-1)/2+3*n掛け算だけで十分です。しかし、これはあまり役に立ちません: の計算の実行時間y=x'*A*xはまだO(n^2)です。

したがって、二次形式の計算は よりも速く行うことはできないと思います。O(n^2)また、これも式 によって達成できるためy=x'*A*x、特別な「二次形式」関数の実際の利点はありません。

===更新===

Matlab の拡張機能として、関数「quadraticform」を C で作成しました。

// 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' A xと比較する簡単なパフォーマンス テストを作成しました。

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)ますが、前者の係数は小さくなります。

quadraticformtest.png

于 2011-12-03T20:05:31.690 に答える
1

MATLAB は、ある種の複合行列式を認識して最適化するのに十分なほど賢く、2 次形式はそれが行う最適化の 1 つであると私は信じています (確実に確認することはできませんが)。

ただし、これは MathWorks が文書化する傾向があるものではありません。なぜなら、a) 通常、スクリプト、コマンド ライン、またはデバッグではなく、関数内でのみ最適化される b) 実際のアプリケーションなど、特定の状況でのみ機能する可能性があるためです。非スパース A c) リリースごとに変更される可能性があるため、信頼してほしくない d) MATLAB を非常に優れたものにしている独自仕様の 1 つです。

y=x'*A*x確認するには、 のタイミングをと比較してみてくださいB=A*x; y=x'*B。を試すこともできますfeature('accel','off')。これにより、これらの種類の最適化のほとんどがオフになります。

最後に、MathWorks サポートに連絡すると、開発者の 1 人に最適化が行われているかどうかを確認してもらうことができる場合があります。

于 2011-12-03T21:02:44.477 に答える