1

概要

mldivide を使用して、3x3 行列と逆 3x3 行列の間で行列乗算を繰り返し実行して、コードの時間効率を向上させようとしています。

バックグラウンド

歩行分析で被験者の下肢に取り付けられたセンサー間 の手と目のキャリブレーションに方向データを使用する前に、ベクトル量子化方法を実装しようとしています...私が従っているアルゴリズムは論文からのものです。必要..."Data Selection for Hand-eye Calibration: A Vector Quantization Approach"

最適化するコード

考えられるすべての「相対移動」(A または B) を解決するより高速な方法を見つけたいと思っていましたが、これには時間がかかりすぎます (C と D は約 2000 要素の長さであるため、A または B のサイズは になります =2000*(2000-1)/2=1999000)。

%C,D is cell array, each cell is 3x3 rotation matrix.
Nframe=length(C);
Nrel = Nframe*(Nframe-1)/2;
A=cell(Nrel,1);
B=cell(Nrel,1);
At=zeros(Nrel,4);
Bt=zeros(Nrel,4);
count = 1;
for i=1:Nframe
    for j=1:Nframe
        if j <= i
            continue;
        else
            % Main place to optimize (avoid looping?)!!
            % In DCM representation (ie. each cell is 3x3 matrix)
            A{count,1} = C{j}\C{i}; %=inv(C{i+x})*C{i}
            B{count,1} = D{j}\D{i};

            %Using the following adds to ~ 4000% to the time (according to tic toc).
            %Represent as axis-angle (ie. each cell -> 1x4 vec) - perhaps compute later
            At(count,:) = SpinConv('DCMtoEV',A{count}); %on Matlab File Exchange
            Bt(count,:) = SpinConv('DCMtoEV',B{count});
            count=count+1;
        end
    end
end

これが質問するのに適切な場所であることを願っています.適用できる以前の解決策を見つけることができませんでした. また、私には実際の経験がないため、大きな行列を扱う場合に計算時間が避けられないかどうかはわかりません。

-------------------------------------------------- ----------------------------------------------

*編集*

行列のプロパティ: 以下にコメントされているように、回転型です。それらは特殊直交群 SO(3) [transpose=inverse] にあります。http://en.wikipedia.org/wiki/Rotation_matrix#Properties_of_a_rotation_matrixを参照してください。

テスト方法: ランダムな回転行列 R を作成するには、次のコードを使用します。

[U,S,V] = svd(randn(3,3)); 
R = U∗V';
if det(R) < 0 
    S(1 ,1) = 1; 
    S(2 ,2) = 1;
    S(3,3) = −1; 
    R = U∗S∗V’;
end

SpinConv : 3x3 方向余弦行列から軸角度表現に変換するために使用しています。それはより複雑で、安定性のために必要以上に変換します (最初にクォータニオンに)。リンクは次のとおりです: http://www.mathworks.com/matlabcentral/fileexchange/41562-spinconv/content/SpinConv.m 行う必要があるのは次のとおりです (SpinConv ではなく、メソッドをすばやく実装しただけです)。

t = (trace(R)-1)/2;
% this is only in case of numerical errors
if t < -1,
    t = -1;
elseif t>1,
    t = 1;
end
theta = acosd(t);
if isequalf(180,theta) || isequalf(0,theta),
    axis = null(R-eye(3));
    axis = axis(:,1)/norm(axis(:,1));
else
    axis = [R(3,2) - R(2,3); R(1,3) - R(3,1); R(2,1) - R(1,2) ];
    axis = axis/(2*sind(theta));
end
At(count,:) = [-axis,theta]; %% NOTE (-ve) as noted in comments of correct answer.

*編集#2 * 代わりに、クォータニオンを使用して3x3マトリックスの使用を避けることができることに気づきました:

したがって、クォータニオンは 1x4 ベクトルです。元のコードは次のように変更できます (else ステートメント内):

A(count,:) = qnorm(qmult(qconj(C(j,:)),C(i,:)));
vec = [q(1) q(2) q(3)]/norm([q(1) q(2) q(3)]);
theta = 2*atan2(norm([q(1) q(2) q(3)]),q(4));
At(count,:)=[vec,theta];

ここで、qconjqmult、およびqnormは四元数演算です。


わかりました、申し訳ありません-それが私が持っているすべての情報と可能性です。

4

2 に答える 2

3

上でコメントしたように、最速の方法は行列のプロパティに大きく依存します。たとえば、一部のアルゴリズムは行列が対称であることから大きな恩恵を受けることができますが、そうでない場合はかなり遅くなります。

したがって、詳細な情報がなければ、いくつかの一般的なステートメントを作成し、ランダム行列でいくつかの方法を比較することしかできませ(通常、逆行列のコンテキストでは適切な比較ができません)。

MATLAB のバージョン (R2011a の JIT または大幅に改善されたもの) によっては、事前に割り当てAB、ループのパフォーマンスを大幅に向上させます。ループ内で配列を動的に拡張することは、通常、非常に非効率的です。

SpinConvこれは外部関数 (MEX か m かは関係ありません) であるため、JIT はこのループをコンパイルできないため、インタープリターの速度によって制限されます。これはかなり低いです。可能であれば、関連する部分をSpinConvループ本体にコピーして貼り付けるだけで、これを回避できます。非常に煩わしいことはわかっています (MATLAB の将来のバージョンでこれが自動化されることを願っています) が、今のところ、JIT にループ構造を理解させてコンパイルさせる唯一の方法です (実際、100 倍以上の係数は珍しくありません)。 )。

そうは言っても、私は2つの異なる方法をテストしました:

  1. Cとの LU 分解を計算して保存Dし、ループで再利用する
  2. Cn \ [C{1} C{2} ... C{n-1}]すべての について連立方程式を解き、n = 2:N並べ替えます。

コード内:

clc
clear all

Nframe = 500;

%// Some random data
C = cellfun(@(~)rand(3), cell(Nframe,1), 'UniformOutput', false); 
D = cellfun(@(~)rand(3), cell(Nframe,1), 'UniformOutput', false);


%// Your original method 
tic

count  = 1;
for i=1:Nframe
    for j=1:Nframe
        if j <= i
            continue;
        else
            %// Main place to optimize (avoid looping?)!!
            %// In DCM representation (ie. each cell is 3x3 matrix)
            A{count,1} = C{j}\C{i}; %=inv(C{i+x})*C{i}
            B{count,1} = D{j}\D{i};

            count=count+1;
        end
    end
end

toc 
A1 = A;



%// First method: compute all LU decompositions and re-use them in the loop
%// ------------------------------------------------------------------------

tic

%// Compute LU decompositions of all C and D
Clu = cell(Nframe, 2);
Dlu = cell(Nframe, 2);
for ii = 1:Nframe
    [Clu{ii,1:2}]  = lu(C{ii});
    [Dlu{ii,1:2}]  = lu(D{ii});
end

%// improvement: pre-allocate A and B
A = cell(Nframe*(Nframe-1)/2, 1);
B = cell(Nframe*(Nframe-1)/2, 1);

%// improvement: don't use i and j as variable names
count  = 1;
for ii = 1:Nframe

    %// improvement: instead of continue if j<=i, just use different range
    for jj = ii+1 : Nframe

        %// mldivide for LU is equal to backwards substitution, which is
        %// trivial and thus fast
        A{count} = Clu{jj,2}\(Clu{jj,1}\C{ii});
        B{count} = Dlu{jj,2}\(Dlu{jj,1}\D{ii});

        count = count+1;

    end
end

toc
A2 = A;



%// Second method: solve all systems simultaneously by concatenation
%// ------------------------------------------------------------------------

tic

% Pre-allocate temporary matrices
Aa = cell(Nframe-1, 1);
Bb = cell(Nframe-1, 1);

for ii = 2:Nframe   
    % Solve Cn \ [C1 C2 C3 ... Cn]
    Aa{ii-1} = C{ii}\[C{1:ii-1}];
    Bb{ii-1} = D{ii}\[D{1:ii-1}];
end
toc

%// Compared to the order in which data is stored in one of the other
%// methods, the order of data in Aa and Bb is different. So, we have to
%// re-order to get the proper order back: 

tic

A = cell(Nframe*(Nframe-1)/2, 1);
B = cell(Nframe*(Nframe-1)/2, 1);
for ii = 1:Nframe-1

     A( (1:Nframe-ii) + (Nframe-1-(ii-2)/2)*(ii-1) ) = ...
         cellfun(@(x) x(:, (1:3) + 3*(ii-1)), Aa(ii:end), 'UniformOutput', false);

     B( (1:Nframe-ii) + (Nframe-1-(ii-2)/2)*(ii-1) ) = ...
         cellfun(@(x) x(:, (1:3) + 3*(ii-1)), Bb(ii:end), 'UniformOutput', false);
end

toc
A3 = A;

% Check validity of outputs
allEqual = all( cellfun(@(x,y,z)isequal(x,y)&&isequal(x,z), A1,A2,A3) )

結果:

Elapsed time is 44.867630 seconds.  %// your original method 
Elapsed time is 1.267333 seconds.   %// with LU decomposition
Elapsed time is 0.183950 seconds.   %// solving en-masse by concatenation
Elapsed time is 1.871149 seconds.   %// re-ordering the output of that

allEqual = 
    1

私はR2010aを使用しているため、元のメソッドの遅さは主に事前に割り当てられていないことがA原因であることに注意してください。Bこの点に関して、新しい MATLAB バージョンのパフォーマンスは向上しますが、事前に割り当てた方がパフォーマンスは向上することに注意してください。

直感的に (そしておそらく他の人が示唆するように)、明示的な逆数を計算できます。

Cinv = cellfun(@inv, C, 'UniformOutput', false);

あるいは

Cinv = cellfun(@(x) [...
    x(5)*x(9)-x(8)*x(6)  x(7)*x(6)-x(4)*x(9)  x(4)*x(8)-x(7)*x(5) 
    x(8)*x(3)-x(2)*x(9)  x(1)*x(9)-x(7)*x(3)  x(7)*x(2)-x(1)*x(8)
    x(2)*x(6)-x(5)*x(3)  x(4)*x(3)-x(1)*x(6)  x(1)*x(5)-x(4)*x(2)] / ...
        (x(1)*x(5)*x(9) + x(4)*x(8)*x(3) + x(7)*x(2)*x(6) - ...
         x(7)*x(5)*x(3) - x(4)*x(2)*x(9) - x(1)*x(8)*x(6)),...
    C, 'UniformOutput', false);

(これはより高速で正確になります)、ループ内で単純に乗算します。おわかりのように、これは en-masse ソルブCn\[C1 C2 ... Cn-1]と LU の両方よりも大幅に遅くなります (ただし行列の性質によって異なります)。また、生成に失敗し allEqual == trueます。違いが小さい場合もありますが、多くの場合 (特に、特異に近い行列やその他の特殊な行列の場合)、違いは非常に大きくなります。

ここのSOに関する他の多くの質問で述べたように、また洗練されたGoogle検索や高度な線形代数の本が教えてくれるように、数値アプリケーションで明示的な逆数を使用すると、通常は遅くなり、常に不正確になり、時には危険ですらあります. 逆は非常に優れた理論的構成要素ですが、その理論を実際に適用する場合にはほとんど役に立ちません。したがって、上記の他の方法のいずれかを使用することをお勧めします。

結論は:

  • 順不同でデータを扱うことができる場合 (後でより複雑なインデックス作成が必要になる可能性があります)、連結によってシステムをまとめて解決するのがはるかに高速です。確かに、データを並べ替える方法は改善できますが、並べ替える必要がある場合は、LU の方が常に高速になると思います。

  • そうでなくても、行列が LU 分解に適している場合は、それを使用してください。これが当てはまるかどうかを確認するには、実際のデータとプロファイルで使用してください。また、LU の追加の出力を試すこともできます (最も注目すべきは、置換行列P、またはスパース行列の場合は、列の並べ替え行列Q)。

  • もちろん、QR 分解の方が適している場合は、 を使用しますqrchol、またはなどについても同じですpcg。さまざまな方法を少し試してみてください。

編集

あなたが言及したように、すべての行列は SO​​(3) 回転行列です。うわー、それは非常に重要な情報です !その場合、逆関数は単なる転置であり、逆関数のどのバリアントよりも 1 桁または 2 桁高速です。また、これらの回転行列を軸角度表現に変換することを示します。次に、カーネルを次の行に沿って何かに変更する必要があります

A = C{ii}.'*C{jj};
B = D{ii}.'*D{jj};

[V,lam] = eig(A);
axis  = V(:,imag(diag(lam))==0);
theta = acos(min(max(-1, (trace(A)-1)/2), 1));
At{count, :} = {axis theta*180/pi};

[V,lam] = eig(B);
axis  = V(:,imag(diag(lam))==0);
theta = acos(min(max(-1, (trace(B)-1)/2), 1));
Bt{count, :} = {axis theta*180/pi};

これは組み込み関数のみを使用するため、かなり効率的です。多くの非組み込み関数 ( 、、、 ) を使用するSpinConvため、少なくともコピー貼り付けよりは優れています。上記の方法は固有値法を使用していることに注意してください。関数で使用される行列式メソッドを使用すると、関数を更新して非組み込み関数を呼び出さないようにすれば、もう少し効率的にすることができます。SpinConvnullisequalfacosdsindSpinConv

注意: そのバージョンのSpinConvの軸の符号が間違っているようです。で計算された軸の符号は、 で計算された軸の符号elseifと反対ですif

于 2013-07-16T07:48:03.320 に答える
1

行列分割に複数回表示されるinv(C{j})ため、計算と保存を試みます。C{j}同上D{j}。それとも、あなたの 3x3 行列は特異ですか?

于 2013-07-15T20:31:44.640 に答える