如何向量化矩阵的行对角化答案

【问题标题】：How to vectorize row-wise diagonalization of a matrix如何向量化矩阵的行对角化
【发布时间】：2012-06-17 14:07:57
【问题描述】：

我有一个 n×m 矩阵，我想将它转换为一个 mn×m 矩阵，结果的每个 m×m 块都包含每一行的对角线。

例如，如果输入是：

[1 2; 3 4; 5 6]

输出应该是：

[1 0; 0 2; 3 0; 0 4; 5 0; 0 6]

当然，我不想自己用for循环一步步组装矩阵。
有没有一种矢量化且简单的方法来实现这一点？

【问题讨论】：

标签： matlab matrix octave vectorization

【解决方案1】：

对于执行此操作的矢量化方法，将对角元素的线性索引创建到结果矩阵中，然后直接分配。

%# create some input data
inArray = [10 11;12 13;14 15];

%# make the index array
[nr,nc]=size(inArray);

idxArray = reshape(1:nr*nc,nc,nr)';
idxArray = bsxfun(@plus,idxArray,0:nr*nc:nr*nc^2-1);

%# create output
out = zeros(nr*nc,nc);
out(idxArray) = inArray(:);

out =

    10     0
     0    11
    12     0
     0    13
    14     0
     0    15

【讨论】：

你能解释一下bsxfun的作用吗？
@ja72：bsxfun 自动扩展数组，以便您可以例如将一个 n×1 数组添加到一个 n×m 数组。开头的句柄（在这种情况下为@plus）告诉执行哪个操作。使用repmat更快更方便。

【解决方案2】：

这是一个简单的矢量化解决方案，假设X 是输入矩阵：

Y = repmat(eye(size(X, 2)), size(X, 1), 1);
Y(find(Y)) = X;

另一种选择是使用sparse，这可以写成一个简洁的单行：

Y = full(sparse(1:numel(X), repmat(1:size(X, 2), 1, size(X, 1)), X'));

【讨论】：

+1：这一定是被接受的解决方案； 2d step 是 Octave 矩阵索引的另一个惊人应用
此解决方案比使用 bsxfun 更可取，因为它不需要函数传递，因此可以在本机代码中更好地优化。另外，考虑使 Y 稀疏。
@dspyz sparse 是个好建议，谢谢。我添加了另一个解决方案:-)

【解决方案3】：

我看到的最简单的方法实际上非常简单，使用简单的索引引用和 reshape 函数：

I = [1 2; 3 4; 5 6];
J(:,[1,4]) = I;
K = reshape(J',2,6)';

如果你检查J，它看起来像这样：

J =
     1     0     0     2
     3     0     0     4
     5     0     0     6

Matrix K 正是我们想要的：

正如 Eitan T 在 cmets 中指出的那样，以上内容是针对示例的，并未涵盖一般解决方案。所以下面是一般的解决方案，其中 m 和 n 如问题中所述。

J(:,1:(m+1):m^2) = I;
K=reshape(J',m,m*n)';

如果您想对其进行测试以使其正常工作，只需使用

I=reshape(1:(m*n),m,n)';

注意：如果 J 已经存在，这可能会导致问题。在这种情况下，您还需要使用

J=zeros(n,m^2);

【讨论】：

-1：这会触发错误：“下标分配维度不匹配。”。即使它会起作用，那也不是正确的输出。
尺寸不匹配是由于拼写错误 - 我说 J(:,1:4)，但它应该说 J(:,[1,4])。我抄错了，仅此而已。错误已得到纠正，它应该可以按预期工作。至于输出，如果你注意的话，我描述的是矩阵 J，它不是最终的输出。 K 是输出，它看起来和 litro 要求的完全一样……至少，它是在 Octave 上的。
它现在可以工作了，所以我删除了反对票。但是请注意，它不是通用的，因为它仅在 I 有两列时才有效。
没错，它不是通用的，就像我写的那样。我将添加一个更通用的解决方案，因为它确实可以很好地概括。
完成。相当有信心它应该按预期工作。不要以为有什么错别字。

【解决方案4】：

这可能不是计算效率最高的解决方案，但这里有一个使用 kron 的 1-liner：

A = [1 2; 3 4; 5 6];
B = diag(reshape(A', 6, 1) * kron(ones(3, 1), eye(2))
% B = 
%     1     0
%     0     2
%     3     0
%     0     4
%     5     0
%     0     6

如果 A 为 n x m，则可以推广：

diag(reshape(A.', n*m, 1)) * kron(ones(n,1), eye(m))

【讨论】：