以下代码行出现在优化包"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
但是，我的问题与算法或特定于包的任何内容无关。我想这更多地与基朱莉娅中的线性代数和因式分解有关。
如果我有一个完整的矩阵 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 ab
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类型。我的第二个问题是：乔列斯基类型是如何出现的？错误消息显示：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
当您遇到在解析过程中使用替换的代码时，准确确定采用的代码路径可能有点棘手，就像'一样。你可以试试
julia> ( J'*J + sqrt(100)*DtD ) -J'fcur
看到另一个替换发生。
我不知道这是否真的是你问题的答案，但我会注意到关于你的例子的三件事。
- 朱莉娅从左到右解析，所以我认为这个例子应该读成
(( J'*J + sqrt(100)*DtD ) ctranspose(-J))*fcur -
我们还没有在求解中实现稀疏右侧，因为即使
Ax=b中的b是稀疏的，那么通常结果也不是稀疏的，因此利用b稀疏的收益是适度的。 -
解决问题的"正确"方法是在
(-Jfcur)周围添加括号，因为这样解决方案是稀疏矩阵向量乘法和稀疏矩阵向量求解，而不是稀疏矩阵矩阵求解和密集矩阵向量乘法。