概括
希望提高我的代码的时间效率,在 3x3 矩阵和逆 3x3 矩阵之间重复执行矩阵乘法 - 使用 mldivide。
背景
我正在尝试实现矢量量化方法,然后在步态分析中使用方向数据在连接到受试者下肢的传感器之间
进行手眼校准"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];
好吧,很抱歉 - 这就是我所拥有的所有信息和可能性。