如何根据每一行将矩阵转换为一堆对角矩阵?
How to convert matrix to a stack of diagonal matrices based on every row?
我有一个矩阵:
A = [1 1 1
2 2 2
3 3 3]
有没有矢量化的获取方式:
B = [1 0 0
0 1 0
0 0 1
2 0 0
0 2 0
0 0 2
3 0 0
0 3 0
0 0 3]
%// example data
data = reshape(1:9,3,3).' %'
n = 3; %// assumed to be known
data =
1 2 3
4 5 6
7 8 9
%// row indices
rows = 1:numel(data);
%// column indices
cols = mod(rows-1,n) + 1;
%// pre-allocation
out = zeros(n*n,n);
%// linear indices
linIdx = sub2ind(size(out),rows,cols);
%// assigning
out(linIdx) = data.'
out =
1 0 0
0 2 0
0 0 3
4 0 0
0 5 0
0 0 6
7 0 0
0 8 0
0 0 9
或者如果您更喜欢节省代码行而不是可读性:
out = zeros(n*n,n);
out(sub2ind(size(out),1:numel(data),mod((1:numel(data))-1,n) + 1)) = data.'
其他两个快速解决方案,但不比其他解决方案快:
%// #1
Z = blockproc(A,[1 size(A,2)],@(x) diag(x.data));
%// #2
n = size(A,2);
Z = zeros(n*n,n);
Z( repmat(logical(eye(n)),n,1) ) = A;
为了比赛 - Benchmark
function [t] = bench()
A = magic(200);
% functions to compare
fcns = {
@() thewaywewalk(A);
@() lhcgeneva(A);
@() rayryeng(A);
@() rlbond(A);
};
% timeit
t = zeros(4,1);
for ii = 1:10;
t = t + cellfun(@timeit, fcns);
end
format long
end
function Z = thewaywewalk(A)
n = size(A,2);
rows = 1:numel(A);
cols = mod(rows-1,n) + 1;
Z = zeros(n*n,n);
linIdx = sub2ind(size(Z),rows,cols);
Z(linIdx) = A.';
end
function Z = lhcgeneva(A)
sz = size(A);
Z = zeros(sz(1)*sz(2), sz(2));
for i = 1 : sz(1)
Z((i-1)*sz(2)+1:i*sz(2), :) = diag(A(i, :));
end
end
function Z = rayryeng(A)
A = A.';
Z = full(sparse(1:numel(A), repmat(1:size(A,2),1,size(A,1)), A(:)));
end
function Z = rlbond(A)
D = cellfun(@diag,mat2cell(A, ones(size(A,1), 1), size(A,2)), 'UniformOutput', false);
Z = vertcat(D{:});
end
ans =
0.322633905428601 %// thewaywewalk
0.550931853207228 %// lhcgeneva
0.254718792359946 %// rayryeng - Winner!
0.898236688657039 %// rlbond
此代码将 A 转换为行向量元胞数组,对每个行向量应用 diag
函数,然后将它们堆叠:
D = cellfun(@diag,mat2cell(A, ones(size(A,1), 1), size(A,2)), 'UniformOutput', false);
B = vertcat(D{:});
结果:
B =
1 0 0
0 1 0
0 0 1
2 0 0
0 2 0
0 0 2
3 0 0
0 3 0
0 0 3
编辑:我想 thewaywewalk 的基准测试只给我留下了可读性论点;)
使用 Beaker 的建议进行编辑:
data = [1 1 1
2 2 2
3 3 3];
sz = size(data);
z = zeros(sz(1)*sz(2), sz(2));
for i = 1 : sz(1)
z((i-1)*sz(2)+1:i*sz(2), :) = diag(data(i, :));
end
这是使用 sparse
and repmat
的另一种方法:
A = [1 2 3; 4 5 6; 7 8 9];
A = A.';
B = full(sparse(1:numel(A), repmat(1:size(A,1),1,size(A,2)), A(:)));
原始矩阵在 A
中,我将其转置以便下一步可以正确展开每个矩阵的行。我使用 sparse
来声明矩阵中的非零值。具体来说,我们看到每行只有一个条目,因此行索引的范围应该从 1 到 A
中的条目数。列从 1 到最后一列波动并重复。 mod
当然是通过 thewaywewalk 的解决方案的方式,但我想使用 repmat
以便这是一个独立于他的方法的解决方案。因此,我们创建了一个向量来访问从 1 到我们拥有的尽可能多的列,并且我们对我们拥有的尽可能多的行重复此操作。这些行和列索引向量将决定非零位置出现的位置。最后,将进入每个非零位置的是 A
的元素,这些元素按照行和列索引向量指定的顺序按行主要顺序展开。
请注意,在repmat
调用中,由于转置操作,调用size
时的行和列被颠倒了。
结果如下,我们得到:
>> B
B =
1 0 0
0 2 0
0 0 3
4 0 0
0 5 0
0 0 6
7 0 0
0 8 0
0 0 9
鉴于上述问题的稀疏性,将矩阵保留为 sparse
形式可能会更快,并且仅在必要时使用 full
进行转换。在两种格式之间进行转换需要时间,因此如果您决定进行基准测试,请考虑到这一点。
另一种方法(对我来说)比上面的任何基准测试方法(对于大矩阵和小矩阵)都快
[m,n] = size(A);
Z(n,m*n) = 0;
for idx = 1:n
Z(idx,((idx-1)*m+1):(idx*m)) = A(:,idx);
end
Z = reshape(Z,m*n,n);
我有一个矩阵:
A = [1 1 1
2 2 2
3 3 3]
有没有矢量化的获取方式:
B = [1 0 0
0 1 0
0 0 1
2 0 0
0 2 0
0 0 2
3 0 0
0 3 0
0 0 3]
%// example data
data = reshape(1:9,3,3).' %'
n = 3; %// assumed to be known
data =
1 2 3
4 5 6
7 8 9
%// row indices
rows = 1:numel(data);
%// column indices
cols = mod(rows-1,n) + 1;
%// pre-allocation
out = zeros(n*n,n);
%// linear indices
linIdx = sub2ind(size(out),rows,cols);
%// assigning
out(linIdx) = data.'
out =
1 0 0
0 2 0
0 0 3
4 0 0
0 5 0
0 0 6
7 0 0
0 8 0
0 0 9
或者如果您更喜欢节省代码行而不是可读性:
out = zeros(n*n,n);
out(sub2ind(size(out),1:numel(data),mod((1:numel(data))-1,n) + 1)) = data.'
其他两个快速解决方案,但不比其他解决方案快:
%// #1
Z = blockproc(A,[1 size(A,2)],@(x) diag(x.data));
%// #2
n = size(A,2);
Z = zeros(n*n,n);
Z( repmat(logical(eye(n)),n,1) ) = A;
为了比赛 - Benchmark
function [t] = bench()
A = magic(200);
% functions to compare
fcns = {
@() thewaywewalk(A);
@() lhcgeneva(A);
@() rayryeng(A);
@() rlbond(A);
};
% timeit
t = zeros(4,1);
for ii = 1:10;
t = t + cellfun(@timeit, fcns);
end
format long
end
function Z = thewaywewalk(A)
n = size(A,2);
rows = 1:numel(A);
cols = mod(rows-1,n) + 1;
Z = zeros(n*n,n);
linIdx = sub2ind(size(Z),rows,cols);
Z(linIdx) = A.';
end
function Z = lhcgeneva(A)
sz = size(A);
Z = zeros(sz(1)*sz(2), sz(2));
for i = 1 : sz(1)
Z((i-1)*sz(2)+1:i*sz(2), :) = diag(A(i, :));
end
end
function Z = rayryeng(A)
A = A.';
Z = full(sparse(1:numel(A), repmat(1:size(A,2),1,size(A,1)), A(:)));
end
function Z = rlbond(A)
D = cellfun(@diag,mat2cell(A, ones(size(A,1), 1), size(A,2)), 'UniformOutput', false);
Z = vertcat(D{:});
end
ans =
0.322633905428601 %// thewaywewalk
0.550931853207228 %// lhcgeneva
0.254718792359946 %// rayryeng - Winner!
0.898236688657039 %// rlbond
此代码将 A 转换为行向量元胞数组,对每个行向量应用 diag
函数,然后将它们堆叠:
D = cellfun(@diag,mat2cell(A, ones(size(A,1), 1), size(A,2)), 'UniformOutput', false);
B = vertcat(D{:});
结果:
B =
1 0 0
0 1 0
0 0 1
2 0 0
0 2 0
0 0 2
3 0 0
0 3 0
0 0 3
编辑:我想 thewaywewalk 的基准测试只给我留下了可读性论点;)
使用 Beaker 的建议进行编辑:
data = [1 1 1
2 2 2
3 3 3];
sz = size(data);
z = zeros(sz(1)*sz(2), sz(2));
for i = 1 : sz(1)
z((i-1)*sz(2)+1:i*sz(2), :) = diag(data(i, :));
end
这是使用 sparse
and repmat
的另一种方法:
A = [1 2 3; 4 5 6; 7 8 9];
A = A.';
B = full(sparse(1:numel(A), repmat(1:size(A,1),1,size(A,2)), A(:)));
原始矩阵在 A
中,我将其转置以便下一步可以正确展开每个矩阵的行。我使用 sparse
来声明矩阵中的非零值。具体来说,我们看到每行只有一个条目,因此行索引的范围应该从 1 到 A
中的条目数。列从 1 到最后一列波动并重复。 mod
当然是通过 thewaywewalk 的解决方案的方式,但我想使用 repmat
以便这是一个独立于他的方法的解决方案。因此,我们创建了一个向量来访问从 1 到我们拥有的尽可能多的列,并且我们对我们拥有的尽可能多的行重复此操作。这些行和列索引向量将决定非零位置出现的位置。最后,将进入每个非零位置的是 A
的元素,这些元素按照行和列索引向量指定的顺序按行主要顺序展开。
请注意,在repmat
调用中,由于转置操作,调用size
时的行和列被颠倒了。
结果如下,我们得到:
>> B
B =
1 0 0
0 2 0
0 0 3
4 0 0
0 5 0
0 0 6
7 0 0
0 8 0
0 0 9
鉴于上述问题的稀疏性,将矩阵保留为 sparse
形式可能会更快,并且仅在必要时使用 full
进行转换。在两种格式之间进行转换需要时间,因此如果您决定进行基准测试,请考虑到这一点。
另一种方法(对我来说)比上面的任何基准测试方法(对于大矩阵和小矩阵)都快
[m,n] = size(A);
Z(n,m*n) = 0;
for idx = 1:n
Z(idx,((idx-1)*m+1):(idx*m)) = A(:,idx);
end
Z = reshape(Z,m*n,n);