A_ldiv_B！稀疏矩阵

Question

在优化包"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

但是，我的问题与算法或任何特定于程序包的内容无关。我想这更多地与 base 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! has no method matching 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

Answer 1

当您运行进入在解析过程中使用替换的代码时（如 ' 的情况），要准确找出采用的代码路径可能有点棘手。你可以试试 julia> ( J'*J + sqrt(100)*DtD ) \ -J'fcur 看到另一个替换发生。

我不知道这是否真的是你问题的答案，但我会注意到关于你的例子的三点。

Julia 从左到右解析，所以我认为这个例子应该这样读 (( J'*J + sqrt(100)*DtD ) \ ctranspose(-J))*fcur
我们没有在求解中实现稀疏右侧，因为即使 Ax=b 中的 b 是稀疏的，通常结果也不是稀疏的，所以利用b 的稀疏度适中。
"right" 解决问题的方法是在 (-Jfcur) 周围添加括号，因为这样解决方案就是稀疏矩阵向量乘法和稀疏矩阵向量求解而不是稀疏矩阵-矩阵求解和密集矩阵-向量乘法。

A_ldiv_B！稀疏矩阵

A_ldiv_B! with sparse matrices

optimization

linear-algebra

julia

factorization

levenberg-marquardt