8

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?

Improve this question

3 Answers 3

Reset to default
13

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.
Improve this answer
Sign up to request clarification or add additional context in comments.

Comments

5

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.';
Improve this answer

5 Comments

You could also similarly do B = 1-eye(N); to avoid conversion to and from a logical matrix, then B(logical(B)) = A.'; B = B.'; to avoid using find. Not sure which is more efficient.
@Wolfie Good point, I think getting rid of conversion and find, and dealing with logical indices would be more efficient. Thanks
Efficiencies now included in my answer :)
@Wolfie tbh I wasn't expecting matrix indexing to be slower than tril triu
Haha didn't know if you'd seen it! Yeah same, that's why I was curious enough to do a benchie! Not sure what happens inside triu/tril but they appear to be compiled functions so something pacey happening in the backend and then a simple addition.
0

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

Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.