1

概括

希望提高我的代码的时间效率,在 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];

其中qconjqmultqnorm是四元数运算。


好吧,很抱歉 - 这就是我所拥有的所有信息和可能性。

4

2 回答 2

3

正如我上面评论的那样,最快的方法很大程度上取决于矩阵的属性。例如,某些算法可以极大地受益于对称矩阵,但如果不是,则相当慢。

所以没有进一步的信息,我只能做一些一般性的陈述,并在随机矩阵上比较一些方法(这通常不会矩阵逆的情况下给出很好的比较)。

根据您的 MATLAB 版本(R2011a 中的 JIT 对其进行了极大的改进),预先分配AB为您提供循环性能的巨大增益;在循环内动态增长数组通常效率很低。

同样是对SpinConv: 的调用,因为这是一个外部函数(MEX 或 m,没关系),JIT 无法编译此循环,因此您受到解释器速度的限制。这是相当低的。如果可能的话,您可以通过简单地将相关部分复制粘贴SpinConv到循环体中来避免这种情况。我知道,这非常烦人(我当然希望这在未来版本的 MATLAB 中是自动化的),但目前它是让 JIT 理解循环结构并编译它的唯一方法(真的,100 或更多的因素并不少见)。

所以,说了这么多,我测试了两种不同的方法:

  1. 计算和存储 和 的 LU 分解,CD在循环中重新使用它们
  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 上,因此原始方法的缓慢主要是由于AB不是预先分配。请注意,在这方面,较新的 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);

(这将更快,更准确),然后在循环内简单地相乘。正如您将看到的,这比整体求解Cn\[C1 C2 ... Cn-1]和 LU 都慢得多(尽管取决于矩阵的性质)。此外,它无法生产 allEqual == true;有时差异很小,但通常(特别是对于近奇异矩阵和其他特殊矩阵),差异很大

正如这里关于 SO 的许多其他问题中所提到的,以及任何精炼的谷歌搜索或高级线性代数书籍都会告诉你的那样,在数值应用中使用显式逆通常会很慢,总是不准确,有时甚至是危险的。逆是一个非常好的理论结构,但在该理论的任何实际应用中几乎没有用。因此,最好使用上述其他方法之一。

综上所述:

  • 如果您可以忍受乱序的数据(以后可能需要更复杂的索引),那么通过串联解决系统整体问题是迄今为止最快的。当然,我重新排序数据的方式可以改进,但我怀疑如果你需要重新排序,LU 仍然会更快。

  • 如果不是这种情况,但您的矩阵适合 LU 分解,请使用它。要确定是否是这种情况,只需在您的真实数据和个人资料中使用它即可。您还可以尝试 LU 的其他输出(最值得注意的是置换矩阵P,或者对于稀疏矩阵,列重新排序矩阵Q)。

  • 当然,如果 QR 分解更合适,使用qr. chol, or等​​也一样pcg。尝试不同的方法。

编辑

正如您所提到的,所有矩阵都是 SO(3) 旋转矩阵。哇,这是重要的信息吗!在这种情况下,逆只是转置,它比逆的任何变体都快一个或两个数量级。此外,您表示要将这些旋转矩阵转换为轴角表示。然后应该将内核更改为类似于以下内容的内容

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,因为SpinConv它使用了很多非内置函数(null, isequalf, acosd, sind)。注意上面的方法使用的是特征值法;如果您使用SpinConv函数中使用的行列式方法,则可以稍微提高效率,前提是您将其更新为不调用非内置函数。

请注意:该版本SpinConv的轴符号似乎不正确;在 中计算的轴的符号与在 中计算的轴的符号elseif相反if

于 2013-07-16T07:48:03.320 回答
1

我会尝试计算和存储inv(C{j}),因为C{j}多次出现在矩阵中。同上D{j}。还是您的 3x3 矩阵是单数的?

于 2013-07-15T20:31:44.640 回答