内存高效中心稀疏SVD/PCA(朱莉娅)？

Question

我有一个300万x 900万的稀疏矩阵,有数十亿个非零项.R和Python不允许稀疏矩阵具有超过MAXINT非零的条目,因此我发现自己使用Julia.

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

有关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_centered版本的子类AtA_or_AAT,或多或少模仿行为,但也存储列方式,并重新定义A_mul_B!功能恰当.

但是,我并没有非常使用朱莉娅,并且已经遇到了一些难以修改的问题.在我再次尝试深入研究之前,我想知道如果这被认为是一个合适的攻击计划,或者如果在这么大的矩阵上进行SVD​​的简单方法,我是否能得到反馈(我没有看到它,但我可能错过了一些东西).

编辑:我没有修改基础Julia,而是尝试编写一个"中心稀疏矩阵"包来保持输入稀疏矩阵的稀疏结构,但是在各种计算中适当地进入列.它的实施范围有限,而且有效.不幸的是,它仍然太慢,尽管为了优化事物做了一些非常广泛的努力.

Answer 1

经过对稀疏矩阵算法的大量研究后，我意识到将乘法分布在减法上的效率要高得多：

如果我们的中心矩阵Ac是由原始nxm矩阵A及其列向量组成的M，其中nx1向量我将调用1。我们乘以mxk矩阵X

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

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

AX可以与通常的稀疏矩阵乘法函数一起执行，M'X是一个稠密向量矩阵内积，并将 1 的向量“广播”（使用 Julia 的术语）到矩阵的每一行AX。大多数语言都有一种在不实现额外内存分配的情况下进行广播的方法。

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