0

我正在尝试获取大型数组中每个向量之间的角度（1896378x4-EDIT：这意味着我需要1.7981e+12角度......太大，但如果有优化代码的空间，无论如何请告诉我）。它太慢了——我还没有看到它完成。以下是我采取的优化步骤：

首先，从逻辑上讲，我（认为我）想要什么（仅Bt=rand(N,4)用于测试）：

    [ro,col]=size(Bt);
    angbtwn = zeros(ro-1); %too long to compute!! total non-zero = ro*(ro-1)/2
    count=1;
    for ii=1:ro-1
        for jj=ii+1:ro
            angbtwn(count) = atan2(norm(cross(Bt(ii,1:3),Bt(jj,1:3))), dot(Bt(ii,1:3),Bt(jj,1:3))).*180/pi;
            count=count+1;
        end
    end

所以，虽然我会尝试对它进行矢量化，并摆脱非内置函数：

[ro,col]=size(Bt);
% angbtwn = zeros(ro-1); %TOO LONG!
for ii=1:ro-1
    allAxes=Bt(ii:ro,1:3);
    repeachAxis = allAxes(ones(ro-ii+1,1),1:3);
    c = [repeachAxis(:,2).*allAxes(:,3)-repeachAxis(:,3).*allAxes(:,2)
        repeachAxis(:,3).*allAxes(:,1)-repeachAxis(:,1).*allAxes(:,3)
        repeachAxis(:,1).*allAxes(:,2)-repeachAxis(:,2).*allAxes(:,1)];
    crossedAxis = reshape(c,size(repeachAxis));
    normedAxis = sqrt(sum(crossedAxis.^2,2));
    dottedAxis = sum(repeachAxis.*allAxes,2);
    angbtwn(1:ro-ii+1,ii) = atan2(normedAxis,dottedAxis)*180/pi;
end
angbtwn(1,:)=[]; %angle btwn vec and itself
%only upper left triangle are values...

仍然太长，甚至预分配......所以我尝试做稀疏，但没有正确实施：

[ro,col]=size(Bt);
%spalloc:
angbtwn = sparse([],[],[],ro,ro,ro*(ro-1)/2);%zeros(ro-1); %cell(ro,1)
for ii=1:ro-1
    ...same
    angbtwn(1:ro-ii+1,ii) = atan2(normedAxis,dottedAxis)*180/pi; %WARNED: indexing = >overhead
    % WHAT? Can't index sparse?? what's the point of spalloc then?
end

因此，如果我的逻辑可以改进，或者如果 sparse 确实是要走的路，而我就是无法正确实施，请告诉我在哪里改进。谢谢你的帮助。

4

2 回答 2

0

您可能希望使用pdistwith 'cosine'distance 来计算1-cos(angbtwn).
这种方法的另一个好处是它不计算n^2值，而是完全.5*(n-1)*n唯一的值:)

于 2013-08-06T15:11:23.360 回答
0

你想得到每对向量之间的角度Bt吗？如果Bt有 200 万个向量，每个向量（显然）需要一个内积来获得之间的角度。我不知道任何类型的优化都会帮助在单台机器上的 MATLAB 中在合理的时间内完成此操作。

在任何情况下，您都可以将此问题转化为单位向量矩阵之间的矩阵乘法：

N=1000;
Bt=rand(N,4); % for testing. A matrix of N (row) vectors of length 4.
[ro,col]=size(Bt);

magnitude = zeros(N,1); % the magnitude of each row vector.
units = zeros(size(Bt)); % the unit vectors

% Compute the unit vectors for the row vectors
for ii=1:ro
  magnitude(ii) = norm(Bt(ii,:));
  units(ii,:) = Bt(ii,:) / magnitude(ii);
end

angbtwn = acos(units * units') * 360 / (2*pi);

但是在对较大的 N 进行矩阵乘法时，您将耗尽内存。

于 2013-08-03T04:18:42.623 回答