2

假设我有一个包含 3 列的矩阵：c1、c2、c3，并且我想创建一个新矩阵，其中每一列都是该矩阵的两列的任何可能乘积。

因此，如果我有一个包含d列的矩阵，我想创建一个包含d+d(d-1)/2+d列的新矩阵。例如，考虑具有 3 列c1, c2,的矩阵c3。我想创建的矩阵应该有列c1, c2, c3, c1xc2, c2xc3, c1xc3, c1^2,c2^2和c3^2。

有没有有效的方法来做到这一点？

4

3 回答 3

4

我很尴尬地发布这个 - 我确信必须有一个更简单的方法（有一个更简单的方法 - 请参阅我在答案底部的 12 月更新），但这将完成这项工作：

A = [1 2 3; 4 5 6];
n =  size(A, 2);

B = A(:, reshape(ones(n, 1) * (1:n), 1, n^2)) .* repmat(A, 1, n);
Soln = [A, B(:, logical(reshape(tril(toeplitz(ones(n, 1))), 1, n^2)'))];

计算效率不高，因为在该B步骤中我实际上计算了我需要的组合数量的两倍（即我得到 c1.*c1、c1.*c2、c1.*c3、c2.*c1、c2.*c2、 c2.*c3, c3.*c1, c3.*c2, c3.*c3)，然后在第二步中，我只取出我需要的列（例如，我删除了 c3.*c1，因为我已经已经得到 c1.*c3 等等）。

更新：刚出去开车，我想到了一个更好的方法。您只需要构造两个形式的索引向量：I1 = [1 1 1 2 2 3]和I2 = [1 2 3 2 3 3]，然后(A(:, I1) .* A(:, I2))将获得您所追求的所有列产品。我现在不在我的电脑旁，但稍后会回来并制定一个构建索引向量的通用方法。我认为使用该结构可以很容易地完成tril(toeplitz)。干杯。几个小时后会更新。

更新： Rody 的第二个解决方案 (+1) 正是我在之前的更新中想到的，所以我不会费心重复他现在在那里所做的事情。尤达的其实也很整洁，所以又+1。

12 月更新：有趣的是，在这里工作之后，我不得不为自己的研究重新审视这个问题（编写 White 的异方差检验）。我现在实际上更喜欢一种新方法，@slayton 在评论中推荐（有点神秘）。具体来说，使用nchoosek. 我的新解决方案如下所示：

T = 20; K = 4;
X = randi(100, T, K);
Index = nchoosek((1:K), 2);
XAll = [X, X(:, Index(:, 1)) .* X(:, Index(:, 2)), X.^2];

nchoosek准确地产生我们需要快速轻松地构建叉积的指数！

于 2012-10-09T02:05:12.590 回答
2

以下是有点体面的：

B = arrayfun(@(x) circshift(A, [0 -x]), 0:size(A,2)-1, 'UniformOutput', false);
B = cat(2, ones(size(A)), B{:});

C = repmat(A, 1, size(A,2)+1) .* B

然而，这将导致矩阵

c1 c2 c3 (c1.*c1) (c2.*c2) (c3.*c3) (c1.*c2) (c2.*c3) (c3.*c1) (c1.*c3) (c2.*c1) (c3.*c2)

其中的顺序与您要求的不同，并且产品不是唯一的。如果您只想要所有独特的产品，请使用：

[sA1, sA2] = size(A);

aa = repmat(1:sA2, sA2,1);

C = [A, A(:,nonzeros(triu(aa))).*A(:,nonzeros(triu(aa.')))]

这导致

c1 c2 c3 (c1.*c1) (c2.*c1) (c2.*c2) (c3.*c1) (c3.*c2) (c3.*c3)

这与您要求的顺序不同，但它仅包含独特的产品。

顺序对您的目的重要吗？

于 2012-10-09T04:38:21.620 回答
1

这是另一种选择。首先，定义一个函数，为给定的列数返回所有可能的对（根据您在问题中的要求）：

cols=@(n)cat(1,num2cell((1:n)'),num2cell((1:n)'*[1,1],2),num2cell(nchoosek(1:n,2),2))

该功能应该是不言自明的。尝试查看输出中的几个小值，n然后亲自查看。有了这个，您可以进行如下操作：

s = RandStream('twister','Seed',1); %for reproducibility 
x = rand(s, 4, 3) %your matrix

%    0.4170    0.1468    0.3968
%    0.7203    0.0923    0.5388
%    0.0001    0.1863    0.4192
%    0.3023    0.3456    0.6852

o = cellfun(@(c)prod(x(:,c),2),cols(size(x,2)),'UniformOutput',0);
out = cat(2,o{:})

%    0.4170    0.1468    0.3968    0.1739    0.0215    0.1574    0.0612    0.1655    0.0582
%    0.7203    0.0923    0.5388    0.5189    0.0085    0.2903    0.0665    0.3881    0.0498
%    0.0001    0.1863    0.4192    0.0000    0.0347    0.1757    0.0000    0.0000    0.0781
%    0.3023    0.3456    0.6852    0.0914    0.1194    0.4695    0.1045    0.2072    0.2368
于 2012-10-09T07:44:00.450 回答