内存有效集中稀疏SVD/PCA（在Julia中）？

我有一个300万×900万的稀疏矩阵，有几十亿个非零条目。 R和Python不允许超过MAXINT非零条目的稀疏矩阵，因此我发现自己使用了Julia。内存有效集中稀疏SVD/PCA（在Julia中）？

虽然用标准偏差缩放这些数据是微不足道的，但贬低当然是一种天真的方式，因为这会创建一个密集的200+太字节矩阵。

从我的阅读发现，该代码的一个关键要素是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的子类，其AtA_or_AAT_centered版本或多或少地模仿行为，但也存储列方法并重新定义A_mul_B！适当地运作。

但是，我不使用茱莉亚非常多，并已经遇到了一些困难，已经修改的东西。在我尝试再次探讨这个问题之前，我想知道如果我认为这是一个适当的攻击计划，或者如果在这样一个大矩阵上做一个简单的SVD方法，我能否得到反馈（我还没有看到它，但我可能错过了什么）。

编辑：我试过编写一个"Centered Sparse Matrix"包来保持输入稀疏矩阵的稀疏结构，而不是修改基本茱莉亚，而是在各种计算中进入适当的列。它的实施受到限制，并且工作成功。不幸的是，它仍然太慢，尽管有一些相当广泛的努力来尝试优化事情。

来源

2017-09-25 James

也许我误解了，但我认为你可以在scikit-learn中使用核心PCA吗？ –

我没有看到如何使用任何方法在scikit-learn中重新对齐数据。我错过了什么吗？ – James

与稀疏矩阵算法多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是致密的矢量矩阵的内积，以及1的“广播”（使用Julia的术语）到AX中间结果的每个行向量。大多数语言都有一种做这种广播的方式，而没有意识到额外的内存分配。

这就是我在我的package中为AcX和Ac'X实现的。然后可以将得到的对象传递给算法，如svds函数，它们只依赖于矩阵乘法和转置乘法。

来源

2017-10-04 22:49:36 James

内存有效集中稀疏SVD/PCA（在Julia中）？

回答

相关问题