在MATLAB中，向矩阵添加一条对角线零

Question

假设我在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; 

你能否提供一个有效的向量化方法？

- TEX

Answer 1

你可以使用矩阵的上三角部分和下三角部分（triu 和 tril）来完成此操作。

然后，这就是一个一行的解决方案：

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

---

编辑：基准测试

这是循环方法、Sardar答案中的两种方法以及我上面的方法之间的比较。

使用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%的处理时间。

- Wolfie

Answer 2

生成一个对角线为零，非对角线为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

编辑：

如同建议的，由于相同的算法，你可以摆脱转换为double和使用find，具体方法见Wolfie。

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

- Sardar Usama

---

你也可以类似地执行 B = 1-eye(N); 来避免转换为逻辑矩阵，然后执行 B(logical(B)) = A.'; B = B.'; 来避免使用 find。不确定哪种方法更高效。 - Wolfie

@Wolfie 很好的建议，我认为摆脱转换和 find，并处理逻辑索引会更有效率。谢谢。 - Sardar Usama

我的答案现在包含了效率优化 :) - Wolfie

@Wolfie 说实话，我没想到矩阵索引比“tril”、“triu”要慢。 - Sardar Usama

哈哈，不知道你有没有看到它！是啊，我也一样，这就是我好奇做一个基准测试的原因！不确定triu/tril内部发生了什么，但它们似乎是编译函数，因此后端会发生一些快速的事情，然后进行简单的加法。 - Wolfie

Answer 3

如果你想在矩阵的对角线上插入任何向量，可以使用普通索引。以下代码片段给出了所需对角线的索引，给定正方形矩阵的大小 n（n x 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);

- Adriaan