Add a diagonal of zeros to a matrix in MATLAB

Question

Suppose I have a matrix A of dimension Nx(N-1) in MATLAB, e.g.

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

I want to transform A into an NxN matrix B, just by adding a zero diagonal, i.e.,

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];

---

This code does what I want:

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; 

Could you suggest an efficient way to vectorise it?

Answer 1

Generate a matrix with zeros at diagonal and ones at non-diagonal indices. Replace the non-diagonal elements with the transpose of A (since MATLAB is column major). Transpose again to get the correct order.

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

Edit:

As suggested by Wolfie for the same algorithm, you can get rid of conversion to double and the use of find with:

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

Answer 2

If you want to insert any vector on a diagonal of a matrix, one can use plain indexing. The following snippet gives you the indices of the desired diagonal, given the size of the square matrix n (matrix is n by n), and the number of the diagonal k, where k=0 corresponds to the main diagonal, positive numbers of k to upper diagonals and negative numbers of k to lower diagonals. ixd finally gives you the 2D indices.

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

Usage:

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

A(idx) = 1:(n-k);

Answer 3

You can do this with upper and lower triangular parts of the matrix (triu and tril).

Then it's a 1 line solution:

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

---

Edit: benchmark

This is a comparison of the loop method, the 2 methods in Sardar's answer, and my method above.

Benchmark code, using timeit for timing and directly lifting code from question and answers:

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

Results:

• Sardar method 1: 2.83 secs
• Sardar method 2: 1.82 secs
• My method: 1.45 secs
• Looping method: 3.80 secs (!)

Conclusions:

• Your desire to vectorise this was well founded, looping is way slower than other methods.
• Avoiding data conversions and find for large matrices is important, saving ~35% processing time between Sardar's methods.
• By avoiding indexing all together you can save a further 20% processing time.