概要
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];
ここで、qconj、qmult、およびqnormは四元数演算です。
わかりました、申し訳ありません-それが私が持っているすべての情報と可能性です。