メモリ効率のよいスパースSVD/PCA（ジュリア）？

Question

私は、数十億の非ゼロエントリを持つ300万x 900万の疎な行列を持っています。 RとPythonは、MAXINT以外のエントリを持つスパース行列を許可していないため、Juliaを使って自分自身を見つけたのです。

このデータを標準偏差でスケーリングするのは簡単ではありませんが、デメメインングはもちろん、200MB以上の高密度のマトリックスを作成するという単純な方法では無駄です。

SVDを行うための関連するコードはジュリアが私の読書からhttps://github.com/JuliaLang/julia/blob/343b7f56fcc84b20cd1a9566fd548130bb883505/base/linalg/arnoldi.jl#L398

で見つけることができますされ、このコードの重要な要素は、AtA_or_AAt構造体と、それらの周りの機能のいくつかは、具体的にA_mul_Bです！便利

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 

のために以下にコピーされたこれは、いくつかのマジックが起こる関数eigs機能、に渡され、出力はSVDに関連する部品に加工されます。

私はこの設定を「オンザフライ・オン・ザ・フライ」タイプの設定にする最良の方法はAtA_or_AAT_centeredバージョンのサブクラスAtA_or_AATのようなことをすることだと思いますが、多かれ少なかれ動作を模倣するだけでなく、 A_mul_B！適切に機能する。

しかし、私はジュリアをあまり使用していないし、すでにいくつかの困難に挑戦しています。私はこれをもう一度やってみる前に、これが適切な攻撃計画と考えられるならフィードバックを得ることができるのか、そういう大きなマトリックスでSVDをやる方が簡単な方法があればそれを見ましたが、私は何かを見逃しているかもしれません）。

編集：基本的なJuliaを変更するのではなく、入力スパース行列のスパース性構造を維持するが、さまざまな計算で適切な場所に列方向に入る"Centered Sparse Matrix"パッケージを作成しようとしました。それは実装されているものに限られており、動作します。残念ながら、物事を最適化しようとするかなりの努力にもかかわらず、まだ遅すぎます。

Answer 1

は疎行列アルゴリズムとはるかにfuddling後、私は、減算上乗算を分配する劇的に、より効率的であったことに気づいた：

我々の中心マトリックスAc元n X mマトリックスAとそのベクトルから形成されている場合列はMを意味し、nx1のベクターを使用します。1と呼ぶことにします。私たちは、m X k行列X

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

乗じあり、私たちは基本的に行われます。愚かな単純な、実際には。

AX通常疎行列乗算機能を行うことができるされ、M'Xは、密ベクトル、行列の内積、及びAX中間結果の各列に1の「ブロードキャスト」（ジュリアの用語を使用する）のベクトルであります。ほとんどの言語は余分なメモリ割り当てを実現することなくその放送を行う方法を持っています。

これは、私がpackageのAcXとAc'Xで実装したものです。結果として得られるオブジェクトは、svds関数のようなアルゴリズムに渡すことができます。この関数は、行列乗算と転置乗算のみに依存します。