Say I have a vector Q = [Q1 Q2 .... QN]
.
I would like to create a matrix A such that the kth "ring" of the matrix is equal to Qk, with the following constraint:
if N
is odd, the central patch is composed of one number, which is QN
For Q = [12 3 27]
this would be :
A =
12 12 12 12 12
12 3 3 3 12
12 3 27 3 12
12 3 3 3 12
12 12 12 12 12
if N
is even, the central patch is a 2x2 patch where QN
gets repeated
for Q = [12 3]
this would be
A =
12 12 12 12
12 3 3 12
12 3 3 12
12 12 12 12
Two for loops
Two for loops work but it is too slow (~13,3s for 5000x5000 matrices) (Code below) :
%% Two for loops :
% Generate random integer vector Q with unique values
N = 5;
n = 15 * N;
Q = randperm(n,N).';
% Double for loop method
if mod(N,2)==1
mSize = 2*N-1;
else
mSize = 2*N;
end
A = zeros(mSize);
for ii=1:(mSize)
for jj=1:(mSize)
IDX = min([ii,jj,mSize-ii+1,mSize-jj+1]);
A(ii,jj) = Q(IDX);
end
end
Faster approach
I have found a faster approach, which is pretty good (~1.46s for 5000x5000 matrices) but there might still be some room for improvement :
if mod(N,2)==1
mSize = 2*N-1;
I_idx = (1:mSize)-N;
A_fast = Q(end-max(abs(I_idx.'),abs(I_idx)));
else
I_idx = [(N-1):-1:0 0:(N-1)];
A_fast = Q(end-max(I_idx.',I_idx));
end
Any ideas?
The logic of the code is lightly simpler if you follow the advice in Wolfie's comment, and compute only one quadrant that you repeat:
I_idx = 1:N;
B = Q(min(I_idx,I_idx.'));
if mod(N,2)==1
B = [B,B(:,end-1:-1:1)]; % same as [B,fliplr(B(:,1:end-1))]
B = [B;B(end-1:-1:1,:)]; % same as [B;flipud(B(1:end-1,:))]
else
B = [B,fliplr(B)];
B = [B;flipud(B)];
end
This is 2-2.5 times as fast depending on whether Q
is even or odd-sized.
Steve's comment suggest building a triangle first, but I don't see that being any faster, due to the complexity of indexing a matrix's upper or lower triangle.
Testing code:
N = 5000;
n = 15 * N;
Q = randperm(n,N).';
tic
if mod(N,2)==1
mSize = 2*N-1;
I_idx = (1:mSize)-N;
A = Q(end-max(abs(I_idx.'),abs(I_idx)));
else
I_idx = [(N-1):-1:0 0:(N-1)];
A = Q(end-max(I_idx.',I_idx));
end
toc
tic
I_idx = 1:N;
B = Q(min(I_idx,I_idx.'));
if mod(N,2)==1
B = [B,B(:,end-1:-1:1)];
B = [B;B(end-1:-1:1,:)];
else
B = [B,fliplr(B)];
B = [B;flipud(B)];
end
toc
isequal(A,B)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With