0

以下代码行出现在优化包“Optim”中的 levenberg-marquardt 算法中:

DtD = diagm(Float64[max(x, MIN_DIAGONAL) for x in sum(J.^2,1)])
delta_x = ( J'*J + sqrt(lambda)*DtD ) \ -J'*fcur

但是,我的问题与算法或包的任何特定内容无关。我想这更多地与基础 julia 中的线性代数和因式分解有关。

如果我有一个完整的矩阵 J,则以下工作:

In  [3]: n = 5
J = rand(n,n)
fcur = rand(n)
DtD = diagm(Float64[max(x, 1e-6) for x in sum(J.^2,1)])
( J'*J + sqrt(100)*DtD ) \ -J'*fcur

Out [3]: 5-element Array{Float64,1}:
 -0.0740316
 -0.0643279
 -0.0665033
 -0.10568  
 -0.0599613

但是,如果 J 是稀疏的,我会收到一个错误:

In  [4]: J = sparse(J)
DtD = diagm(Float64[max(x, 1e-6) for x in sum(J.^2,1)])
( J'*J + sqrt(100)*DtD ) \ -J'*fcur

Out [4]: ERROR: `A_ldiv_B!` has no method matching A_ldiv_B!(::Cholesky{Float64}, ::SparseMatrixCSC{Float64,Int64})
while loading In[4], in expression starting on line 3

 in \ at linalg/generic.jl:233

据我所知(作为初学者,我对 julia 的了解有限),发生此错误是因为 julia 尝试先计算( J'*J + sqrt(100)*DtD ) \ -J'我的第一个问题是我如何知道 julia 在实现上述代码时所采用的路径?我知道@which但我不知道如何使用它来到达 A_ldiv_B!因为这应该从\(A,B)A_ldiv_B 开始,然后以某种方式结束!:

In  [6]: a = ( J'*J + sqrt(100)*DtD ); b = -J'; @which a\b

Out [6]: \(A::Union(SubArray{T,2,A<:DenseArray{T,N},I<:(Union(Range{Int64},Int64)...,)},DenseArray{T,2}),B::Union(SubArray{T,2,A<:DenseArray{T,N},I<:(Union(Range{Int64},Int64)...,)},SubArray{T,1,A<:DenseArray{T,N},I<:(Union(Range{Int64},Int64)...,)},DenseArray{T,2},DenseArray{T,1})) at linalg/dense.jl:409

另请注意:

In  [7]: typeof(a)

Out [7]: Array{Float64,2}

In  [8]: typeof(b)

Out [8]: Array{Float64,2}

这使得这更加令人困惑,因为上面没有 Cholesky 类型。我的第二个问题是:Cholesky 类型是如何出现的?错误消息说:A_ldiv_B!没有方法匹配 A_ldiv_B!(:: Cholesky{Float64} , ::SparseMatrixCSC{Float64,Int64})

我偶然发现的另一个有趣的点是,如果稀疏矩阵的大小为 (2,2),则不会发生上述错误:

In  [9]: n = 2
J = sparse(rand(n,n))
fcur = rand(n)
DtD = diagm(Float64[max(x, 1e-6) for x in sum(J.^2,1)])
( J'*J + sqrt(100)*DtD ) \ -J'*fcur

Out [9]: 2-element Array{Float64,1}:
 -0.0397989
 -0.052129 

最后,我可以通过插入括号来解决这个问题-J'*fcur,这似乎是作者的意图。但我很困惑。任何想法都非常感谢。谢谢!

In  [12]: n = 5
J = sparse(rand(n,n))
fcur = rand(n)
DtD = diagm(Float64[max(x, 1e-6) for x in sum(J.^2,1)])
( J'*J + sqrt(100)*DtD ) \ (-J'*fcur)

Out [12]: 5-element Array{Float64,1}:
 -0.0536266
 -0.0692286
 -0.0673166
 -0.0616349
 -0.0559813
4

1 回答 1

1

当您遇到在解析期间使用替换的代码时,准确地确定采用哪个代码路径可能有点棘手,例如'. 您可以尝试 julia> ( J'*J + sqrt(100)*DtD ) \ -J'fcur 查看另一个替换发生。

我不知道这是否真的是您问题的答案,但我会注意关于您的示例的三件事。

  1. Julia 从左到右解析,所以我认为这个例子应该像 (( J'*J + sqrt(100)*DtD ) \ ctranspose(-J))*fcur
  2. 我们还没有在求解中实现稀疏右手边,因为即使binAx=b是稀疏的,通常结果也不是稀疏的,因此利用稀疏性的收益b是适度的。

  3. 解决问题的“正确”方法是在周围加上括号,(-Jfcur)因为这样解决方案是稀疏矩阵向量乘法和稀疏矩阵向量求解,而不是稀疏矩阵矩阵求解和密集矩阵向量乘法。

于 2015-03-11T21:04:19.813 回答