A_ldiv_B! с разреженными матрицами

Question

Следующие строчки кода появляются в алгоритме Левенберга-Марквардта в пакете оптимизации "Оптим":

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

Итак, насколько я понимаю (с моим ограниченным знанием Юлии как новичка), эта ошибка происходит, потому что Юлия пытается вычислить ( J'*J + sqrt(100)*DtD ) \ -J' первый. Мой первый вопрос: как я могу узнать, какой путь выберет Юлия при реализации приведенного выше кода? Я в курсе @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}

Что делает это еще более запутанным, так как в приведенном выше нет никакого типа Холецкого. Мой второй вопрос: как появляется тип Холецкого? Сообщение об ошибке говорит: 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

Answer 1

Может быть немного сложно выяснить, какой именно путь к коду взят, когда вы работаете с кодом, который использует подстановки во время синтаксического анализа, как в случае с ', Вы могли бы попробовать julia> ( J'*J + sqrt(100)*DtD ) \ -J'fcur чтобы увидеть другую замену.

Я не знаю, действительно ли это ответ на ваш вопрос, но я отмечу три вещи в вашем примере.

1. Юля разбирает слева направо, поэтому я думаю, что пример должен читаться как (( J'*J + sqrt(100)*DtD ) \ ctranspose(-J))*fcur

2. Мы не реализовали разреженную правую часть в решениях, потому что даже когда b в Ax=b как правило, результат не является разреженным, поэтому выгода от использования разреженности b скромный.

3. "Правильный" способ решения проблемы - добавить круглые скобки (-Jfcur) потому что тогда решением является разреженное матричное-векторное умножение и разреженное матричное-векторное решение вместо разреженного матрично-матричного решения и плотного матрично-векторного умножения.