MATLAB：使用另一个矩阵作为列键/索引处理 LARGE 矩阵的特定列

Question

我有一个矩阵 A，它是 6 行 x 40 列。它填充有随机数。 我使用 nchoosek(1:40, 6) 创建了一个矩阵 B，其中包含列的所有可能 linear 组合的索引；一种可能的组合是 (1,2,3,4,5,6) 和 (1,2,3,4,5,40)。这个矩阵的大小是 3838380 x 6.

然而，矩阵 B 只给出了我想要制作的矩阵的索引，它也是 3838380 x 6，并且包含 A 的实际值，使用 B 的值进行索引。例如，B 的第 4 行是 [1 2 3 4 5 9].

对于矩阵 B 的所有 3838380 行，我需要使用 A 创建一个矩阵 C。这个新矩阵 C 使用 B 的一行作为索引来检索 A 的列，生成一个 6x6 矩阵：

[1     2     3     4     5     9] --> [0.4178    0.6562   -0.8633    0.7979   -0.9162    0.8720
    0.4864    0.0149   -0.8301   -0.2927   -0.7153   -0.2214
    0.7994   -0.2677   -0.8633   -0.7596   -0.8468   -0.7657
   -0.8695   -0.5467   -0.1804    0.1382    0.4811   -0.5192
   -0.3282    0.0697   -0.7532    0.7501   -0.0869    0.3698
   -0.9913   -0.4210   -0.1140   -0.3029    0.3365    0.6785].

C 用于求解 Cx = d 中的 x，其中 d 是一个 6x1 向量。

我正在使用 For 循环执行此操作 3838380 次：

for j = 1:length(B) %
    C = A(:,B(j,:))
    x = C\d % d is a 6x1 vector - x returns a 1x6 vector as well
end

目前大约需要 3 分钟。我想知道是否有更快或矢量化的方式，也许创建一个 6 x 3838380 x 6 的 3 维矩阵（每个条目都是一个 6x6 矩阵）？我仍然需要单独处理每个 6x6 矩阵并将返回的向量存储到另一个矩阵。

失败/我该怎么办？ 我实际上尝试创建 6 x 3838380 x 6 矩阵，完全使用 C = A(:,B)。它不适用于我的计算机，也不适用于我的大脑——我无法弄清楚分解这个巨大的矩阵与我之前所做的有什么不同。

我觉得有一种方法可以不用for循环来做到这一点

干杯！

Answer 1

这可以通过创建一个大小为 6×6×3838380 的 3D 数组来矢量化，其中包含所有 6×6 矩阵作为“页面”，然后使用

pagemldivide

（在 R2022a 中引入）一次求解所有线性系统:

% Example data
A = rand(6, 40);
d = rand(size(A,1), 1);
B = nchoosek(1:size(A,2), size(A,1));

% Non-vectorized approach
tic
x = NaN(size(A,1), size(B,1)); % initiallize
for jj = 1:length(B) %
    C = A(:, B(jj,:));
    x(:, jj) = C\d;
end
toc

% Vectorized approach
tic
CC = reshape(A(:, B.'), size(A,1), size(A,1), []);
xx = permute(pagemldivide(CC, d), [1 3 2]);
toc

% Check
isequal(x, xx)

矢量化版本似乎确实更快。我在 R2022b 中得到了这些结果（使用 Matlab Online）：

Elapsed time is 10.310424 seconds.
Elapsed time is 0.812846 seconds.
ans =
  logical
   1