内存高效中心稀疏 SVD/PCA（在 Julia 中）？答案

【问题标题】：Memory Efficient Centered Sparse SVD/PCA (in Julia)?内存高效中心稀疏 SVD/PCA（在 Julia 中）？
【发布时间】：2018-03-06 01:35:03
【问题描述】：

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

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

根据我的阅读，这段代码的一个关键元素是 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！功能正常。

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

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

【问题讨论】：

也许我理解错了，但我认为你可以在 scikit-learn 中做核心 PCA 之外的事情？
我没有看到如何使用任何方法将 scikit-learn 中的数据重新居中。我错过了什么吗？

标签： julia sparse-matrix pca svd

【解决方案1】：

在大量使用稀疏矩阵算法后，我意识到将乘法分配到减法上会大大提高效率：

如果我们的中心矩阵Ac 是由原始nxm 矩阵A 和它的列向量表示M 形成的，则nx1 向量是我只会打电话给1。我们乘以 mxk 矩阵 X

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

我们基本上完成了。实际上非常简单。

AX 可以用通常的稀疏矩阵乘法函数进行，M'X 是一个密集的向量矩阵内积，而 1 的“广播”（使用 Julia 的术语）的向量到每一行AX 中间结果。大多数语言都有一种广播方式，无需意识到额外的内存分配。

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

【讨论】：