Ошибка "недостаточно памяти" для mvregress в matlab
Я пытаюсь использовать mvregress с данными, имеющимися у меня с размерностью в несколько сотен. (3~4). Используя 32 ГБ оперативной памяти, я не могу вычислить бета-версию, и я получаю сообщение "недостаточно памяти". Я не смог найти никаких ограничений использования 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
2 ответа
Хорошо, я думаю, у меня есть решение для вас, сначала короткая версия:
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 на 1001000, которая имеет 4*10^11 элементов. Теперь у вас их четыреста, и каждый элемент занимает 8 байтов для двойной точности, то есть общий объем около 1,281 петабайта памяти (или 1,138 Pebibytes, если вы предпочитаете), что недосягаемо даже для ваших больших 32 Gibibyte,
Видя, что одна из ваших матриц, глазная, содержит в основном нули, а полученная матрица содержит все возможные комбинации элементов, большинство из них также будут равны нулю. В частности, для таких случаев MATLAB предлагает формат разреженной матрицы, который экономит много памяти в зависимости от количества нулевых элементов в матрице, сохраняя только ненулевые элементы. Вы можете преобразовать полную матрицу в разреженное представление с sparse(X)или вы получаете матрицу глаз непосредственно с помощью speye(n), что я и сделал выше. Свойство sparse распространяется на результат, для которого теперь у вас должно быть достаточно памяти (у меня есть 1/4 доступной памяти, и это работает).
Тем не менее, остается проблема, упомянутая Мэтью Ганном в комментарии. Я получаю сообщение об ошибке:
Ошибка при использовании mvregress (строка 260). Недостаточно данных для оценки моделей с полными или наименьшими квадратами.
Предисловие
Если ваши регрессоры одинаковы для каждого уравнения регрессии и вас интересует оценка OLS, вы можете заменить вызов mvregress простым вызовом \,
Появляется в звонке mlrtrain у вас была ошибка транспонирования матрицы (так как исправлено). На языке mvregress n - количество наблюдений, d - количество переменных результата. Вы генерируете матрицу Y, которая d на n. Но ТОГДА, когда вы должны называть 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на матрице данных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
Мы можем построить матрицу данных 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 численно эквивалентны. Они эквивалентны, даже если сигма Крайне отличается от нуля.