Вопрос по sparse-matrix, linear-algebra, python, least-squares, matrix-inverse – @eat вы правы, я переключил V и U. Я думаю, что все остальное правильно, но я согласен с точкой вашего ответа, что псевдообратный, вероятно, не лучший способ решения ее проблемы.

6

отаю с данными из нейровизуализации, и из-за большого объема данных я хотел бы использовать разреженные матрицы для моего кода (scipy.sparse.lil_matrix или csr_matrix).

В частности, мне нужно будет вычислить псевдообратную матрицу для решения задачи наименьших квадратов. Я нашел метод sparse.lsqr, но он не очень эффективен. Есть ли способ для вычисления псевдообратного Мура-Пенроуза (соответствует pinv для нормальных матриц).

Размер моей матрицы A составляет около 600'000x2000, и в каждой строке матрицы у меня будет от 0 до 4 ненулевых значений. Размер матрицы A задается волокнистым пучком воксела х (волоконные тракты белого вещества), и мы ожидаем, что максимум 4 тракта пересекутся в вокселе. Мы ожидаем, что в большинстве вокселей белого вещества будет как минимум 1 тракт, но я скажу, что около 20% линий могут быть нулями.

Вектор b не должен быть разреженным, фактически b содержит меру для каждого вокселя, которая в общем случае не равна нулю.

Мне нужно минимизировать ошибку, но есть также некоторые условия на вектор х. Когда я попробовал модель на меньших матрицах, мне никогда не нужно было ограничивать систему, чтобы удовлетворить эти условия (в общем 0

Это поможет? Есть ли способ избежать принятия псевдообратного А?

Спасибо

Обновление 1 июня: Спасибо еще раз за помощь. Я ничего не могу показать вам о моих данных, потому что код на python доставляет мне некоторые проблемы. Однако, чтобы понять, как я мог выбрать хороший k, я попытался создать функцию тестирования в Matlab.

Код выглядит следующим образом:

F=zeros(100000,1000);

for k=1:150000
    p=rand(1);
    a=0;
    b=0;
    while a<=0 || b<=0
    a=random('Binomial',100000,p);
    b=random('Binomial',1000,p);
    end
    F(a,b)=rand(1);
end

solution=repmat([0.5,0.5,0.8,0.7,0.9,0.4,0.7,0.7,0.9,0.6],1,100);
size(solution)
solution=solution';
measure=F*solution;
%check=pinvF*measure;
k=250;
F=sparse(F);
[U,S,V]=svds(F,k);
s=svds(F,k);
plot(s)
max(max(U*S*V'-F))
for s=1:k
    if S(s,s)~=0
        S(s,s)=1/S(s,s);
    end
end

inv=V*S'*U';
inv*measure
max(inv*measure-solution)

У вас есть представление о том, что должно быть k по сравнению с размером F? Я взял 250 (более 1000), и результаты не являются удовлетворительными (время ожидания приемлемое, но не короткое). Также теперь я могу сравнить результаты с известным решением, но как вообще выбрать k? Я также приложил график 250 отдельных значений, которые я получаю, и их квадраты нормализованы. Я не знаю точно, как лучше сделать скриншот в матлабе. Теперь я продолжаю с большим k, чтобы увидеть, будет ли значение неожиданно намного меньше.

Еще раз спасибо, Дженнифер

Я ничего не знаю о реализации lsqr в sparse.linalg, но вряд ли найдется более быстрый способ вычисления псевдообратного; если так, то не будет ли lsqr просто вызывать это, а затем делать умножение? РЕДАКТИРОВАТЬ: Или проблема в том, что вам нужно сделать много таких вычислений и хотеть псевдообращения, чтобы вам не пришлось пересчитывать? bnaul

Ваш Ответ

2   ответа
6

scipy.sparse.linalg.

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

Вы можете описать немного более подробно вашу конкретную проблему (dot(A, x)= b+ e). По крайней мере, укажите:

«типичный» размерA«типичный» процент ненулевых записей вAнаименьших квадратов подразумевает, чтоnorm(e) свернут, но, пожалуйста, укажите, включен ли ваш основной интересx_hat или наb_hat, гдеe= b- b_hat а такжеb_hat= dot(A, x_hat)

Обновить: Если у вас есть представление о рангеA (и это намного меньше, чем количество столбцов), вы можете попробоватьнаименьших квадратов метод. Вот простая реализация, гдеk это число первых сингулярных значений и векторов, которые нужно использовать (то есть «эффективный» ранг).

from scipy.sparse import hstack
from scipy.sparse.linalg import svds

def tls(A, b, k= 6):
    """A tls solution of Ax= b, for sparse A."""
    u, s, v= svds(hstack([A, b]), k)
    return v[-1, :-1]/ -v[-1, -1]
@ Дженнифер: Во-вторых, есть некоторые теоретические причины, считайте, чтоsum(s)(из полногоsvd) объясняет все вариации в[A b]´, so you should find aгдеsum(s[:k]) объяснил бы разумное количество изменений в[A b]´ (consider this signal) andsum (s [k:]) `объяснит остальное (рассмотрим этот шум). Вы можете изучить линейную алгебру позадиsvd более подробно, но на практике вы могли бы построитьs[:k] для некоторого выполнимого k, и если вы в состоянии определить точку изменения там, используйте это какk, для более подробной информации Googlescree plot. eat
Есть ли лучшая возможность решить систему? Jennifer
@ Дженнифер: Во-первых, практические причины диктуют, что k должно быть относительно небольшим, обратите внимание, чтоu а такжеv скорее всего плотные. eat
@ Дженнифер: я обновила свой ответ. Существует реализация простого метода, который вы можете попробовать. Спасибо eat
Спасибо за ваш ответ, но в обоих случаях я не понимаю, как бороться с тем, что в SVDS k должно быть меньше ранга. Jennifer
4

я думаю, вы могли бы сделать это довольно легко, используяМур-Пенроуз СВД представительство, Найдите SVD с помощью scipy.sparse.linalg.svds, замените Sigma на псевдообратную, а затем умножьте V * Sigma_pi * U ', чтобы найти псевдообратную матрицу.

@eat вы правы, я переключил V и U. Я думаю, что все остальное правильно, но я согласен с точкой вашего ответа, что псевдообратный, вероятно, не лучший способ решения ее проблемы. bnaul

Похожие вопросы