Skip to content

Commit

Permalink
General tweaks
Browse files Browse the repository at this point in the history
  • Loading branch information
@ojwoodford
ojwoodford authored Oct 26, 2023
1 parent d4e3d98 commit bfa26ce
Show file tree
Hide file tree
Showing 9 changed files with 51 additions and 31 deletions.
2 changes: 1 addition & 1 deletion examples/rosenbrock.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ NLLSsolver.ndeps(::Rosenbrock) = static(1) # Residual depends on 1 variable
NLLSsolver.nres(::Rosenbrock) = static(2) # Residual has length 2
NLLSsolver.varindices(::Rosenbrock) = SVector(1) # There's only one variable
NLLSsolver.getvars(::Rosenbrock, vars::Vector) = (vars[1]::NLLSsolver.EuclideanVector{2, Float64},)
NLLSsolver.computeresidual(res::Rosenbrock, x) = SVector(res.a - x[1], res.b * (x[1] ^ 2 - x[2]))
NLLSsolver.computeresidual(res::Rosenbrock, x) = SVector(res.a * (1 - x[1]), res.b * (x[1] ^ 2 - x[2]))

function constructrosenbrockprob()
# Create the problem
Expand Down
12 changes: 8 additions & 4 deletions src/callbacks.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,21 +5,25 @@ nullcallback(cost, unusedargs...) = (cost, 0)

# Print out per-iteration results
function printoutcallback(cost, problem, data, trailingargs...)
prevcost = data.bestcost
if data.iternum == 1
prevcost = data.startcost
# First iteration, so print out column headers and the zeroth iteration (i.e. start) values
println("iter cost cost_change |step|")
@printf("% 4d % 8e % 4.3e % 3.2e\n", 0, data.bestcost, 0, 0)
@printf("% 4d % 8e % 4.3e % 3.2e\n", 0, prevcost, 0, 0)
end
@printf("% 4d % 8e % 4.3e % 3.2e\n", data.iternum, cost, data.bestcost-cost, norm(data.linsystem.x))
@printf("% 4d % 8e % 4.3e % 3.2e\n", data.iternum, cost, prevcost-cost, norm(data.linsystem.x))
return cost, 0
end
function printoutcallback(cost, data, trradius::Float64)
prevcost = data.bestcost
if data.iternum == 1
prevcost = data.startcost
# First iteration, so print out column headers and the zeroth iteration (i.e. start) values
println("iter cost cost_change |step| tr_radius")
@printf("% 4d % 8e % 4.3e % 3.2e % 2.1e\n", 0, data.bestcost, 0, 0, trradius)
@printf("% 4d % 8e % 4.3e % 3.2e % 2.1e\n", 0, prevcost, 0, 0, trradius)
end
@printf("% 4d % 8e % 4.3e % 3.2e % 2.1e\n", data.iternum, cost, data.bestcost-cost, norm(data.linsystem.x), trradius)
@printf("% 4d % 8e % 4.3e % 3.2e % 2.1e\n", data.iternum, cost, prevcost-cost, norm(data.linsystem.x), trradius)
return cost, 0
end

Expand Down
28 changes: 20 additions & 8 deletions src/iterators.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,9 @@ using SparseArrays

negate!(x) = @.(x = -x)

# Default preoptimization - do nothing, return lowest cost possible
preoptimization(::Any, unusedargs...) = -Inf

# Iterators assume that the linear problem has been constructed

# Newton optimization (undamped-Hessian form)
Expand Down Expand Up @@ -32,15 +35,14 @@ mutable struct DoglegData{T}
return new{typeof(data.linsystem.x)}(0.0, similar(data.linsystem.x))
end
end
gettr(dd::DoglegData) = dd.trustradius
settr!(dd::DoglegData, tr) = dd.trustradius = tr
function reset!(dd::DoglegData{T}, ::NLLSProblem, data::NLLSInternal) where T<:Vector
dd.trustradius = 0.0
settr!(dd, 0.0)
resize!(dd.cauchy, length(data.linsystem.x))
return
end
function reset!(dd::DoglegData{T}, ::NLLSProblem, data::NLLSInternal) where T<:StaticVector
dd.trustradius = 0.0
return
end
reset!(dd::DoglegData{T}, ::NLLSProblem, data::NLLSInternal) where T<:StaticVector = settr!(dd, 0.0)

function iterate!(doglegdata::DoglegData, data, problem::NLLSProblem, options::NLLSOptions)::Float64
hessian, gradient = gethessgrad(data.linsystem)
Expand Down Expand Up @@ -122,13 +124,23 @@ mutable struct LevMarData
return new(0.0)
end
end
reset!(lmd::LevMarData, ::NLLSProblem, ::NLLSInternal) = lmd.lambda = 0.0
gettr(lmd::LevMarData) = lmd.lambda
settr!(lmd::LevMarData, tr) = lmd.lambda = tr
reset!(lmd::LevMarData, ::NLLSProblem, ::NLLSInternal) = settr!(lmd, 0.0)

function initlambda(hessian)
m = zero(eltype(hessian))
for i in indices(hessian, 1)
@inbounds m = max(m, abs(hessian[i,i]))
end
return m * 1e-6
end

function iterate!(levmardata::LevMarData, data, problem::NLLSProblem, options::NLLSOptions)::Float64
@assert levmardata.lambda >= 0.
hessian, gradient = gethessgrad(data.linsystem)
if levmardata.lambda == 0
levmardata.lambda = tr(hessian) ./ (size(hessian, 1) * 1e6)
levmardata.lambda = initlambda(hessian)
end
lastlambda = 0.
mu = 2.
Expand All @@ -145,7 +157,7 @@ function iterate!(levmardata::LevMarData, data, problem::NLLSProblem, options::N
data.timecost += @elapsed_ns cost_ = cost(problem.varnext, problem.costs)
data.costcomputations += 1
# Check for exit
if !(cost_ > data.bestcost) || (maximum(abs, data.linsystem.x) < options.dstep)
if !(cost_ > data.bestcost) || maximum(abs, data.linsystem.x) < options.dstep
# Success (or convergence) - update lambda
uniformscaling!(hessian, -lastlambda)
stepquality = (cost_ - data.bestcost) / (0.5 * fast_bAb(hessian, data.linsystem.x) + dot(gradient, data.linsystem.x))
Expand Down
4 changes: 2 additions & 2 deletions src/linearsystem.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,9 @@ struct MultiVariateLSsparse
hessian = SparseMatrixCSC{Float64, Int}(sparseindices.m, sparseindices.n, sparseindices.colptr, sparseindices.rowval, Vector{Float64}(undef, length(sparseindices.nzval)))
sparseindices = sparseindices.nzval
else
x = Vector{Float64}()
x = Vector{Float64}(undef, 1)
sparseindices = Vector{Int}()
hessian = spzeros(0, 0)
hessian = spzeros(1, 1)
end
ldlfac = ldl_analyze(hessian)
return new(A, zeros(Float64, blen), x, blockindices, boffsets, hessian, sparseindices, ldlfac)
Expand Down
7 changes: 5 additions & 2 deletions src/marginalize.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ function marginalize!(to::MultiVariateLS, from::MultiVariateLSsparse, fromblock
marginalize!(to, from, range, blocksz)
end
end
return to
end

function initcrop!(to::MultiVariateLSsparse, from::MultiVariateLSsparse, fromblock=length(to.A.rowblocksizes)+1)
Expand All @@ -100,6 +101,7 @@ function initcrop!(to::MultiVariateLSsparse, from::MultiVariateLSsparse, fromblo
to.b .= view(from.b, 1:lastindex(to.b))
endind = from.A.indicestransposed.nzval[from.A.indicestransposed.colptr[fromblock]] - 1
view(to.A.data, 1:endind) .= view(from.A.data, 1:endind)
return to
end

function initcrop!(to::MultiVariateLSdense, from::MultiVariateLSsparse, fromblock=length(to.A.rowblockoffsets))
Expand All @@ -115,6 +117,7 @@ function initcrop!(to::MultiVariateLSdense, from::MultiVariateLSsparse, frombloc
block(to.A, row, col, lenr, lenc) .= reshape(view(from.A.data, (0:lenr*lenc-1) .+ from.A.indicestransposed.nzval[colind]), lenr, lenc)
end
end
return to
end

function constructcrop(from::MultiVariateLSsparse, fromblock, forcesparse=false)
Expand Down Expand Up @@ -150,10 +153,10 @@ function constructcrop(from::MultiVariateLSsparse, fromblock, forcesparse=false)
A = BlockSparseMatrix{Float64}(start-1, cropsparsity, blocksizes, blocksizes)

# Construct the sparse linear system
return MultiVariateLSsparse(A, from.blockindices[1:findfirst(isequal(fromblock), from.blockindices)])
return MultiVariateLSsparse(A, from.blockindices[1:findfirst(isequal(fromblock), from.blockindices)-1])
end
end

# Construct a dense linear system
return MultiVariateLSdense(toblocksizes, from.blockindices)
return MultiVariateLSdense(toblocksizes, from.blockindices[1:findfirst(isequal(fromblock), from.blockindices)-1])
end
21 changes: 11 additions & 10 deletions src/optimize.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,18 +78,19 @@ end

# The meat of an optimization
function optimizeinternal!(problem::NLLSProblem, options::NLLSOptions, data, iteratedata, callback)
# Timing initializations
# Do any preoptimization for the iterator
data.startcost = preoptimization(iteratedata, problem, options, data)::Float64
# Other initializations
fails = 0
data.iternum = 0
stoptime = data.starttime + options.maxtime
data.timeinit += Base.time_ns() - data.starttime
# Initialize the linear problem
data.timegradient += @elapsed_ns data.bestcost = costgradhess!(data.linsystem, problem.variables, problem.costs)
data.timegradient += @elapsed_ns cost = costgradhess!(data.linsystem, problem.variables, problem.costs)
data.gradientcomputations += 1
data.startcost = data.bestcost
data.bestcost = cost
data.startcost = max(cost, data.startcost)
# Do the iterations
fails = 0
cost = data.bestcost
converged = 0
data.iternum = 0
while true
data.iternum += 1
# Call the per iteration solver
Expand Down Expand Up @@ -158,11 +159,11 @@ function optimizeinternal!(problem::NLLSProblem, options::NLLSOptions, data, ite
end

# Optimizing variables one at a time (e.g. in alternation)
function optimizesinglesinternal!(problem::NLLSProblem, options::NLLSOptions, data::NLLSInternal{LST}, iteratedata, allcosts::CostStruct, costindices, indices, first) where {LST<:UniVariateLS}
function optimizesinglesinternal!(problem::NLLSProblem, options::NLLSOptions, data::NLLSInternal{LST}, iteratedata, allcosts::CostStruct, costindices, varindices, first) where {LST<:UniVariateLS}
iternum = data.iternum
while first <= length(indices)
while first <= length(varindices)
# Bail out if the variable size changes
ind = indices[first]
ind = varindices[first]
if nvars(problem.variables[ind]) != length(data.linsystem.b)
break
end
Expand Down
2 changes: 1 addition & 1 deletion src/utils.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@ function fast_bAb(A::SparseMatrixCSC, b::Vector)
return total
end

sparse_dense_decision(d, nnz) = (nnz * 32) < (25 * d * (d - 40)) # Threshold nnz = 25/32 * (d^2 - 40d)
sparse_dense_decision(d, nnz) = (nnz * 64) < (25 * d * (d - 40)) # Threshold nnz (for lower triangle) = 25/64 * (d^2 - 40d)

function block_sparse_nnz(sparsity, blocksizes)
# Compute the number of non-zeros in a block sparse matrix
Expand Down
2 changes: 1 addition & 1 deletion test/functional.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ NLLSsolver.ndeps(::RosenbrockA) = static(1)
NLLSsolver.nres(::RosenbrockA) = 1
NLLSsolver.varindices(::RosenbrockA) = SVector(1)
NLLSsolver.getvars(::RosenbrockA, vars::Vector) = (vars[1]::Float64,)
NLLSsolver.computeresidual(res::RosenbrockA, x) = res.a - x
NLLSsolver.computeresidual(res::RosenbrockA, x) = res.a * (1 - x)
Base.eltype(::RosenbrockA) = Float64
const rosenbrockrobustifier = NLLSsolver.Scaled(NLLSsolver.Huber2oKernel(1.6), 1.0)
NLLSsolver.robustkernel(::RosenbrockA) = rosenbrockrobustifier
Expand Down
4 changes: 2 additions & 2 deletions test/marginalize.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ using NLLSsolver, SparseArrays, StaticArrays, Test, Random
@test hess_d ≈ hess_s
@test to_d.b == to_s.b

# Check that the reduced systems gives the correct variable update
# Check that the reduced systems give the correct variable update
@test hessian \ gradient ≈ gtupdate
@test hess_d \ to_d.b ≈ gtupdate
@test hess_s \ to_s.b ≈ gtupdate
Expand All @@ -79,7 +79,7 @@ using NLLSsolver, SparseArrays, StaticArrays, Test, Random
@test hess_d ≈ hess_s
@test to_d.b == to_s.b

# Check that the reduced systems gives the correct variable update
# Check that the reduced systems give the correct variable update
@test hess_d \ to_d.b ≈ gtupdate
@test hess_s \ to_s.b ≈ gtupdate
end

2 comments on commit bfa26ce

@ojwoodford
Copy link
Owner Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/94136

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v3.3.2 -m "<description of version>" bfa26ce56ca887e2cc94aa714e77ee6601d0d7c0
git push origin v3.3.2

Please sign in to comment.