MATLAB:使用块向量高效生成块矩阵
MATLAB: efficient generating of block matrices using a block vector
假设
x = [x1;x2; ...; xn]
其中每个 xi 是长度为 l(i)
的列向量。我们可以设置L = sum(l)
,总长度x
。我想根据 x:
生成 2 个矩阵
我们称它们为 A
和 B
。例如,当 x
仅作为 2 个块 x1
和 x2
则:
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];
在问题的符号中,A
总是L
被L
和B
总是L
被n
。我可以使用循环生成给定 x
的 A
和 B
,但这很乏味。是否有一种聪明的(无循环)方式来生成 A
和 B
。我使用的是 MATLAB 2018b,但如有必要,您可以使用较早版本的 MATLAB。
以下应该有效。在这种情况下,我对元胞数组的转换效率很低,因此可能会有更高效的实现。
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';
编辑
在 运行 一些基准测试之后,我发现基于直接循环的实现速度大约是上述基于单元格的方法的两倍。
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
这是另一种方法:
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{:});
我觉得既短又快:
B = x .* (repelem((1:numel(l)).',l)==(1:numel(l)));
A = B * B.';
如果你有大数据最好使用稀疏矩阵:
B = sparse(1:numel(x), repelem(1:numel(l), l), x);
A = B * B.';
假设
x = [x1;x2; ...; xn]
其中每个 xi 是长度为 l(i)
的列向量。我们可以设置L = sum(l)
,总长度x
。我想根据 x:
我们称它们为 A
和 B
。例如,当 x
仅作为 2 个块 x1
和 x2
则:
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];
在问题的符号中,A
总是L
被L
和B
总是L
被n
。我可以使用循环生成给定 x
的 A
和 B
,但这很乏味。是否有一种聪明的(无循环)方式来生成 A
和 B
。我使用的是 MATLAB 2018b,但如有必要,您可以使用较早版本的 MATLAB。
以下应该有效。在这种情况下,我对元胞数组的转换效率很低,因此可能会有更高效的实现。
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';
编辑
在 运行 一些基准测试之后,我发现基于直接循环的实现速度大约是上述基于单元格的方法的两倍。
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
这是另一种方法:
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{:});
我觉得既短又快:
B = x .* (repelem((1:numel(l)).',l)==(1:numel(l)));
A = B * B.';
如果你有大数据最好使用稀疏矩阵:
B = sparse(1:numel(x), repelem(1:numel(l), l), x);
A = B * B.';