多个常量到矩阵并在matlab中将它们转换为块对角矩阵

Question

我有a1 a2 a3.它们是常数.我有一个矩阵A.我想要做的是获得a1*A,a2*A,a3*A三个矩阵.然后我想将它们转换成对角块矩阵.对于三个常数的情况,这很容易.我可以让b1 = a1*A,b2 = a2*A,b3 = a3*A,然后在matlab中使用blkdiag(b1,b2,b3).

如果我有n个常数,a1 ......一个怎么办？我怎么能在没有任何循环的情况下做到这一点？我知道这可以通过kronecker产品完成,但这非常耗时,你需要做很多不必要的0*常量.

谢谢.

Answer 1

讨论和代码

这可能是一种方法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;

使用方法:使用上面列出的功能代码,让我们假设你有a1,a2,a3,...,an存储在向量a,然后做这样的事情out = bsxfun_linidx(A,a),以在所需的输出out.

标杆

本节将本答案中列出的方法与运行时性能的其他答案中列出的其他两种方法进行比较或基准测试.

其他答案转换为功能形式,如此 -

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

和,

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

这是另一种方法:

1. 使用以下产品作为3D阵列计算产品bsxfun;
2. 转换为每个单元格中包含一个产品(矩阵)的单元阵列;
3. 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