【问题标题】:Add a diagonal of zeros to a matrix in MATLAB将零对角线添加到MATLAB中的矩阵
【发布时间】:2018-09-21 22:49:49
【问题描述】:

假设我在 MATLAB 中有一个维度为 Nx(N-1) 的矩阵 A,例如

N=5;
A=[1  2  3  4;
   5  6  7  8;
   9  10 11 12;
   13 14 15 16;
   17 18 19 20 ];

我想将A 转换为NxN 矩阵B,只需添加一个零对角线,即,

B=[ 0  1   2   3   4;
    5  0   6   7   8;
    9  10  0   11  12;
    13 14  15  0   16;
    17 18  19  20  0];

这段代码做我想做的事:

B_temp = zeros(N,N); 
B_temp(1,:) = [0 A(1,:)];
B_temp(N,:) = [A(N,:) 0];
for j=2:N-1
    B_temp(j,:)= [A(j,1:j-1) 0 A(j,j:end)];
end
B = B_temp; 

您能否提出一种有效的矢量化方法?

【问题讨论】:

    标签: matlab matrix vectorization


    【解决方案1】:

    您可以使用矩阵的上下三角形部分(triutril)来执行此操作。

    那就是 1 行解决方案:

    B = [tril(A,-1) zeros(N, 1)] + [zeros(N,1) triu(A)];
    

    编辑:基准测试

    这是循环方法、Sardar's answer中的2种方法和我上面的方法的对比。

    基准代码,使用timeit 进行计时并直接从问答中提取代码:

    function benchie()
        N = 1e4; A = rand(N,N-1); % Initialise large matrix
        % Set up anonymous functions for input to timeit
        s1 = @() sardar1(A,N); s2 = @() sardar2(A,N); 
        w =  @() wolfie(A,N); u = @() user3285148(A,N);
        % timings
        timeit(s1), timeit(s2), timeit(w), timeit(u)
    end
    function sardar1(A, N) % using eye as an indexing matrix
        B=double(~eye(N)); B(find(B))=A.'; B=B.';
    end
    function sardar2(A,N) % similar to sardar1, but avoiding slow operations
        B=1-eye(N); B(logical(B))=A.'; B=B.';
    end
    function wolfie(A,N) % using triangular parts of the matrix
        B = [tril(A,-1) zeros(N, 1)] + [zeros(N,1) triu(A)];
    end
    function user3285148(A, N) % original looping method
        B = zeros(N,N); B(1,:) = [0 A(1,:)]; B(N,:) = [A(N,:) 0];
        for j=2:N-1; B(j,:)= [A(j,1:j-1) 0 A(j,j:end)]; end
    end
    

    结果:

    • Sardar 方法 1:2.83 秒
    • Sardar 方法 2:1.82 秒
    • 我的方法:1.45 秒
    • 循环方法:3.80 秒 (!)

    结论:

    • 您希望将其矢量化是有根据的,循环比其他方法慢得多。
    • 避免对大型矩阵进行数据转换和find 很重要,Sardar 的方法之间可以节省约 35% 的处理时间。
    • 通过避免同时编制索引,您可以进一步节省 20% 的处理时间。

    【讨论】:

      【解决方案2】:

      生成一个矩阵,其中对角线处为 0,非对角线处为 1。用 A 的转置替换非对角元素(因为 MATLAB 是主要列)。再次转置以获得正确的顺序。

      B = double(~eye(N));  %Converting to double since we want to replace with double entries
      B(find(B)) = A.';     %Replacing the entries
      B = B.';              %Transposing again to get the matrix in the correct order
      

      编辑:

      作为suggestedWolfie的相同算法,您可以摆脱对double的转换和find的使用:

      B = 1-eye(N);
      B(logical(B)) = A.'; 
      B = B.';
      

      【讨论】:

      • 您也可以类似地使用B = 1-eye(N); 来避免与逻辑矩阵之间的转换,然后使用B(logical(B)) = A.'; B = B.'; 来避免使用find。不确定哪个更有效。
      • @Wolfie 好点,我认为摆脱转换和find,处理逻辑索引会更有效。谢谢
      • 效率现在包含在我的答案中:)
      • @Wolfie tbh 我没想到矩阵索引会比 tril triu
      • 哈哈不知道你有没有看过!是的,同样,这就是为什么我很好奇做一个长凳!不确定triu/tril 内部会发生什么,但它们似乎是已编译的函数,因此后端发生了一些快速的事情,然后是一个简单的添加。
      【解决方案3】:

      如果你想在矩阵的对角线上插入任何向量,可以使用普通索引。以下 sn-p 为您提供所需对角线的索引,给定方阵 n 的大小(矩阵为 n by n),以及对角线的数量 k,其中 k=0对应于主对角线,正数k 对应上对角线,负数k 对应下对角线。 ixd 终于给你二维索引了。

      function [idx] = diagidx(n,k)
      % n size of square matrix
      % k number of diagonal
      if k==0 % identity
          idx = [(1:n).' (1:n).']; % [row col]
      elseif k>0 % Upper diagonal
          idx = [(1:n-k).' (1+k:n).'];
      elseif k<0 % lower diagonal
          idx = [(1+abs(k):n).' (1:n-abs(k)).'];
      end
      end
      

      用法:

      n=10;
      k=3;
      A = rand(n);
      idx = diagidx(n,k);
      
      A(idx) = 1:(n-k);
      

      【讨论】:

        猜你喜欢
        • 2019-03-10
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        相关资源
        最近更新 更多