2nd Order Newton's Method with Sparse Hessian

General Usage
7

Thanks for the suggestions! I ended up going with MUMPS.jl and adding a scalar multiple of the identity. What pushed me towards this package: I can do one symbolic factorization, then repeated numeric factorization/solve calls. The following should give the reader the general idea. (The full code is PR #148 in PowerFlows.jl.)

using MUMPS, MPI # may need MPIPreferences too.

@enum MUMPS_JOB begin
    INIT = -1
    CLEANUP = -2
    ANALYZE = 1
    FACTOR = 2
    SOLVE = 3
    SAVE = 7
    RESTORE = 8
    DELETE = -3
    FACTOR_CLEANUP = -4
end

function mumps_job!(mumps::Mumps, job::MUMPS_JOB)
    MUMPS.set_job!(mumps, Integer(job))
    MUMPS.invoke_mumps!(mumps)
    return
end

const β = 10.0^-3
"""Modify the Hessian to be positive definite, to guarantee that the solution 
to `H Δx = -∇f` is a direction of descent. Currently using algorithm 3.3 from 
Nocedal & Wright: add a multiple of the identity."""
function modify_hessian!(H::SparseMatrixCSC, mumps::Mumps)
    minDiagElem = minimum(H[i, i] for i in axes(H, 1))
    if minDiagElem > 0.0
        τ = 0.0
    else
        τ = -minDiagElem + β
        for i in axes(H, 1)
            H[i, i] += τ
        end
    end
    MUMPS.associate_matrix!(mumps, H)
    mumps_job!(mumps, FACTOR)
    while mumps.infog[12] > 0 # while matrix isn't positive definite.
        mumps_job!(mumps, FACTOR_CLEANUP)
        τ_old = τ
        τ = max(2 * τ, β)
        for i in axes(H, 1)
            H[i, i] += τ - τ_old # now try H + τ*I
        end
        MUMPS.associate_matrix!(mumps, H)
        mumps_job!(mumps, FACTOR)
    end
end

# setup code: run once
if !MPI.Initialized()
    MPI.Init()
end
icntl = deepcopy(MUMPS.default_icntl)
icntl[4] = 1 # report errors only
icntl[13] = 1 # due to our _modify_hessian! strategy,
# we need to know the exact number of negative pivots.
mumps = Mumps{Float64}(MUMPS.mumps_symmetric, icntl, MUMPS.default_cntl32)
MPI.add_finalize_hook!(() -> MUMPS.finalize(mumps))
# [initialize sparse structure of H here]
MUMPS.associate_matrix!(mumps, H)
mumps_job!(mumps, ANALYZE) # symbolic factorization.

# run each iteration.
# [compute H and grad at x here]
modify_hessian!(H, mumps)
MUMPS.associate_rhs!(mumps, grad)
mumps_job!(mumps, SOLVE)
MUMPS.mumps_solve!(Δx, mumps)
Δx .*= -1
# retain the symbolic factorization, so we can FACTOR and SOLVE next iteration.
mumps_job!(mumps, FACTOR_CLEANUP)
# [do a line search or otherwise decide on step size c]
x += c * Δx

It’d be lovely if one of the various packages integrated something like this: a 2nd order newton solver, reusing the symbolic factorization of the sparse hessian.

edit: well, it turns out the default MPI backend for windows is currently broken. Relevant issue

1 Like