SVD/PCA épars centrés sur la mémoire (dans Julia)?

Question

J'ai une matrice clairsemée de 3 millions x 9 millions avec plusieurs milliards d'entrées non nulles. R et Python n'autorisent pas les matrices clairsemées avec plus de MAXINT entrées non nulles, donc pourquoi je me suis retrouvé en utilisant Julia. Bien que la mise à l'échelle de ces données avec l'écart-type soit triviale, l'avilissement est bien entendu une façon naïve de créer une matrice dense de 200+ téraoctets.

Le code pertinent pour faire SVD est julia se trouve à https://github.com/JuliaLang/julia/blob/343b7f56fcc84b20cd1a9566fd548130bb883505/base/linalg/arnoldi.jl#L398

De ma lecture, un élément clé de ce code est le struct AtA_or_AAt et plusieurs des fonctions autour de celles-ci, en particulier A_mul_B !. Copié ci-dessous pour votre commodité

struct AtA_or_AAt{T,S} <: AbstractArray{T, 2} 
    A::S 
    buffer::Vector{T} 
end 

function AtA_or_AAt(A::AbstractMatrix{T}) where T 
    Tnew = typeof(zero(T)/sqrt(one(T))) 
    Anew = convert(AbstractMatrix{Tnew}, A) 
    AtA_or_AAt{Tnew,typeof(Anew)}(Anew, Vector{Tnew}(max(size(A)...))) 
end 

function A_mul_B!(y::StridedVector{T}, A::AtA_or_AAt{T}, x::StridedVector{T}) where T 
    if size(A.A, 1) >= size(A.A, 2) 
     A_mul_B!(A.buffer, A.A, x) 
     return Ac_mul_B!(y, A.A, A.buffer) 
    else 
     Ac_mul_B!(A.buffer, A.A, x) 
     return A_mul_B!(y, A.A, A.buffer) 
    end 
end 
size(A::AtA_or_AAt) = ntuple(i -> min(size(A.A)...), Val(2)) 
ishermitian(s::AtA_or_AAt) = true 

Ceci est passé dans la fonction GIEs, où la magie se produit, et la sortie est ensuite traitée pour les composants pertinents pour SVD. Je pense que la meilleure façon de faire cela pour une installation de type 'centrer à la volée' est de faire quelque chose comme la sous-classe AtA_ou_AAT avec une version AtA_ou_AAT_centered qui imite plus ou moins le comportement mais stocke aussi les moyens de colonne, et redéfinit le A_mul_B! fonctionner correctement. Cependant, je n'utilise pas beaucoup Julia et j'ai déjà eu du mal à modifier certaines choses. Avant d'essayer de replonger dans le sujet, je me demandais si je pouvais avoir un retour si cela serait considéré comme un plan d'attaque approprié, ou s'il y a simplement une façon beaucoup plus facile de faire de la SVD sur une telle matrice (je n'ai pas vu, mais j'ai peut-être manqué quelque chose). Edit: Plutôt que de modifier la base de Julia, j'ai essayé d'écrire un paquet "Centered Sparse Matrix" qui conserve la structure éparse de la matrice d'entrée clairsemée, mais qui entre dans les moyens de colonne le cas échéant dans divers calculs. C'est limité dans ce qu'il a mis en place, et ça fonctionne. Malheureusement, c'est encore trop lent, malgré quelques efforts assez importants pour essayer d'optimiser les choses.

Answer 1

Après beaucoup fuddling avec l'algorithme de matrice creuse, je réalise que la distribution de la multiplication sur la soustraction était considérablement plus efficace:

Si notre matrice centrée Ac est formée à partir de la matrice n x m originale A et son vecteur de colonne signifie M, avec un n x 1 vecteur de ceux que je vais simplement appeler 1. Nous multiplions par une matrice m x kX

Ac := (A - 1M') 
AcX = X 
    = AX - 1M'X 

Et nous sommes fondamentalement fait. Stupidement simple, en fait.

AX est peut être effectuée avec la fonction de multiplication de la matrice creuse habituelle, M'X est un vecteur-matrice dense produit scalaire et le vecteur des « émissions » de 1 (pour reprendre la terminologie de Julia) à chaque ligne du résultat intermédiaire AX . La plupart des langues ont un moyen de faire cette diffusion sans réaliser l'allocation de mémoire supplémentaire.

C'est ce que j'ai implémenté dans mon package pour AcX et Ac'X. L'objet résultant peut ensuite être transmis à des algorithmes, tels que la fonction svds, qui dépendent uniquement de la multiplication matricielle et de la multiplication transposée.