Multiple matrix constant and converting them to block diagonal matrix in matlab

Question

I have a1 a2 a3. These are constants. I have a matrix A. What I want to do is get a1 * A, a2 * A, a3 * A three matrices. Then I want to transfer them to the matrix of diagonal blocks. For the case of three constants, this is easy. I can enable b1 = a1 * A, b2 = a2 * A, b3 = a3 * A, and then use blkdiag (b1, b2, b3) in matlab.

What if I have n constants, a1 ... an. How can I do this without any loops? I know that this can be done using the kronecker product, but it is very time consuming and you need to make a lot of unnecessary 0 * constant.

Thank.

Answer 1

Discussion and code

This may be one approach with , which makes it easier to function-encoded - bsxfun(@plus linear indexing

function out = bsxfun_linidx(A,a)
%// Get sizes
[A_nrows,A_ncols] = size(A);
N_a = numel(a);

%// Linear indexing offsets between 2 columns in a block & between 2 blocks
off1 = A_nrows*N_a;
off2 = off1*A_ncols+A_nrows;

%// Get the matrix multiplication results
vals = bsxfun(@times,A,permute(a,[1 3 2])); %// OR vals = A(:)*a_arr;

%// Get linear indices for the first block
block1_idx = bsxfun(@plus,[1:A_nrows]',[0:A_ncols-1]*off1);  %//'

%// Initialize output array base on fast pre-allocation inspired by -
%// http://undocumentedmatlab.com/blog/preallocation-performance
out(A_nrows*N_a,A_ncols*N_a) = 0; 

%// Get linear indices for all blocks and place vals in out indexed by them
out(bsxfun(@plus,block1_idx(:),(0:N_a-1)*off2)) = vals;

return;

How to use: To use the above function code, assume that you have a1, a2, a3, ...., anstored in a vector a, then do something like this out = bsxfun_linidx(A,a)to get the desired result out.

Benchmarking

This section compares or compares the approach specified in this answer compared to the other two approaches listed in the other answers for execution at runtime.

Other answers have been converted to function forms, for example:

function B = bsxfun_blkdiag(A,a)
B = bsxfun(@times, A, reshape(a,1,1,[])); %// step 1: compute products as a 3D array
B = mat2cell(B,size(A,1),size(A,2),ones(1,numel(a))); %// step 2: convert to cell array
B = blkdiag(B{:}); %// step 3: call blkdiag with comma-separated list from cell array

and

function out = kron_diag(A,a_arr)
out = kron(diag(a_arr),A);

a a, <

• a 500 x 500 a 1 x 10
• a 200 x 200 a 1 x 50
• a 100 x 100 a 1 x 100
• a 50 x 50 a 1 x 200

-

%// Datasizes
N_a = [10  50  100 200];
N_A = [500 200 100 50];

timeall = zeros(3,numel(N_a)); %// Array to store runtimes
for iter = 1:numel(N_a)

    %// Create random inputs
    a = randi(9,1,N_a(iter));
    A = rand(N_A(iter),N_A(iter));

    %// Time the approaches
    func1 = @() kron_diag(A,a);
    timeall(1,iter) = timeit(func1); clear func1

    func2 = @() bsxfun_blkdiag(A,a);
    timeall(2,iter) = timeit(func2); clear func2

    func3 = @() bsxfun_linidx(A,a);
    timeall(3,iter) = timeit(func3); clear func3
end

%// Plot runtimes against size of A
figure,hold on,grid on
plot(N_A,timeall(1,:),'-ro'),
plot(N_A,timeall(2,:),'-kx'),
plot(N_A,timeall(3,:),'-b+'),
legend('KRON + DIAG','BSXFUN + BLKDIAG','BSXFUN + LINEAR INDEXING'),
xlabel('Datasize (Size of A) ->'),ylabel('Runtimes (sec)'),title('Runtime Plot')

%// Plot runtimes against size of a
figure,hold on,grid on
plot(N_a,timeall(1,:),'-ro'),
plot(N_a,timeall(2,:),'-kx'),
plot(N_a,timeall(3,:),'-b+'),
legend('KRON + DIAG','BSXFUN + BLKDIAG','BSXFUN + LINEAR INDEXING'),
xlabel('Datasize (Size of a) ->'),ylabel('Runtimes (sec)'),title('Runtime Plot')

, ,

: , , bsxfun, , !

Answer 2

:

• , bsxfun;
• () ;
• blkdiag , , .

A , A . B

B = bsxfun(@times, A, reshape(a,1,1,[])); %// step 1: compute products as a 3D array
B = mat2cell(B,size(A,1),size(A,2),ones(1,numel(a))); %// step 2: convert to cell array
B = blkdiag(B{:}); %// step 3: call blkdiag with comma-separated list from cell array

Answer 3

It uses a method kronthat seems to be faster and more memory efficient than a Divakar-based solution bsxfun. I'm not sure if this is different from your method, but the time seems pretty good. It might be worth doing some testing between the various methods that will be more effective for you.

A=magic(4);

a1=1;
a2=2;
a3=3;

kron(diag([a1 a2 a3]),A)