Thread Subject: construct diagonal block matrix

I want to build a diagonal matrix such as

1 2 0 0 0 0
0 0 3 4 0 0
0 0 0 0 5 6

with a given (arbitrary) matrix,
1 2
3 4
5 6

without using loops and cell arrays (conversion takes time) blkdiag works with only parameters (a,b,c,d...) a,b,c,d... are row vectors.
suggestions?

Tim

Hi Tim,

I have a solution. Let

M = [1 2; 3 4; 5 6];

Then, the matrix you are looking for can be computed with the following two lines:
A = [diag(M(:,1)) diag(M(:,2))];
B = A(:,reshape(reshape(1:6,3,2)',1,6))

Hope this helps,
Danny

"Tim Yang" <dlISCool@gmail.com> wrote in message <h0uj19$i02$1@fred.mathworks.com>...
> I want to build a diagonal matrix such as
>
> 1 2 0 0 0 0
> 0 0 3 4 0 0
> 0 0 0 0 5 6
>
> with a given (arbitrary) matrix,
> 1 2
> 3 4
> 5 6
>
> without using loops and cell arrays (conversion takes time) blkdiag works with only parameters (a,b,c,d...) a,b,c,d... are row vectors.
> suggestions?
>
> Tim

"Tim Yang" <dlISCool@gmail.com> wrote in message <h0uj19$i02$1@fred.mathworks.com>...
> I want to build a diagonal matrix such as
>
> 1 2 0 0 0 0
> 0 0 3 4 0 0
> 0 0 0 0 5 6
>
> with a given (arbitrary) matrix,
> 1 2
> 3 4
> 5 6
>
> without using loops and cell arrays (conversion takes time) blkdiag works with only parameters (a,b,c,d...) a,b,c,d... are row vectors.
> suggestions?
>

spdiags(M,[0 1],3,4)

% Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h6kbdj$b2k$1@fred.mathworks.com>...
> "Tim Yang" <dlISCool@gmail.com> wrote in message <h0uj19$i02$1@fred.mathworks.com>...
> > I want to build a diagonal matrix such as
> >
> > 1 2 0 0 0 0
> > 0 0 3 4 0 0
> > 0 0 0 0 5 6
> >
> > with a given (arbitrary) matrix,
> > 1 2
> > 3 4
> > 5 6
> >
> > without using loops and cell arrays (conversion takes time) blkdiag works with only parameters (a,b,c,d...) a,b,c,d... are row vectors.
> > suggestions?
> >
>
> spdiags(M,[0 1],3,4)
>
> % Bruno

Bruno,

That gets close, but is not quite what the OP wants:

full(spdiags(M,[0 1],3,4))

returns

1 2 0 0
0 3 4 0
0 0 5 6

Danny

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h6kbdj$b2k$1@fred.mathworks.com>...

Oopps my soltution won't do. I didn't see you want block shift.

Bruno

Here's a solution that meets the stated requirements, but I don't know if it saves anything. It uses kron, which is generally a slow operation

A=[1 2
   3 4
   5 6];

[m,n]=size(A);

map=logical(kron(speye(m),ones(1,n)));
[M,N]=size(map);

B=spalloc(M,N,numel(A));
At=A.';
B(map)=At(:);

B=full(B),

M = [1 2; 3 4; 5 6];

I=repmat((1:size(M,1)),size(M,2),1);
J=1:numel(M);
Mt=M.';

% sparse
S=sparse(I(:),J(:),Mt(:))
% or full directly
A=accumarray([I(:),J(:)],Mt(:))

% Why should we bother with full matrix in algebra?

Bruno

>
> That gets close, but is not quite what the OP wants:
>
>
> Danny

Yes, yes; Thanks Danny, I saw it after I send the post. I should read more carefully.

Bruno

"Matt " <xys@whatever.com> wrote in message <h6kcbr$c3k$1@fred.mathworks.com>...
> Here's a solution that meets the stated requirements, but I don't know if it saves anything. It uses kron, which is generally a slow operation
>
> A=[1 2
> 3 4
> 5 6];
>
> [m,n]=size(A);
>
> map=logical(kron(speye(m),ones(1,n)));
> [M,N]=size(map);
>
> B=spalloc(M,N,numel(A));
> At=A.';
> B(map)=At(:);
>
> B=full(B),

Matt,

Using a stupid tic-toc, I find that your method takes about .01 sec on my machine. As you though, not that fast. My solution,

M = [1 2; 3 4; 5 6]
A = [diag(M(:,1)) diag(M(:,2))];
B = A(:,reshape(reshape(1:6,3,2)',1,6))

Takes about one-tenth the time and produces the correct result.

Danny

"Matt " <xys@whatever.com> wrote in message <h6kcbr$c3k$1@fred.mathworks.com>...
> Here's a solution that meets the stated requirements, but I don't know if it saves anything. It uses kron, which is generally a slow operation

Actually, it does seem to save something. I did the following timing test on a larger data set. Note that my proposed solution constructs the final matrix B in sparse form, which saves some memory allocation ops. Conversely, blkdiag only returns its output in full form. However even when converting to full in method 1, I still find it to be considerably faster.


A=rand(1000,2);

%method 1
tic
[m,n]=size(A);
map=logical(kron(speye(m),ones(1,n)));
[M,N]=size(map);

B=spalloc(M,N,numel(A));
At=A.';
B(map)=At(:);

%B=full(B);

toc
%Elapsed time is 0.004230 seconds.

%Using blkdiag
tic
 C=num2cell(A,2);
 BB=blkdiag(C{:});

toc
%Elapsed time is 0.050917 seconds.

This should be pretty fast:

A(:,1,:)=diag(M(:,1));
A(:,2,:)=diag(M(:,2));
A=reshape(A,size(M,1),[])

% Bruno

"Matt " <xys@whatever.com> wrote in message <h6kcvo$mkl$1@fred.mathworks.com>...
> "Matt " <xys@whatever.com> wrote in message <h6kcbr$c3k$1@fred.mathworks.com>...
> > Here's a solution that meets the stated requirements, but I don't know if it saves anything. It uses kron, which is generally a slow operation
>
> Actually, it does seem to save something. I did the following timing test on a larger data set. Note that my proposed solution constructs the final matrix B in sparse form, which saves some memory allocation ops. Conversely, blkdiag only returns its output in full form. However even when converting to full in method 1, I still find it to be considerably faster.
>
>
> A=rand(1000,2);
>
> %method 1
> tic
> [m,n]=size(A);
> map=logical(kron(speye(m),ones(1,n)));
> [M,N]=size(map);
>
> B=spalloc(M,N,numel(A));
> At=A.';
> B(map)=At(:);
>
> %B=full(B);
>
> toc
> %Elapsed time is 0.004230 seconds.
>
> %Using blkdiag
> tic
> C=num2cell(A,2);
> BB=blkdiag(C{:});
>
> toc
> %Elapsed time is 0.050917 seconds.

Matt - You're right, I included my solution in the test code:

clc;

A=rand(1000,2);

%method 1
tic
[m,n]=size(A);
map=logical(kron(speye(m),ones(1,n)));
[M,N]=size(map);

B=spalloc(M,N,numel(A));
At=A.';
B(map)=At(:);

B1=full(B);

toc
%Elapsed time is 0.017191 seconds.

%Using blkdiag
tic
 C=num2cell(A,2);
 B2=blkdiag(C{:});

toc
%Elapsed time is 0.035645 seconds.

% Using diag
tic
[m,n]=size(A);
A = [diag(A(:,1)) diag(A(:,2))];
B3 = A(:,reshape(reshape(1:m*n,m,n)',1,m*n));
toc
%Elapsed time is 0.063309 seconds.

All three methods give the same result, but Matt's method 1 is quickest.

Danny

% Sorry make it:

A(:,2,:)=diag(M(:,2));
A(:,1,:)=diag(M(:,1));
A=reshape(A,size(M,1),[]);

% Bruno

"Danny Klein" <djklein.dont.include.this.part@gmail.com> wrote in message <h6kcug$k3q$1@fred.mathworks.com>...

> Using a stupid tic-toc, I find that your method takes about .01 sec on my machine. As you though, not that fast. My solution,
>
> M = [1 2; 3 4; 5 6]
> A = [diag(M(:,1)) diag(M(:,2))];
> B = A(:,reshape(reshape(1:6,3,2)',1,6))
>
> Takes about one-tenth the time and produces the correct result.

Sorry, that's not what I'm finding, at least not on larger data sets. I've posted below a timing test of all 4 methods, including Bruno's, which proves to be the fastest. Yours appears to be slowest, but that might be because it's still using full arrays rather than sparse. I think Bruno's would be the hard to beat however. It's the most direct construction and features the least re-ordering of matrix elements....


A=rand(1000,2);


%Matt - Elapsed time is 0.005380 seconds.
tic

[m,n]=size(A);
map=logical(kron(speye(m),ones(1,n)));
[M,N]=size(map);

B=spalloc(M,N,numel(A));
At=A.';
B(map)=At(:);

B=full(B);

toc
clear B


%blkdiag - Elapsed time is 0.036094 seconds.
tic

 C=num2cell(A,2);
 BB=blkdiag(C{:});

toc
clear C

%Danny - Elapsed time is 0.056357 seconds.
tic
[m,n]=size(A);
M = [diag(A(:,1)) diag(A(:,2))];
J=[1:m;(1:m)+m];
BBB = [M(:,J(:))];
toc
clear I J M BBB

%Bruno - Elapsed time is 0.000348 seconds.
tic
I=repmat((1:size(A,1)),size(A,2),1);
J=1:numel(A);
At=A.';

% sparse
BBBB=sparse(I(:),J(:),At(:));
toc;

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h6ke5n$8sf$1@fred.mathworks.com>...
> This should be pretty fast:
>
> A(:,1,:)=diag(M(:,1));
> A(:,2,:)=diag(M(:,2));
> A=reshape(A,size(M,1),[])
>
> % Bruno

Looks like an improvement to the solution I proposed, unfortunately it is not as fast as Matt's approach:

tic
Z(:,1,:)=diag(A(:,1));
Z(:,2,:)=diag(A(:,2));
B4=reshape(Z,size(A,1),[]);
toc

Gives elapsed time of 0.052722 seconds on my machine whereas Matt's method 1 gives elapsed time of 0.022288 seconds.

Danny

"Danny Klein" <djklein.dont.include.this.part@gmail.com> wrote in message <h6kepr$knh$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h6ke5n$8sf$1@fred.mathworks.com>...
> > This should be pretty fast:
> >
> > A(:,1,:)=diag(M(:,1));
> > A(:,2,:)=diag(M(:,2));
> > A=reshape(A,size(M,1),[])
> >
> > % Bruno
>
> Looks like an improvement to the solution I proposed, unfortunately it is not as fast as Matt's approach:

Danny please reverse the two DIAG lines (avoid reallocation, which is a big drawback).

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <h6kek8$9o4$1@fred.mathworks.com>...
> % Sorry make it:
>
> A(:,2,:)=diag(M(:,2));
> A(:,1,:)=diag(M(:,1));
> A=reshape(A,size(M,1),[]);

Bruno- I think your original proposal would have to be faster than this for large data sets (but don't have time to check right now). Your original proposal works entirely in the sparse matrix domain, thus saving a lot of memory allocation ops, as compared to the above which constructs large full diagonal matrices. Unfortunately, your above proposal can't be modified to work in the sparse domain, because 3D arrays can't be represented as type sparse.

"Matt " <xys@whatever.com> wrote in message <h6kf7g$hq4$1@fred.mathworks.com>...

> Bruno- I think your original proposal would have to be faster than this for large data sets (but don't have time to check right now). Your original proposal works entirely in the sparse matrix domain, thus saving a lot of memory allocation ops, as compared to the above which constructs large full diagonal matrices. Unfortunately, your above proposal can't be modified to work in the sparse domain, because 3D arrays can't be represented as type sparse.

100% agree Matt. And I already tested it.

Bruno

For small problems, it would seem that my original approach works faster than the sparse methods you guys have suggested. I suppose the added overhead of sparse manipulation doesn't pay off. For the size of the matrix in the original post by Tim (3x2), my method gives about 0.000998 sec whereas both of Matt's give about 0.0017 sec. All of these times are super small, so it probably doesn't matter. I think Tim was just looking for _any_ way to compute the answer.

"Danny Klein" <djklein.dont.include.this.part@gmail.com> wrote in message <h6kgbv$mt$1@fred.mathworks.com>...
> For the size of the matrix in the original post by Tim (3x2), my method gives about 0.000998 sec whereas both of Matt's give about 0.0017 sec. All of these times are super small, so it probably doesn't matter. I think Tim was just looking for _any_ way to compute the answer.
----------------

Then it would have made sense just to go with blkdiag()