Matlab, обратная большой матрице

Вот уравнение, которое я пытаюсь решить:

h = (X'*X)^-1*X'*y

где X — матрица, а y — вектор ((X'X)^-1 является инверсией X-транспонирования, умноженного на X). Я закодировал это в Matlab как:

h = (X'*X)\X'*y

что я считаю правильным. Проблема в том, что X составляет около 10000x10000, и попытка вычислить инверсию приводит к сбою Matlab даже на самом мощном компьютере, который я могу найти (16 ядер, 24 ГБ ОЗУ). Есть ли способ разделить это или библиотеку, предназначенную для таких больших инверсий?

Спасибо.


matlab matrix-inverse
person Jordan    schedule 12.03.2013    source источник
comment
Насколько велика матрица? Можете ли вы показать мне вектор y и матрицу X. Или хотя бы дать мне представление о том, на что это похоже?   -  person Derek W    schedule 13.03.2013
comment
Какова размерность X и y? Вы можете захотеть вычислить приближение обратного, если размерность очень высока. Кстати, ты не скучаешь по некоторым *s?   -  person sfotiadis    schedule 13.03.2013
comment
Отредактированный вопрос с информацией. X составляет около 10 000x10 000. Размерность вектора y совпадает (длиной также около 10 000).   -  person Jordan    schedule 13.03.2013
comment
Вы знаете что-нибудь о структуре матрицы? например, если его редкость и т. д.? Кстати, я считаю, что вместо '\' вы можете использовать функцию inv(), поскольку '\' решает задачу наименьших квадратов (псевдоинверсия)   -  person amas    schedule 13.03.2013
comment
Является редким? вам нужно сохранить точность double? у вас есть доступ к gpu?   -  person bla    schedule 13.03.2013
comment
X — это особая структура, она вроде как 3000-диагональная (первый столбец будет состоять из 3000 ненулевых, затем 4000 нулевых, второй столбец будет из 1 нуля, затем 3000 ненулевых, затем 3999 нулевых и т. д. Чтобы исправить мой предыдущий комментарий, X на самом деле 7000x10000, поэтому X'*X равно 10000x10000.   -  person Jordan    schedule 13.03.2013
comment
Это, вероятно, будет вам полезно... math.stackexchange.com/questions/16940/   -  person bla    schedule 13.03.2013

Ответы (4)


arrow_upward
0
arrow_downward

Я сгенерировал случайную матрицу X размером 10 000 на 10 000 и случайный вектор y размером 10 000 на 1.

Я просто разбил свои вычисления шаг за шагом. (Код показан ниже)

  1. Вычислил транспонирование и сохранил его в матрице K
  2. Затем я вычислил матрицу A, умножив K на X.
  3. Вычисленный вектор b путем умножения K на вектор y
  4. Наконец, я использовал оператор обратной косой черты для A и b, чтобы решить

У меня не было проблем с вычислением. Это заняло некоторое время, но разделение операций на наименьшие возможные группы помогло предотвратить перегрузку компьютера. Однако это может быть состав матрицы, которую вы используете (например, разреженные, десятичные числа и т. д.).

X = randi(2000, [10000, 10000]);
y = randi(2000, 10000, 1);
K = X';
A = K*X;
b = K*y;
S = A\b;
person Derek W    schedule 13.03.2013
comment
randi создает одно значение в диапазоне от 1 до 10 000, а не матрицу 10 000 x 10 000. Окончательное решение должно быть вектором. - person Jordan; 13.03.2013
comment
Моя ошибка! Вместо этого полностью посмотрел синтаксис randint! Я обновлю свой ответ. - person Derek W; 13.03.2013
comment
Спасибо. Однако есть небольшая проблема с тем, как мне нужно реализовать код: мне нужно применить один и тот же X ко многим различным значениям для вектора y, и я не хочу каждый раз инвертировать. Итак, вас нельзя представить до самого последнего шага. - person Jordan; 13.03.2013

arrow_upward
1
arrow_downward

Это похоже на псевдоинверсию. Возможно, вы ищете только

h = X \ y;
person Ben Voigt    schedule 13.03.2013

arrow_upward
0
arrow_downward

Если в вашем распоряжении несколько машин, и вы можете преобразовать свою проблему в форму h = X\y, предложенную @Ben, вы можете использовать распределенные массивы. Эта демонстрация показывает, как вы можете это сделать.

person Edric    schedule 13.03.2013

arrow_upward
0
arrow_downward

Джордан, Ваше уравнение точно определяет «обратную матрицу Мура-Пенроуза».

Проверьте: http://mathworld.wolfram.com/Moore-PenroseMatrixInverse.html

Непосредственное использование h = X \ y; должно помочь. Или проверьте Matlab pinv(X)*y

person Dazhao Liu    schedule 24.06.2013