julia - 内存高效中心稀疏 SVD/PCA(在 Julia 中)？

Question

我有一个 300 万 x 900 万的稀疏矩阵，其中包含数十亿个非零条目。 R 和 Python 不允许包含超过 MAXINT 个非零条目的稀疏矩阵，因此我发现自己使用 Julia。

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

做 svd 的相关代码是 julia 可以在 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_Center 版本执行类似子类 AtA_or_AAT 的操作，该版本或多或少地模仿行为，但也存储列均值，并重新定义 A_mul_B!适当发挥作用。

但是，我不太使用 Julia，并且已经在修改内容时遇到了一些困难。在我再次尝试深入研究之前，我想知道我是否可以得到反馈，这是否被认为是一个合适的攻击计划，或者是否有一种更简单的方法在如此大的矩阵上进行 SVD(我还没有看过，但我可能错过了一些东西)。

编辑:我没有修改基本的 Julia，而是尝试写一个 "Centered Sparse Matrix"保持输入稀疏矩阵的稀疏结构的包，但在各种计算中适当的时候进入列意味着。它的实现受到限制，并且有效。不幸的是，它仍然太慢，尽管进行了相当广泛的努力来尝试优化事物。

Answer 1

在对稀疏矩阵算法进行了大量摸索之后，我意识到将乘法分布在减法上的效率要高得多:

如果我们的中心矩阵 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函数，只依赖于矩阵乘法和转置乘法。