内存高效集中稀疏 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 的相关组件中。

我认为使 'centering on the fly' 类型设置工作的最佳方法是使用 AtA_or_AAT_centered 版本做一些类似子类 AtA_or_AAT 的事情，该版本或多或少模仿了行为而且还存储列均值，并重新定义 A_mul_B!正常运作。

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

编辑：我没有修改基础 Julia，而是尝试编写一个 "Centered Sparse Matrix" 包来保持输入稀疏矩阵的稀疏结构，但在各种计算中适当地输入列均值。它在实施方面受到限制，并且有效。不幸的是，它仍然太慢了，尽管付出了相当大的努力来尝试优化事物。

Answer 1

在大量使用稀疏矩阵算法之后，我意识到将乘法分布在减法上的效率要高得多：

如果我们的居中矩阵 Ac 是由原始 nxm 矩阵 A 及其列向量表示 M 形成的，具有nx1 个向量，我将称之为 1。我们乘以 mxk 矩阵 X

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

我们基本上完成了。其实很简单。

AX可以用通常的稀疏矩阵乘法函数进行，M'X是稠密向量-矩阵内积，1的向量"broadcasts"（用Julia 的术语）到 AX 中间结果的每一行。大多数语言都有一种方法可以在不实现额外内存分配的情况下进行广播。

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

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

Memory Efficient Centered Sparse SVD/PCA (in Julia)?

julia

sparse-matrix

pca

svd