MATLAB: efficient generating of block matrices using a block vector

user_6046501 0
Viewed 180

Suppose

x = [x1;x2; ...; xn]

where each xi is a column vector with length l(i). We can set L = sum(l), the total length of x. I would like to generate 2 matrices based on x:

enter image description here

Let's call them A and B. For example, when x only as 2 blocks x1 and x2 then:

A = [x1*x1' zeros(l(1),l(2)); zeros(l(2),l(1)), x2*x2'];

B = [x1 zeros(l(1),1); zeros(l(2),1), x2];

In the notation of the problem, A is always L by L and B is L by n. I can generate A and B given x using loops but it is tedious. Is there a clever (loop-free) way to generate A and B. I am using MATLAB 2018b but you can assume earlier version of MATLAB if necessary.

linear-algebra matlab loops performance
3 Answers

The following should work. In this case I do an inefficient conversion to cell arrays so there may be a more efficient implementation possible.

cuml = [0; cumsum(l(:))];
get_x = @(idx) x((1:l(idx))+cuml(idx));
x_cell = arrayfun(get_x, 1:numel(l), 'UniformOutput', false);
B = blkdiag(x_cell{:});
A = B*B';

Edit

After running some benchmarks I found a direct loop based implementation to be about twice as fast as the cell based approach above.

A = zeros(sum(l));
B = zeros(sum(l), numel(l));
prev = 0;
for idx = 1:numel(l)
    xidx = (1:l(idx))+prev;
    A(xidx, xidx) = x(xidx,1) * x(xidx,1)';
    B(xidx, idx) = x(idx,1);
    prev = prev + l(idx);
end
user_27900470

Here's an alternative approach:

s = repelem(1:numel(l), l).';
t = accumarray(s, x, [], @(x){x*x'});
A = blkdiag(t{:});
t = accumarray(s, x, [], @(x){x});
B = blkdiag(t{:});
user_25869220

I think it is both short and fast:

B = x .* (repelem((1:numel(l)).',l)==(1:numel(l)));
A = B * B.';

If you have large data It is better to use sparse matrix:

B = sparse(1:numel(x), repelem(1:numel(l), l), x);
A = B * B.';
user_65797440
Related