"out of memory" matlab 中的 mvregress 错误
"out of memory" error for mvregress in matlab
我正在尝试将 mvregress 与我拥有的数百维数据一起使用。 (3~4)。使用 32 gb 的 ram，我无法计算 beta，我收到 "out of memory" 消息。我找不到任何 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
好的，我想我有一个解决方案，首先是简短的版本：
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 和 p×q 矩阵的结果大小为 mxp×nxq，在您的情况下为 400×1001 和 1000×1000，这构成了一个 400000×1001000 矩阵，它具有4*10^11 个元素。现在你有 400 个，每个元素占用 8 个字节以实现双精度，总内存大小约为 1.281 PB（如果你愿意，也可以是 1.138 Pebibytes），即使你的大 32 Gibibyte 也无法达到.
看到你的一个矩阵，眼睛一，主要包含零，并且生成的矩阵包含所有可能的元素乘积组合，其中大部分也将是零。特别是对于这种情况，MATLAB 提供了稀疏矩阵格式，它通过仅存储非零元素来节省大量内存，具体取决于矩阵中零元素的数量。您可以使用 sparse(X) 将完整矩阵转换为稀疏表示，或者使用 speye(n) 直接获得眼图矩阵，这就是我在上面所做的。稀疏 属性 传播到结果，您现在应该有足够的内存（我有 1/4 的可用内存，并且它有效）。
然而，剩下的就是Matthew Gunn在评论中提到的问题。我收到一条错误消息：
使用 mvregress 时出错（第 260 行）
数据不足，无法估计完整模型或最小二乘模型。
前言
如果您的回归量在每个回归方程中都相同并且您对 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 dim 大小的数据矩阵 X。- 一个
nsamp × ny 大小的结果变量矩阵 Y
- 您想要对数据矩阵
X 的 Y 的每一列进行回归的结果。也就是说，我们正在进行多元回归，但有一个共同的数据矩阵 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
如果你想用一个向量构造一个新的数据矩阵，你会这样做：
X_withones = [ones(nsamp, 1), X];
进一步澄清一些困惑
假设我们想要运行回归
y_i = \sum_j x_{ij} + e_i i=1...n, j=1...k
我们可以构建数据矩阵 n x k 数据矩阵 X 和 n x 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的第一列是在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 与零有很大不同，它们也是等价的。
我正在尝试将 mvregress 与我拥有的数百维数据一起使用。 (3~4)。使用 32 gb 的 ram，我无法计算 beta，我收到 "out of memory" 消息。我找不到任何 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
好的，我想我有一个解决方案，首先是简短的版本：
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 和 p×q 矩阵的结果大小为 mxp×nxq，在您的情况下为 400×1001 和 1000×1000，这构成了一个 400000×1001000 矩阵，它具有4*10^11 个元素。现在你有 400 个，每个元素占用 8 个字节以实现双精度，总内存大小约为 1.281 PB（如果你愿意，也可以是 1.138 Pebibytes），即使你的大 32 Gibibyte 也无法达到.
看到你的一个矩阵，眼睛一，主要包含零，并且生成的矩阵包含所有可能的元素乘积组合，其中大部分也将是零。特别是对于这种情况，MATLAB 提供了稀疏矩阵格式，它通过仅存储非零元素来节省大量内存，具体取决于矩阵中零元素的数量。您可以使用 sparse(X) 将完整矩阵转换为稀疏表示，或者使用 speye(n) 直接获得眼图矩阵，这就是我在上面所做的。稀疏 属性 传播到结果，您现在应该有足够的内存（我有 1/4 的可用内存，并且它有效）。
然而，剩下的就是Matthew Gunn在评论中提到的问题。我收到一条错误消息：
使用 mvregress 时出错（第 260 行） 数据不足，无法估计完整模型或最小二乘模型。
前言
如果您的回归量在每个回归方程中都相同并且您对 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
说明
我假设你有：
nsampdim大小的数据矩阵X。- 一个
nsamp×ny大小的结果变量矩阵Y - 您想要对数据矩阵
X的Y的每一列进行回归的结果。也就是说，我们正在进行多元回归，但有一个共同的数据矩阵 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
如果你想用一个向量构造一个新的数据矩阵，你会这样做：
X_withones = [ones(nsamp, 1), X];
进一步澄清一些困惑
假设我们想要运行回归
y_i = \sum_j x_{ij} + e_i i=1...n, j=1...k
我们可以构建数据矩阵 n x k 数据矩阵 X 和 n x 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的第一列是在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 与零有很大不同，它们也是等价的。