Add incomplete LU preconditioner

[MaximilianJHuber]

I have had a good experience using an incomplete LU as a preconditioner for some linear high-dimensional fixed-effect models. I tried to feed it to LeastSquaresOptim, but seem to fail. What am I doing wrong?

A sparse, not diagonally dominant MWE with 1000 variables and observations:

using LeastSquaresOptim, IncompleteLU, SparsityDetection, SparseArrays, SparseDiffTools, BenchmarkTools

f! = function(y,x)
    for i in 2:length(x)-1
        y[i] = x[i-1] - 2x[i] + x[i+1]
    end
    y[1] = x[1] - 10
    y[end] = x[end-1] - 2x[end]
nothing
end

x = rand(1000)
y = similar(x)

# detect sparsity and color the Jacobian
sparsity_pattern = jacobian_sparsity(f!, y, x)
jac = Float64.(sparse(sparsity_pattern))
colors = matrix_colors(jac)
g! = (jac, x) -> forwarddiff_color_jacobian!(jac, f!, x, colorvec = colors)

result = optimize!(
    LeastSquaresProblem(deepcopy(x), y, f!, jac, g!), 
    LevenbergMarquardt(LeastSquaresOptim.LSMR()),
    show_trace = false)

result.iterations, result.ssr

yields (1000, 1.0825449074177966e-5). It works, but is slow.

I thought I could inject the LU via:

preconditioner! = function(P, x, J, λ)
    P = ilu(J'J .+ diagm(λ), τ = 0.001).U
end

import LinearAlgebra.ldiv!
function ldiv!(Y::Vector{Float64}, A::SparseMatrixCSC{Float64, Int64}, B::Vector{Float64})
    Y .= A \ B
end

result = optimize!(
    LeastSquaresProblem(deepcopy(x), y, f!, jac, g!), 
    LevenbergMarquardt(LeastSquaresOptim.LSMR(preconditioner!, ilu(jac'jac, τ = 0.001).U)),
    show_trace = false)

result.iterations, result.ssr

but it yields (1000, 1.3643893415608515e6).