复制向量将它们向右移动
Replicate vectors shifting them to the right
在 Matlab 中,我在 2x249 矩阵中有两个单行 (1x249) 向量,我必须通过多次复制它们来创建矩阵 A
,每次将 2 个位置的向量向右移动.我想用零填充左侧的条目。有没有聪明的方法来做到这一点?目前,我正在使用 for 循环和 circshift,并且我在每次迭代中添加我将新行添加到 A
,但这可能是非常低效的。
代码(myMat
是我要平移的矩阵):
A = [];
myMat = [1 0 -1 zeros(1,246); 0 2 0 -2 zeros(1,245)];
N = 20;
for i=1:N-1
aux = circshift(myMat,[0,2*(i-1)]);
aux(:,1:2*(i-1)) = 0;
A =[A; aux];
end
您可能已经知道,Matlab 中的循环效率不高。
我知道 Mathworks 一直在说 JIT 不再如此
编译,但我还没有体验过快速循环。
我把你构造矩阵A的方法放在一个函数中:
function A = replvector1(myMat,shift_right,width,N)
pre_alloc = true; % make implementation faster using pre-allocation yes/no
% Pad myMat with zeros to make it wide enough
myMat(1,width)=0;
% initialize A
if pre_alloc
A = zeros(size(myMat,1)*(N-1),width);
else
A = [];
end
% Fill A
for i=1:N-1
aux = circshift(myMat,[0,shift_right*(i-1)]);
aux(:,1:min(width,shift_right*(i-1))) = 0;
A(size(myMat,1)*(i-1)+1:size(myMat,1)*i,:) =aux;
end
你的矩阵运算看起来很像克罗内克积,但是
块矩阵具有重叠的列范围,因此直接克罗内克产品
不管用。相反,我构建了以下函数:
function A = replvector2(myMat,shift_right,width,N)
[i,j,a] = find(myMat);
i = kron(ones(N-1,1),i) + kron([0:N-2]',ones(size(i))) * size(myMat,1);
j = kron(ones(N-1,1),j) + kron([0:N-2]',ones(size(j))) * shift_right;
a = kron(ones(N-1,1),a);
ok = j<=width;
A = full(sparse(i(ok),j(ok),a(ok),(N-1)*size(myMat,1),width));
您可以通过删除分号并查看中间值来遵循该算法
结果。
下面的主程序运行是您的示例,可以很容易地修改为
运行 相似例子:
% inputs (you may vary them to see that it always works)
shift_right = 2;
width = 249;
myMat1 = [ 1 0 -1 0 ;
0 2 0 -2 ];
N = 20;
% Run your implementation
tic;
A = replvector1(myMat,shift_right,width,N);
disp(sprintf('\n original implementation took %e sec',toc))
% Run the new implementation
tic;
B = replvector2(myMat,shift_right,width,N);
disp(sprintf(' new implementation took %e sec',toc))
disp(sprintf('\n norm(B-A)=%e\n',norm(B-A)))
我采用了 Nathan 的代码 (),并添加了另一个可能的实现 (replvector3
)。
我的想法源于你并不真的需要循环转变。您需要右移并在左侧添加零。如果您从一个预先分配的数组开始(这确实是您及时获得大赢家的地方,其余的都是花生),那么您已经有了零。现在您只需将 myMat
复制到正确的位置。
这些是我看到的时间 (MATLAB R2017a):
OP's, with pre-allocation: 1.1730e-04
Nathan's: 5.1992e-05
Mine: 3.5426e-05
^ shift by one on purpose, to make comparison of times easier
这是完整的副本,复制粘贴到 M 文件中 运行:
function so
shift_right = 2;
width = 249;
myMat = [ 1 0 -1 0 ;
0 2 0 -2 ];
N = 20;
A = replvector1(myMat,shift_right,width,N);
B = replvector2(myMat,shift_right,width,N);
norm(B(:)-A(:))
C = replvector3(myMat,shift_right,width,N);
norm(C(:)-A(:))
timeit(@()replvector1(myMat,shift_right,width,N))
timeit(@()replvector2(myMat,shift_right,width,N))
timeit(@()replvector3(myMat,shift_right,width,N))
% Original version, modified to pre-allocate
function A = replvector1(myMat,shift_right,width,N)
% Assuming width > shift_right * (N-1) + size(myMat,2)
myMat(1,width) = 0;
M = size(myMat,1);
A = zeros(M*(N-1),width);
for i = 1:N-1
aux = circshift(myMat,[0,shift_right*(i-1)]);
aux(:,1:shift_right*(i-1)) = 0;
A(M*(i-1)+(1:M),:) = aux;
end
% Nathan's version
function A = replvector2(myMat,shift_right,width,N)
[i,j,a] = find(myMat);
i = kron(ones(N-1,1),i) + kron((0:N-2)',ones(size(i))) * size(myMat,1);
j = kron(ones(N-1,1),j) + kron((0:N-2)',ones(size(j))) * shift_right;
a = kron(ones(N-1,1),a);
ok = j<=width;
A = full(sparse(i(ok),j(ok),a(ok),(N-1)*size(myMat,1),width));
% My trivial version with loops
function A = replvector3(myMat,shift_right,width,N)
% Assuming width > shift_right * (N-1) + size(myMat,2)
[M,K] = size(myMat);
A = zeros(M*(N-1),width);
for i = 1:N-1
A(M*(i-1)+(1:M),shift_right*(i-1)+(1:K)) = myMat;
end
在 Matlab 中,我在 2x249 矩阵中有两个单行 (1x249) 向量,我必须通过多次复制它们来创建矩阵 A
,每次将 2 个位置的向量向右移动.我想用零填充左侧的条目。有没有聪明的方法来做到这一点?目前,我正在使用 for 循环和 circshift,并且我在每次迭代中添加我将新行添加到 A
,但这可能是非常低效的。
代码(myMat
是我要平移的矩阵):
A = [];
myMat = [1 0 -1 zeros(1,246); 0 2 0 -2 zeros(1,245)];
N = 20;
for i=1:N-1
aux = circshift(myMat,[0,2*(i-1)]);
aux(:,1:2*(i-1)) = 0;
A =[A; aux];
end
您可能已经知道,Matlab 中的循环效率不高。 我知道 Mathworks 一直在说 JIT 不再如此 编译,但我还没有体验过快速循环。
我把你构造矩阵A的方法放在一个函数中:
function A = replvector1(myMat,shift_right,width,N)
pre_alloc = true; % make implementation faster using pre-allocation yes/no
% Pad myMat with zeros to make it wide enough
myMat(1,width)=0;
% initialize A
if pre_alloc
A = zeros(size(myMat,1)*(N-1),width);
else
A = [];
end
% Fill A
for i=1:N-1
aux = circshift(myMat,[0,shift_right*(i-1)]);
aux(:,1:min(width,shift_right*(i-1))) = 0;
A(size(myMat,1)*(i-1)+1:size(myMat,1)*i,:) =aux;
end
你的矩阵运算看起来很像克罗内克积,但是 块矩阵具有重叠的列范围,因此直接克罗内克产品 不管用。相反,我构建了以下函数:
function A = replvector2(myMat,shift_right,width,N)
[i,j,a] = find(myMat);
i = kron(ones(N-1,1),i) + kron([0:N-2]',ones(size(i))) * size(myMat,1);
j = kron(ones(N-1,1),j) + kron([0:N-2]',ones(size(j))) * shift_right;
a = kron(ones(N-1,1),a);
ok = j<=width;
A = full(sparse(i(ok),j(ok),a(ok),(N-1)*size(myMat,1),width));
您可以通过删除分号并查看中间值来遵循该算法 结果。
下面的主程序运行是您的示例,可以很容易地修改为 运行 相似例子:
% inputs (you may vary them to see that it always works)
shift_right = 2;
width = 249;
myMat1 = [ 1 0 -1 0 ;
0 2 0 -2 ];
N = 20;
% Run your implementation
tic;
A = replvector1(myMat,shift_right,width,N);
disp(sprintf('\n original implementation took %e sec',toc))
% Run the new implementation
tic;
B = replvector2(myMat,shift_right,width,N);
disp(sprintf(' new implementation took %e sec',toc))
disp(sprintf('\n norm(B-A)=%e\n',norm(B-A)))
我采用了 Nathan 的代码 (replvector3
)。
我的想法源于你并不真的需要循环转变。您需要右移并在左侧添加零。如果您从一个预先分配的数组开始(这确实是您及时获得大赢家的地方,其余的都是花生),那么您已经有了零。现在您只需将 myMat
复制到正确的位置。
这些是我看到的时间 (MATLAB R2017a):
OP's, with pre-allocation: 1.1730e-04
Nathan's: 5.1992e-05
Mine: 3.5426e-05
^ shift by one on purpose, to make comparison of times easier
这是完整的副本,复制粘贴到 M 文件中 运行:
function so
shift_right = 2;
width = 249;
myMat = [ 1 0 -1 0 ;
0 2 0 -2 ];
N = 20;
A = replvector1(myMat,shift_right,width,N);
B = replvector2(myMat,shift_right,width,N);
norm(B(:)-A(:))
C = replvector3(myMat,shift_right,width,N);
norm(C(:)-A(:))
timeit(@()replvector1(myMat,shift_right,width,N))
timeit(@()replvector2(myMat,shift_right,width,N))
timeit(@()replvector3(myMat,shift_right,width,N))
% Original version, modified to pre-allocate
function A = replvector1(myMat,shift_right,width,N)
% Assuming width > shift_right * (N-1) + size(myMat,2)
myMat(1,width) = 0;
M = size(myMat,1);
A = zeros(M*(N-1),width);
for i = 1:N-1
aux = circshift(myMat,[0,shift_right*(i-1)]);
aux(:,1:shift_right*(i-1)) = 0;
A(M*(i-1)+(1:M),:) = aux;
end
% Nathan's version
function A = replvector2(myMat,shift_right,width,N)
[i,j,a] = find(myMat);
i = kron(ones(N-1,1),i) + kron((0:N-2)',ones(size(i))) * size(myMat,1);
j = kron(ones(N-1,1),j) + kron((0:N-2)',ones(size(j))) * shift_right;
a = kron(ones(N-1,1),a);
ok = j<=width;
A = full(sparse(i(ok),j(ok),a(ok),(N-1)*size(myMat,1),width));
% My trivial version with loops
function A = replvector3(myMat,shift_right,width,N)
% Assuming width > shift_right * (N-1) + size(myMat,2)
[M,K] = size(myMat);
A = zeros(M*(N-1),width);
for i = 1:N-1
A(M*(i-1)+(1:M),shift_right*(i-1)+(1:K)) = myMat;
end