3

我正在尝试将 mvregress 与我拥有的数百维数据一起使用。(3~4)。使用 32 gb 的 ram，我无法计算 beta，并且收到“内存不足”消息。我找不到任何使用 mvregress 的限制阻止我将它应用于具有这种维度的向量，我做错了什么吗？有没有办法通过我的数据使用多元线性回归？

这是一个出错的例子：

dim=400;
nsamp=1000;
dataVariance = .10;
noiseVariance = .05;
mixtureCenters=randn(dim,1);
X=randn(dim, nsamp)*sqrt(dataVariance ) + repmat(mixtureCenters,1,nsamp);
N=randn(dim, nsamp)*sqrt(noiseVariance ) + repmat(mixtureCenters,1,nsamp);
A=2*eye(dim);
Y=A*X+N;
%without residual term:
A_hat=mvregress(X',Y');
%wit residual term:
[B, y_hat]=mlrtrain(X,Y)

在哪里

function [B, y_hat]=mlrtrain(X,Y)
[n,d] = size(Y);
Xmat = [ones(n,1) X];
Xmat_sz=size(Xmat);
Xcell = cell(1,n);
for i = 1:n
    Xcell{i} = [kron([Xmat(i,:)],eye(d))];
end
[beta,sigma,E,V] = mvregress(Xcell,Y);
B = reshape(beta,d,Xmat_sz(2))';
y_hat=Xmat * B ;
end

错误是：

Error using bsxfun
Out of memory. Type HELP MEMORY for your options.

Error in kron (line 36)
   K = reshape(bsxfun(@times,A,B),[ma*mb na*nb]);

Error in mvregress (line 319)
            c{j} = kron(eye(NumSeries),Design(j,:));

这是 whos 命令的结果：

whos
  Name                  Size                Bytes  Class     Attributes

  A                   400x400             1280000  double              
  N                   400x1000            3200000  double              
  X                   400x1000            3200000  double              
  Y                   400x1000            3200000  double              
  dataVariance          1x1                     8  double              
  dim                   1x1                     8  double              
  mixtureCenters      400x1                  3200  double              
  noiseVariance         1x1                     8  double              
  nsamp                 1x1                     8  double
4

2 回答 2

3

好的，我想我有一个解决方案，首先是简短的版本：

dim=400;
nsamp=1000;
dataVariance = .10;
noiseVariance = .05;
mixtureCenters=randn(dim,1);
X=randn(dim, nsamp)*sqrt(dataVariance ) + repmat(mixtureCenters,1,nsamp);
N=randn(dim, nsamp)*sqrt(noiseVariance ) + repmat(mixtureCenters,1,nsamp);
A=2*eye(dim);
Y=A*X+N;

[n,d] = size(Y);
Xmat = [ones(n,1) X];
Xmat_sz=size(Xmat);
Xcell = cell(1,n);
for i = 1:n
    Xcell{i} = kron(Xmat(i,:),speye(d));
end
[beta,sigma,E,V] = mvregress(Xcell,Y);
B = reshape(beta,d,Xmat_sz(2))';
y_hat=Xmat * B ;

奇怪的是，我无法访问函数的工作区，它没有出现在调用堆栈中。这就是我将函数放在脚本之后的原因。

以下是将来可能对您有所帮助的解释：查看kron定义，当插入 m 乘 n 和 ap 乘 q 矩阵时，结果的大小为 mxp 乘 nxq，在您的情况下为 400 乘 1001 和 1000 乘 1000，这使得400000 x 1001000 矩阵，它有 4*10^11 个元素。现在你有 400 个，每个元素占用 8 个字节以实现双精度，即总大小约为 1.281 PB 的内存（或 1.138 Pebibytes，如果你愿意的话），即使你的 32 Gibibyte 也遥不可及.

看到您的矩阵之一，眼睛一，大部分包含零，并且生成的矩阵包含所有可能的元素乘积组合，它们中的大多数也将为零。特别是对于这种情况，MATLAB 提供了稀疏矩阵格式，它通过仅存储非零元素来节省大量内存，具体取决于矩阵中零元素的数量。您可以使用 将完整矩阵转换为稀疏表示sparse(X)，或者使用 直接获得眼睛矩阵speye(n)，这就是我在上面所做的。sparse 属性传播到结果，您现在应该有足够的内存（我有 1/4 的可用内存，它可以工作）。

然而，剩下的是 Matthew Gunn 在评论中提到的问题。我收到一条错误消息：

使用 mvregress 时出错（第 260 行） 数据不足，无法估计完整或最小二乘模型。

于 2016-02-03T16:47:29.147 回答
2

前言

如果您的回归量在每个回归方程中都相同，并且您对 OLS 估计感兴趣，则可以将 mvregress 调用替换为对\.

它在给您的调用中出现mlrtrain了矩阵转置错误（已更正）。在 mvregress 的语言中，n 是观察的数量，d 是结果变量的数量。您生成一个 d 乘以 n 的矩阵 Y。但是 THEN 当你应该调用 mlrtrain(X', Y') 而不是 mlrtrain(X, Y) 时。

如果以下不是具体的，您正在寻找什么，我建议您准确定义您要估计的内容。

如果我是你我会写什么

这里说了这么多完全没有根据，我发布了如果我是你我会写的代码。我已经减少了维度以显示在您的特殊情况下与简单调用的等效性\。我还以更标准的方式编写了一些东西（即让观察沿着行向下运行并且不产生矩阵转置错误）。

dim=5;              % These can go way higher but only if you use my code 
nsamp=20;           % rather than call mvregress
dataVariance = .10;
noiseVariance = .05;
mixtureCenters=randn(dim,1);

X = randn(nsamp, dim)*sqrt(dataVariance )  + repmat(mixtureCenters', nsamp, 1); %'
E = randn(nsamp, dim)*sqrt(noiseVariance);   %noise should be mean zero
B = 2*eye(dim);
Y = X*B+E;
% without constant:
B_hat  = mvregress(X,Y);     %<-------- slow, blows up with high dimension
B_hat2 = X \ Y;              %<-------- fast, fine with higher dimensions

norm(B_hat - B_hat2)         % show numerical equivalent if basically 0

% with constant:
B_constant_hat  = mlrtrain(X,Y)       %<-------- slow, blows up with high dimension
B_constant_hat2 = [ones(nsamp, 1), X] \ Y; % <-- fast, and fine with higher dimensions

norm(B_constant_hat - B_constant_hat2)         % show numerical equivalent if basically 0

解释

我假设你有：

  • 一个nsamp按dim大小排列的数据矩阵X。
  • nsamp按大小排列的ny结果变量矩阵Y
  • Y您希望对on data matrix的每一列进行回归得到结果X。也就是说，我们正在做多元回归，但有一个公共数据矩阵 X。

也就是说，我们正在估计：

y_{ij} = \sum_k b_k * x_{ik} + e_{ijk} for i=1...nsamp, j = 1...ny, k=1...dim

如果您尝试做与此不同的事情，则需要清楚地说明您要做什么！

Y要回归X你可以这样做：

[beta_mvr, sigma_mvr, resid_mvr] = mvregress(X, Y);

这似乎非常缓慢。对于每个回归使用相同数据矩阵的情况，以下内容应与 mvregress 匹配。

beta_hat  = X \ Y;            % estimate beta using least squares
resid     = Y - X * beta_hat;     % calculate residual

如果你想用一个向量构造一个新的数据矩阵，你会这样做：

X_withones = [ones(nsamp, 1), X];

进一步澄清一些困惑

假设我们要运行回归

y_i = \sum_j x_{ij} + e_i  i=1...n, j=1...k

我们可以通过 k 个数据矩阵 X 和一个 n x 1 的结果向量 y 构建数据矩阵 n。OLS 估计值bhat = pinv(X' * X) * X' * y也可以在 MATLAB 中使用bhat = X \ y.

如果您想多次执行此操作（即对同一数据矩阵 X 运行多元回归），您可以构建一个结果矩阵 Y，其中每列代表一个单独的结果变量。Y = [ya, yb, yc, ...]。简单地说，OLS 解决方案B = pinv(X'*X)*X'*Y可以计算为B = X \ Y。B 的第一列是在 X 上回归 Y(:,1) 的结果。B 的第二列是在 X 上回归 Y(:,2) 的结果，依此类推......在这些条件下，这相当于调用 B = mvregress(X, Y)

更多测试代码

如果回归量相同并且通过简单的 OLS 进行估计，则多元回归和方程通过普通最小二乘法方程之间存在等价性。

d = 10;
k = 15;
n = 100;

C = RandomCorr(d + k, 1);  %Use any method you like to generate a random correlation matrix
s = randn(d+k , 1) * 10;
S = (s * s') .* C;         % generate covariance matrix

mu = randn(d+k,1);

data = mvnrnd(ones(n, 1) * mu', S);

Y = data(:,1:d);
X = data(:,d+1:end);

[b1, sigma] = mvregress(X, Y);
b2 = X \ Y;

norm(b1 - b2)

你会注意到 b1 和 b2 在数值上是等价的。即使 sigma 与零相差甚远，它们也是等价的。

于 2016-02-03T21:04:28.297 回答