是否可以为矩阵分解预先分配数组？

Question

我的问题不是F = svd(A)，可以先为 SVD 结构分配适当的内存，然后再做F .= svd(A)吗？

我的想法是这样的：

function main()
    F = Vector{SVD}(undef,10)
    # how to preallocate F?
    test(F)
end

function test(F::Vector{SVD})
    for i in 1:10
        F .= svd(rand(3,3))
    end
end

Answer 1

你的代码几乎可以工作。但你可能想要的是：

using LinearAlgebra

function main()
    F = Vector{SVD}(undef, 10)
    test(F)
end

function test(F::Vector{SVD})
    for i in 1:10
        F[i] = svd(rand(3, 3))
    end
    return F
end

你在for循环中的那一行是这样的：

F .= svd(rand(3,3))

它对每个循环执行相同的操作，因为您没有索引到F. 特别是，此操作试图将单个SVD对象广播到F循环的每次迭代中的所有字段中。（并且广播操作失败，因为默认情况下structs 被视为具有length方法的可迭代对象，但SVD没有length方法。）

但是，我建议不要在这种情况下预先分配向量。首先，让我们看一下 的类型F：

using LinearAlgebra

function main()
    F = Vector{SVD}(undef, 10)
    test(F)
end

function test(F::Vector{SVD})
    for i in 1:10
        F[i] = svd(rand(3, 3))
    end
    return F
end

这个向量的问题在于它是由抽象类型参数化的。手册的性能提示章节中有一节建议不要这样做。SVD是一个抽象类型，因为它的参数类型没有被指定。为了使其具体，您需要指定参数的类型，如下所示：

F .= svd(rand(3,3))

如您所见，当您使用像SVD. 此外，如果您这样做，您的代码将不会像它应有的那样通用。

解决此类问题的更好方法是使用映射、广播或列表理解。然后将自动推断正确的输出类型。这里有些例子：

列表理解

julia> typeof(Vector{SVD}(undef, 10))
Array{SVD,1}

地图

julia> SVD{Float64,Float64,Array{Float64,2}}
SVD{Float64,Float64,Array{Float64,2}}

julia> Vector{SVD{Float64,Float64,Array{Float64,2}}}(undef, 2)
2-element Array{SVD{Float64,Float64,Array{Float64,2}},1}:
 #undef
 #undef

广播

julia> [svd(rand(3, 3)) for _ in 1:2]
2-element Array{SVD{Float64,Float64,Array{Float64,2}},1}:
 SVD{Float64,Float64,Array{Float64,2}}([-0.6357040496635746 -0.2941425771794837 -0.7136949667270628; -0.45459999623274916 -0.6045700314848496 0.654090147040599; -0.6238743500629883 0.7402534845042064 0.2506104028424691], [1.4535849689665463, 0.7212190827260345, 0.05010669163393896], [-0.5975505057447164 -0.588792736048385 -0.5442945039782142; 0.7619724725128861 -0.6283345569895092 -0.15682358121595258; -0.2496624605679292 -0.5084474392397449 0.8241054891903787])
 SVD{Float64,Float64,Array{Float64,2}}([-0.5593632049776268 0.654338345992878 -0.5088753618327984; -0.6687620652652163 -0.7189576326033171 -0.18936003428293915; -0.4897653570633183 0.23439550227070827 0.8397551092645418], [1.8461274187259178, 0.21226179692488983, 0.14194607536315287], [-0.29089551972856004 -0.7086270946133293 -0.6428276887173754; -0.9203610429640889 0.023709029028269546 0.390350397126212; 0.2613720474647311 -0.7051847436823973 0.6590896221923739])

此外，映射、广播和列表推导应该与预分配向量一样有效。如果您正在做一个简单的映射，那么使用映射、广播或列表推导式通常更容易且更易读。预分配向量是我保留用于从头开始编写自定义算法的工具。

最后一点。在大多数情况下，类型参数被视为实现细节，而不是类型的公共 API 的一部分。因此，最好使用不依赖于固定类型参数类型的泛型编程方法。当然，这条经验法则也有一些例外，例如Array{T,N}和Dict{K,V}。