Fastest way to build matrix with concentric "rings" of values-CodePudding

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?

CodePudding user response：

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)

CodePudding user response：

I came up with a solution using repmat, then flipping along a diagonal to get a quarter of the solution, finally flipping and reversing twice to get the full output matrix.

function A = flip_it_and_reverse_it(Q)
N = length(Q);
QQ = repmat(Q(:), 1, N);
quarter_A = triu(QQ)   triu(QQ, 1).';
half_A = [quarter_A, quarter_A(:, end-1:-1:1)];
A = [half_A; half_A(end-1:-1:1, :)];
end

There may be improvements that can be made to get faster flips/reverses with some clever transposes.

For the even case in your updated question, the indices in the lines starting half_A and A should be end:-1:1 instead of end-1:-1:1.

Running some quick timings, it looks like my solution is comprable (sometimes slightly slower) to your faster approach:

N = 5000;
n = 15 * N;
Q = randperm(n,N).';

disp('double loop')
tic
double_loop(Q);
disp(toc)

disp('faster approach')
tic
faster_approach(Q);
disp(toc)

disp('flip_it_and_reverse_it')
tic
flip_it_and_reverse_it(Q);
disp(toc)

Results:

double loop
   14.4767

faster approach
    1.8137

flip_it_and_reverse_it
    1.6556

Note: sometimes faster_approach wins, sometimes flip - I've got some other jobs running on my laptop.