How to convert matrix to a stack of diagonal matri

问题:

I have a matrix:

A = [1 1 1
     2 2 2
     3 3 3]

Is there a vectorized way of obtaining:

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]

回答1:

Here's another way using 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(:)));

The original matrix is in A, and I transpose it so I can unroll the rows of each matrix properly for the next step. I use sparse to declare what is non-zero in a matrix. Specifically, we see that there is only one entry per row and so the row indices should range from 1 up to as many entries as there are in A. The columns fluctuate from 1 up to the last column and repeat. mod is certainly the way to go via thewaywewalk's solution, but I wanted to use repmat so that this is an independent solution from his approach. As such, we create a vector for accessing the columns that goes from 1 up to as many columns as we have, and we repeat this for as many rows as we have. These row and column index vectors are is going to dictate where the non-zero locations will appear. Finally, what will go into each non-zero location are the elements of A unrolled in row major order, following the order dictated by the row and column index vectors.

Take note that in the repmat call, the rows and columns when calling size are reversed due to the transpose operation.

The result thus follows and we get:

>> 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

Given the sparsity of the above problem, it may be faster to leave the matrix in sparse form and only convert using full if necessary. There will be time spent to convert between the two formats so take that into consideration if you decide to benchmark.

回答2:

This is one solution using mod and sub2ind:

%// 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

---

Or if you prefer saving lines of code, instead of readability:

out = zeros(n*n,n);
out(sub2ind(size(out),1:numel(data),mod((1:numel(data))-1,n) + 1)) = data.'

---

Two other fast solutions, but not faster than the others:

%// #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;

---

For the sake of competition - 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

回答3:

Edit: I guess thewaywewalk's benchmark leaves me only with a readability argument ;)

Edit using Beaker's suggestion:

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

回答4:

An alternative, which (for me) is faster than any of the methods benchmarked above (for large and small matrices) is

[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);

回答5:

This code converts A to a cell array of row vectors, applies the diag function to each, and then stacks them:

D = cellfun(@diag,mat2cell(A, ones(size(A,1), 1), size(A,2)), 'UniformOutput', false);
B = vertcat(D{:});

Results:

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