matlab代码优化-在一个大单元阵列中将所有3x3矩阵相乘

Question

摘要

为了提高我的代码的时间效率，重复地在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

希望这是正确的地方问，我没有找到一个以前的解决方案，我可以申请。而且，我没有真正的经验，所以我不确定在处理大型矩阵时，计算时间是否是不可避免的。

*EDIT*

矩阵属性：旋转，如下所示-它们是‘好’，而不是单数。它们属于特殊正交群，因此(3) transpose=inverse。请参阅矩阵

方法来测试：要创建随机旋转矩阵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矢量。可以将原始代码更改为(inside 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和Q范数是四元数操作。

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

Answer 1

如前所述，最快的方法强烈地依赖于矩阵的属性。例如，一些算法可以从矩阵的对称性中受益很大，但如果不对称，则会比较慢。

因此，如果没有进一步的信息，我只能做一些一般性的陈述，并比较一些关于随机矩阵的方法(在矩阵逆的情况下，这通常不能给出一个很好的比较)。

根据您的MATLAB版本( R2011a中的JIT或相当大的改进)，预分配A和B在循环性能上会给您带来很大的提高；在循环中动态地增长数组通常效率很低。

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

所以，说了这些之后，我测试了两种不同的方法：

1. 计算和存储C和D的LU分解，并在循环中重用它们。
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);

(这将更快、更准确)，然后简单地在循环中进行乘法。正如您将看到的那样，这要比Cn\[C1 C2 ... Cn-1]和LU都慢得多(尽管这取决于矩阵的性质)。而且，它不能产生allEqual == true；有时差异很小，但是通常(特别是对于近奇异矩阵和其他特殊情况)，差异是巨大的。

正如很多其他问题中所提到的，任何经过改进的Google搜索或高级线性代数书都会告诉你，在数值应用中使用显式逆通常是缓慢的，总是不准确的，有时甚至是危险的。逆是一个很好的理论构造，但在任何实际应用中都是无用的。因此，最好使用上述其他方法之一。

最后：

• 如果您可以忍受数据无序(这可能需要以后进行更复杂的索引)，那么通过连接来解决整个系统是最快的。当然，我重新排序数据的方式是可以改进的，但我怀疑如果您需要重新排序，LU仍然会更快。
• 如果不是这样，但是您的矩阵适合LU分解，请使用它。要确定是否是这样，只需将其用于您的真实数据和配置文件。您还可以尝试LU的附加输出(最显著的是置换矩阵P，或者对于稀疏矩阵，列重排序矩阵Q)。
• 当然，如果QR分解更合适，则使用qr。对于chol，或pcg等也是一样，用不同的方法做一些实验。

编辑

正如您提到的，所有的矩阵都是如此(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中计算的轴的符号相反。

Answer 2

我会尝试计算和存储inv(C{j})，因为C{j}多次出现在矩阵除法中。D{j}也是。还是你的3x3矩阵是单数的？