如果我的问题是正确的，那么问题就出来了$a_{ijk}$和$b_{ijkl}$来形成

C(i,j,k,l) = A(i,j,B(i,j,l,k));

4个嵌套循环

如果没有冒号符号，这将需要 4 个嵌套循环变量i,j,k,l：

C(n1,n2,N,n3) = 0.0;          % equivalent to C = zeros(n1,n2,N,n3);
for i = 1:n1,
  for j = 1:n2,
    for k = 1:N,
      for l = 1:n3,
        C(i,j,k,l) = A(i,j,B(i,j,l,k));
      end,
    end,
  end,
end,

请注意数组的预分配C：在循环中具有可变大小的数组for通常是效率低下的主要原因。事实上matlab，在评估 for 循环方面非常好（与 python 相比），前提是 lhs 在迭代期间不会改变大小。

2个嵌套循环

现在让我们进行简单的矢量化，消除 和 上的k循环l。这里的主要问题k,l是交换了C，B所以我们必须对permute数组维度：

C(n1,n2,n3,N) = 0.0;
for i = 1:n1,
  for j = 1:n2,
    C(i,j,:) = A(i,j,B(i,j,:));
  end,
end,
C = permute(C, [1 2 4 3]);

（请记住，数组中的尾随冒号将所有剩余维度压缩为一个：size(C(i,j,:)) == [1 1 n3*N]；保留正确的步幅，因此C(i,j,:)可以作为 lhs）

没有循环

编辑：更正错误并插入最终解决方案。

现在是问题的难点：如何消除i和上的循环j。这个想法很简单：找到一个整数数组indx，使得C == A(indx): 实际上A(indx)是一个具有相同形状的数组，indx并且从中提取的元素被A视为单个列。

为了构造这个数组，我们必须记住matlab遵循 fortran 约定：

size(A) == [n1 n2 n3]
A(i,j,k) == A(i + (j-1)*n1 + (k-1)*n1*n2)

所以这里是如何用 2 个嵌套循环构造这个矩阵：

indx(n1,n2,n3,N) = 0.0;
for j=1:n2,
  for i=1:n1,
    indx(i,j,:) = i+(j-1)*n1+(B(i,j,:)-1)*n1*n2;
  end,
end,

或有两条指令

indx = (B-1)*n1*n2;
indx(:) = indx(:) + repmat(1:n1*n2, n3*N)'

无循环代码读取

indx = (B-1)*n1*n2;
indx(:) = indx(:) + repmat(1:n1*n2, 1, n3*N)';
indx = permute(indx, [1 2 4 3]);
C = A(indx);

结论

在我看来，最有效的计算方式C应该是2 个嵌套循环版本。在无循环版本中，我们在计算上付出了很大的努力，indx而不是让matlab他自己做指针运算，所以我不鼓励这样做，除非同一个数组B，因此同一个向量indx可以多次用于不同的A数组。

如果是这种情况，我仍然建议indx通过两个嵌套循环来计算，而不是无循环版本。（这里戈德里克先知是对的：几乎不可能理解背后的逻辑。）

让我强调一个重要的一点：如果数组C是预先分配的，两个嵌套循环在 中将非常快matlab，并且进入无循环版本将没有实质性优势，除非我们有一种快速计算indx向量的方法。（例如，通过使用整数算术。）

希望这个冗长的回复会有用。

编辑 了以前错误的无循环版本：

C(n1,n2,n3,N) = 0.0;
indx =  ... % to be defined later
C(:) = A(indx);
C = permute(C, [1 2 4 3]);

size(C) == [n1 n2 n3 N]
C(i,j,k,l) == C(i + (j-1)*n1 + (k-1)*n1*n2 + (l-1)*n1*n2*n3 )

C(i,j,k,l) = A(i,j,B(i,j,k,l))

C(i,j,k,l) = A(i + (j-1)*n1 + (B(i,j,k,l)-1)*n1*n2 );

C(i,j,:) = A(i + (j-1)*n1 + (B(i,j,:)-1)*n1*n2);

indx = [];
for j=1:n2,
  for i=1:n1,
    indx = [ indx, reshape( i+(j-1)*n1+(B(i,j,:)-1)*n1*n2, 1, n3*N ) ];
  end,
end,

我假设您想对操作进行矢量化以提高效率。虽然矢量化分配可能是可能的（我个人已经尝试了几件事并且到目前为止都失败了）我将不得不争辩说你最好不要在这个应用程序中使用矢量化 matlab。

我的论点的原因是您的代码的清晰度将消失。即使是单行分配

C(m,l,:,k)=A(m,l,B(m,l,k,:))

需要一点时间来回想一下，在一行中执行三个循环所需的复杂向量表示法将完全掩盖其他阅读代码的人的分配过程（这包括一个月后的你）。

如果您必须从矢量化中看到效率提升，您可能会通过使用简短的 fortran 或 C 函数获得更好的性能和更清晰的内容。上面的循环可以用最少的努力来实现，两种语言在从你的 matlab 脚本中调用时可能会在性能上击败矢量化 matlab。

我知道这并不能直接回答您的问题，但希望在没有人为您的循环找到可接受的矢量化的情况下有所帮助。

感谢您的所有回复。我实际上在周末解决了这个问题，但没有机会查看这个论坛（坦率地说，我对这个社区的人们有多么乐于助人感到非常惊讶。我实际上只是假设我的问题会被置之不理鉴于那天晚上没有给出答案。）无论如何，我使用了以下命令，我认为这与第一个答案相同，但我没有机会彻底回顾慷慨长的回应。再次感谢两位响应者的时间和见解。我会阅读您的建议并将它们合并到我的代码中。

A(repmat(reshape(1:n1*n2,n1,n2),[1,1,n3,N])+n1*n2*(B-1);