logo标签列表

Matlab中mvregress出现“内存不足”错误

matlablinear-regressionmultivariate-testing
3
我正在尝试使用mvregress处理我的数据，其维度为几百个（3~4）。使用32GB的内存，我无法计算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
Anoosh，从你对Daniel问题的回答中，不清楚你是否查找了正确的值。它不是您在任务管理器中看到的系统内存量。请在MATLAB命令窗口中键入“memory”命令，它应该类似于以下内容：
memory 最大可能数组：12699 MB（1.332e + 010字节）* 所有数组可用的内存：12699 MB（1.332e + 010字节）* MATLAB使用的内存：710 MB（7.445e + 008字节） 物理内存（RAM）：8098 MB（8.491e + 009字节）
 - Ingo Schalk-Schupp
哦，我不知道。你能否提供在出现错误时“whos”的结果？这样我们就可以得到所需空间的数量级线索。我猜它是在函数内发生的，对吗？那么我认为n等于400，对吗？可能有一些地方可以节省内存。 - Ingo Schalk-Schupp
@Douba 谢谢！那是一些有用的信息。我希望你能找到问题所在！ - PickleRick
1
线性回归模型的一个经典假设是严格外生性，即 E[e|X] = 0。考虑到您的设置中 N 具有非零均值，除非将问题重新框定为满足严格外生性的问题（例如通过添加每个混合物的指示变量来捕获混合物特定的平均值），否则您将无法一致地估计 A。 - Matthew Gunn
1
如果您想让 n = nsamp 和 d = dim，则几乎肯定要使用 [B，y_hat] = mlrtrain（X'，Y'） 而不是当前的 [B，y_hat] = mlrtrain（X，Y）。 - Matthew Gunn
显示剩余12条评论
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 x n和p x q矩阵时的结果大小为mxp x nxq，在您的情况下为400 x 1001和1000 x 1000，这使得400000 x 1001000矩阵具有4 * 10 ^ 11个元素。现在你有四百个，每个元素占用8个字节的双精度，总大小约为1.281 PB（或1.138 Pebibytes，如果您喜欢），即使使用您宏伟的32 GiB也无法达到。
由于您的矩阵之一——眼睛矩阵——主要包含零，而生成的矩阵包含所有可能的元素乘积组合，其中大多数也将是零。对于这种情况，MATLAB提供了稀疏矩阵格式，它通过仅存储非零元素来节省大量内存，具体取决于矩阵中零元素的数量。您可以使用sparse（X）将完整矩阵转换为稀疏表示，或者直接使用speye（n）获取眼睛矩阵，这就是我在上面所做的。稀疏属性传播到结果，现在您应该有足够的内存了（即使我只有您可用内存的1/4，它也有效）。
然而，马修·冈在评论中提到的问题仍然存在。我收到一个错误消息，内容如下： 使用mvregress（第260行）时出错 估计完整或最小二乘模型的数据不足。
那是一个聪明的技巧！但不幸的是，没有找到答案。因此，基本上我们需要更多的观察才能找到解决方案。如果我们增加观察大小，又需要更多的内存。有没有办法找出解决这个问题需要多少内存？如果我们直接使用A_hat=mvregress(X',Y');呢？它可以被优化吗？ - PickleRick
2
@Anoosh 你可以考虑使用降维技术，而不是保留所有的400个回归器。本论文考虑了在具有分类预测变量的多元回归中的充分维度约减。它可能涉及到更通用的例子。 - Oleg
2
不用客气。也许对于你正在尝试解决的实际数学/算法问题，你可以在适当的论坛上提出另一个问题。内存问题本身似乎已经解决了，对吧？ - Ingo Schalk-Schupp
2
几乎可以确定存在矩阵转置错误，也就是说，在执行Xcell操作之前应该有一个X = X'; Y = Y'。n应该等同于nsamp，d应该等同于dim。这就是为什么它会显示数据不足以估计模型，并且添加观测值似乎会使情况变得更糟（因为在某种程度上，你混淆了观测值和回归器）。 - Matthew Gunn
2

前言

如果你的回归变量在每个回归方程中都相同，并且你对OLS估计感兴趣，可以用一个简单的\调用取代mvregress。

在调用mlrtrain时出现了矩阵转置错误（已经纠正）。在mvregress的语言中，n是观测值的数量，d是结果变量的数量。你生成了一个d乘以n的矩阵Y。但是当你调用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 x dim的数据矩阵X。
  • 一个大小为nsamp x ny的结果变量矩阵Y
  • 你想得到每个Y列在数据矩阵X上的回归结果。也就是说，我们正在进行多元回归，但有一个共同的数据矩阵X。

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

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

如果您想做一些与此不同的事情，您需要清楚地说明您要做什么！
要在X上回归Y，您可以执行以下操作：
[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

如果您想用一个全为1的向量构建一个新的数据矩阵，您可以执行以下操作：

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

一些人可能会感到困惑，这里进行进一步的澄清

假设我们想要运行回归分析

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

我们可以通过一个 n 行 k 列的数据矩阵 X 和一个 n 行 1 列的结果向量 y 来构建数据矩阵。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 的第一列是将 Y(:,1) 回归到 X 上的结果。B 的第二列是将 Y(:,2) 回归到 X 上的结果，以此类推... 在这些条件下，这相当于调用 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 非常不为零的情况下仍然等价。
1
如果我没弄错，你的方法被称为正规方程：beta = pinv(X' * X) * X' * y; 我不确定的是，这个方法是否等同于多元回归。因为我的响应向量y是相关的。 - PickleRick
1
@Anoosh 如果你正在进行普通最小二乘法（即不尝试通过使用残差的协方差结构来提高效率），那么这就是你得到的结果。 - Matthew Gunn
2
多元回归分析，不应与“多变量”混淆，它估计所有Y_i和X_i的联合概率分布，因此从根本上区别于只有一个y和许多回归器X的简单线性回归。参考https://en.wikipedia.org/wiki/Linear_regression#Simple_and_multiple_regression或Greene的计量经济学书籍。这就是为什么没有人质疑OP的意图，他的问题听起来很合理。此外，我们没有足够的信息来判断该方法是否适当，也不应在SO上提出这样的问题。 - Oleg
1
此外，Anoosh的代码片段基本上是从mvregress()中提取出来的。我实际上很惊讶TheMathworks一开始没有考虑构建稀疏矩阵（请参见mvregress的内部代码），并且应该作为该函数的增强功能转发。 - Oleg
1
你可能没有注意到Anoosh提供了一个带有iid样本的虚拟示例，而他的真实数据是相关的。同样，在Greene中有很好的解释，多元回归恰好是n个线性回归，当且仅当Y_i不相关时。我还建议您避免使用色彩丰富的语言。 - Oleg
显示剩余5条评论
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接