获取所有矩阵行列的所有双因子产品

Question

我有一个n×m的矩阵。设该矩阵的某一行为x。每一行代表一个数字的功能x1, x2, x3, ...

现在，我想收到上述对角线x * x'，即元素：x1*x2, x1*x3, x2*x3，...但不x1*x1。另外，如果我有x1*x2，我不需要x2*x1。

我想将这些列与产品添加到我的矩阵。鉴于我以前有过m列，我将为这些产品添加更多列，即：(m^2 + m)/2 - m更多列。

这应该为我的矩阵的每一行完成。

我已经在Matlab中找到了解决方案。然而，它似乎很慢，我想知道是否有更多矢量化解决方案可供Matlab使用，可以更快地执行。

我目前的解决方案使用一个包来获得高于上对角线元素的矢量：https://de.mathworks.com/matlabcentral/fileexchange/23391-triangular-and-diagonal-indexing

矩阵M会给我我的矩阵在对角线之上的元素M(itriu(size(M),1))。例如，如果我抛出[1 2 3; 4 5 6; 7 8 9]，我将得到2 3 6。

我的代码如下：

function [ X_out ] = permutateFeatures(X_in) 
%PERMUTATEFEATURES given a matrix with m features in the columns 
% and n samples in the rows, return a [n (m^2 + m)/2] matrix 
% where each additional column contains a element-wise product of two of 
% the original columns 

n = size(X_in, 1); 
m = size(X_in, 2); 

X_out = [X_in zeros(n, (m^2 + m)/2 - m)]; 

for i = 1:n 
    outerProduct = X_out(i,1:m)' * X_out(i,1:m); 
    X_out(i,:) = [X_in(i,:) outerProduct(itriu(size(outerProduct),1))']; 
end 

end 

有没有更有效的解决方案？

Answer 1

这里的一个矢量化溶液 -

[r,c] = find(triu(true(size(X_in,2)),1)); 
out = [X_in X_in(:,r).*X_in(:,c)]; 

运行测试

时序代码 -

% Setup input array 
% (as stated in comments : m is mostly <20. n goes into the millions) 
X_in = randi(5,[50000,20]); 

disp('--------------------------- Original Solution') 
tic, 
n = size(X_in, 1); 
m = size(X_in, 2); 
X_out = [X_in zeros(n, (m^2 + m)/2 - m)]; 
for i = 1:n 
    outerProduct = X_out(i,1:m)' * X_out(i,1:m); 
    X_out(i,:) = [X_in(i,:) outerProduct(itriu(size(outerProduct),1))']; 
end 
toc 

disp('--------------------------- Proposed Solution') 
tic, 
[r,c] = find(triu(true(size(X_in,2)),1)); 
out = [X_in X_in(:,r).*X_in(:,c)]; 
toc, 

计时 -

--------------------------- Original Solution 
Elapsed time is 8.618389 seconds. 
--------------------------- Proposed Solution 
Elapsed time is 0.131146 seconds. 

巨大的加速有60x+！

Answer 2

这里的矩阵乘法是向量化的，这是计算的很大一部分。如果你愿意，你可以向量化VEC 1及VEC 2的创作为好，但只有一点效率刍议获得：

vec1=[]; 
vec2=[]; 
for i = 1:n 
    vec1=[vec1 i*ones(1,n-i)]; 
    vec2=[vec2 (i+1):n]; 
end 
X_out2=[X_in X_in(:,vec1).*X_in(:,vec2)]; 

为rand(1000,1000)老办法，这一次执行

Elapsed time is 24.709988 seconds. 
Elapsed time is 6.753230 seconds. 

在我的机器上，使用相同的解决方案。