메모리 효율이 좋은 스파 스 SVD/PCA (줄리아에서)?

Question

수십억 개의 0이 아닌 항목이 포함 된 3 백만 x 900만의 희소 행렬이 있습니다. R과 Python은 MAXINT가 아닌 항목이 0 인 희소 행렬을 허용하지 않기 때문에 Julia를 사용하는 이유를 발견하게되었습니다.

이 데이터를 표준 편차로 스케일링하는 것은 쉬운 일이 아니지만, 비열한 것은 물론 200+ 테라 바이트 밀도의 고밀도 매트릭스를 생성 할 수있는 순진한 방식으로 진행되지 않습니다.

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에 대한 관련 구성 요소에서 처리 아래 복사.

나는이 작업을 'flying'타입 셋업에 사용하는 가장 좋은 방법은 AtA_or_AAT_centered 버전을 사용하여 서브 클래스 AtA_or_AAT와 같은 일을하는 것이라고 생각한다. 그러나 동작을 모방하거나 컬럼을 저장하고 다시 정의한다. A_mul_B! 적절하게 기능하십시오.

그러나 저는 Julia를 사용하지 않고 이미 문제를 수정하는 데 어려움을 겪고 있습니다. 이것을 다시 시도하기 전에 이것이 적절한 공격 계획으로 간주된다면 의견을 얻을 수 있을지, 아니면 그러한 대규모 매트릭스에서 SVD를 수행하는 훨씬 쉬운 방법이 있다면 궁금합니다. 그것을 보았다. 그러나 나는 무엇인가 놓쳤을지도 모른다).

편집 : 기본 줄리아를 수정하는 대신 입력 희소 매트릭스의 희박 구조를 유지하는 "Centered Sparse Matrix" 패키지를 작성하려고 시도했지만 다양한 계산에서 적절한 경우 해당 열을 입력합니다. 그것은 구현 된 것으로 제한되어 있으며 작동합니다. 불행히도, 물건을 최적화하려고 시도한 꽤 광범위한 노력에도 불구하고 여전히 너무 느립니다.

Answer 1

는 스파 스 매트릭스 알고리즘을 많이 fuddling 후, 나는 뺄셈보다 곱셈을 배포하는 것은 극적으로보다 효율적인 것을 깨달았다 :

우리의 중심 매트릭스 Ac은 원래 n X m 매트릭스 A와 자사의 벡터로 형성되는 경우 열은 M을 의미하고 n x 1 벡터는 단지 1이라고 부릅니다. 우리는 m X k 매트릭스 X

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

곱이며 우리는 기본적으로 수행됩니다. 어리석은 간단한, 실제로.

AX 통상 희소 행렬 곱셈 기능을 수행 할 수있다, M'X 치밀한 벡터 행렬 내적 및 AX 중간 결과의 각 행에 하나의 '방송'(줄리아의 용어를 사용하는)의 벡터이다 . 대부분의 언어는 추가 메모리 할당을 인식하지 않고 해당 브로드 캐스팅을 수행 할 수있는 방법이 있습니다.

이것은 내 package의 AcX와 Ac'X에서 구현 한 것입니다. 그런 다음 결과 객체를 svds 함수와 같은 알고리즘에 전달할 수 있습니다.이 함수는 행렬 곱셈 및 곱셈 곱하기에만 의존합니다.