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

отаю с данными из нейровизуализации, и из-за большого объема данных я хотел бы использовать разреженные матрицы для моего кода (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, чтобы увидеть, будет ли значение неожиданно намного меньше.

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

Ответы на вопрос(2)

Ваш ответ на вопрос