用同心 "rings" 值构建矩阵的最快方法

Question

假设我有一个向量 Q = [Q1 Q2 .... QN]。

我想创建一个矩阵 A，使矩阵的第 k 个“环”等于 Qk，具有以下约束：

若N为奇数，则中心块由一个数字组成，即QN

对于 Q = [12 3 27] 这将是：

 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

如果 N 是偶数，中央补丁是一个 2x2 补丁，其中 QN 重复

对于 Q = [12 3] 这将是

A =

12    12    12    12    
12     3     3    12 
12     3     3    12  
12    12    12    12

两个for循环

两个 for 循环有效，但速度太慢（5000x5000 矩阵约 13.3 秒）（代码如下）：

%% 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

更快的方法

我找到了一种更快的方法，非常好（对于 5000x5000 矩阵约为 1.46 秒）但可能仍有一些改进空间：

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

有什么想法吗？

Answer 1

我想出了一个使用repmat的解决方案，然后沿着对角线翻转得到四分之一的解决方案，最后翻转和反转两次得到完整的输出矩阵。

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

可以通过一些巧妙的转置来进行改进以提高速度flips/reverses。

对于更新问题中的偶数情况，以 half_A 和 A 开头的行中的索引应为 end:-1:1 而不是 end-1:-1:1。

运行一些快速的计时，看起来我的解决方案与您更快的方法相当（有时稍慢）：

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)

结果：

double loop
   14.4767

faster approach
    1.8137

flip_it_and_reverse_it
    1.6556

注意：有时 faster_approach 获胜，有时 flip - 我的笔记本电脑上还有一些其他工作运行。

Answer 2

如果您遵循建议，并且只计算您重复的一个象限，代码的逻辑会稍微简单一些：

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

这是快 2-2.5 倍，具体取决于 Q 是偶数还是 odd-sized。

建议先构建一个三角形，但由于索引矩阵的上三角或下三角的复杂性，我认为这样不会更快。

测试代码：

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)

用同心 "rings" 值构建矩阵的最快方法

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

performance

matlab

matrix