Estoy trabajando con datos de neuroimagen y, debido a la gran cantidad de datos, me gustaría utilizar matrices dispersas para mi código (scipy.sparse.lil_matrix o csr_matrix).

En particular, necesitaré calcular el pseudoinverso de mi matriz para resolver un problema de mínimos cuadrados. He encontrado el método sparse.lsqr, pero no es muy eficiente. ¿Existe algún método para calcular el pseudoinverso de Moore-Penrose (correspondiente a pinv para matrices normales)?

El tamaño de mi matriz A es de aproximadamente 600'000x2000 y en cada fila de la matriz tendré de 0 a 4 valores distintos de cero. El tamaño de la matriz A viene dado por el paquete de fibras voxel x (tractos de fibra de materia blanca) y esperamos que se crucen un máximo de 4 tractos en un voxel. En la mayoría de los vóxeles de la materia blanca esperamos tener al menos 1 tracto, pero diré que alrededor del 20% de las líneas podrían ser ceros.

El vector b no debe ser escaso, en realidad b contiene la medida para cada vóxel, que en general no es cero.

Necesitaría minimizar el error, pero también hay algunas condiciones en el vector x. Cuando probé el modelo en matrices más pequeñas, nunca tuve que restringir el sistema para satisfacer estas condiciones (en general 0

¿Es eso de alguna ayuda? ¿Hay alguna manera de evitar tomar el pseudoinverso de A?

Gracia

Update 1 de junio: Gracias de nuevo por la ayuda. Realmente no puedo mostrarle nada sobre mis datos, porque el código en Python me da algunos problemas. Sin embargo, para entender cómo podría elegir una buena k, intenté crear una función de prueba en Matlab.

El código es el siguiente:

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)

¿Tienes alguna idea de lo que debería ser k en comparación con el tamaño de F? He tomado 250 (más de 1000) y los resultados no son satisfactorios (el tiempo de espera es aceptable, pero no corto). También ahora puedo comparar los resultados con la solución conocida, pero ¿cómo podría uno elegir k en general? También adjunté la gráfica de los 250 valores individuales que obtengo y sus cuadrados se normalizaron. No sé exactamente cómo hacer mejor un diagrama de pantalla en matlab. Ahora procedo con una k más grande para ver si de repente el valor será mucho más pequeño.

Gracias de nuevo, Jennifer