Skip to content

Commit

Permalink
IPM/HSD: use logical slicing instead of elementwise multiplication
Browse files Browse the repository at this point in the history
We now do basically this for a logical vector l: `dot(a[l], b[l])`,
instead of `dot(a .* l,  b .* l) as before.

This is as suggested here:
ds4dm#122 (comment)

Helps with the performance of the BigFloat arithmetic.

A Julia script and a Unix shell script were used to conduct an
experiment for assesing the impact of this commit and the previous
commit on performance. The scripts and the resulting CSV file follow.

The benchmark experiment is conducted and the CSV created by removing
sources of system load on a computer and running the shell script four
times: once with the `init` command and three times with the `run`
command, each time checking out a different commit in the Tulip git
repo.

Unix (Bourne) shell script:
```sh

set -u

command=$1

julia_opts='-O3 --min-optlevel=3 --heap-size-hint=5G --depwarn=error --warn-overwrite=yes'

script=tulip_benchmark.jl

case "$command" in

init_csv)
  printf '%s,%s,%s,%s\n' 'Tulip version' estimator 'measurement type' value
  ;;

run)
  tulip_version=$2
  $PATH_TO_JULIA_BIN $julia_opts "$script" "$tulip_version"
  ;;

*)
  printf '%s\n' error 2>&1
  exit 1
  ;;

esac
```

Julia script:
```julia
const benchmark_seconds = 500
const polynomial_degree = 20
setprecision(BigFloat, 12 * 2^7)
using BenchmarkTools
import
  FindMinimaxPolynomial,  # v0.2.3
  Tulip,
  MathOptInterface
const FMP = FindMinimaxPolynomial
const MMX = FMP.Minimax
const PPTI = FMP.PolynomialPassingThroughIntervals
const NE = FMP.NumericalErrorTypes
const to_poly = FMP.ToSparsePolynomial.to_sparse_polynomial
const mmx = MMX.minimax_polynomial
const error_type_relative = NE.RelativeError()
const MOI = MathOptInterface
const itv_max_err = FMP.ApproximateInfinityNorm.interval_max_err

function make_lp()
  lp = Tulip.Optimizer{BigFloat}()

  # Remove iteration limit just in case
  MOI.set(lp, MOI.RawOptimizerAttribute("IPM_IterationsLimit"), 2000)

  # Disable presolve, speeds things up
  #MOI.set(lp, MOI.RawOptimizerAttribute("Presolve_Level"), 0)

  lp
end

const itv = (-big"2.0"^-3, big"45.0")

odd_monomials(n::Int) = 1:2:n

sind_mmx(n::Int) =
  mmx(
    make_lp,
    sind,
    (itv,),
    odd_monomials(n),

    # Small factor to have less variance in the results
    initial_perturb_factor = 1//(2^20),

    # We're benchmarking LP, so disable other stuff
    worst_segments_density = 5,
    worst_segments_breadth_limit = 2,
    worst_segments_depth_ratio = 1/2,

    # Exit right after the first step
    exit_condition = true,
  )

function report(estimator; benchmark, benchmark_name)
  b = estimator(benchmark)
  println("$benchmark_name,$estimator,time,$(b.time)")
  println("$benchmark_name,$estimator,gctime,$(b.gctime)")
  println("$benchmark_name,$estimator,memory,$(b.memory)")
  println("$benchmark_name,$estimator,allocs,$(b.allocs)")
end

function report(;benchmark, benchmark_name)
  for quantile in (minimum, median, maximum)
    report(
      quantile,
      benchmark = benchmark,
      benchmark_name = benchmark_name,
    )
  end
end

report(
  benchmark = (@benchmark sind_mmx(polynomial_degree) seconds=benchmark_seconds),
  benchmark_name = first(ARGS),
)
```

CSV results:
```csv
Tulip version,estimator,measurement type,value
v0.9.5,minimum,time,2.63759609e8
v0.9.5,minimum,gctime,3.4505713e7
v0.9.5,minimum,memory,495299296
v0.9.5,minimum,allocs,3629874
v0.9.5,median,time,4.2594161e8
v0.9.5,median,gctime,1.3250321e8
v0.9.5,median,memory,495299296
v0.9.5,median,allocs,3629874
v0.9.5,maximum,time,4.55935021e8
v0.9.5,maximum,gctime,1.40426286e8
v0.9.5,maximum,memory,495299296
v0.9.5,maximum,allocs,3629874
MutableArithmetics for IPM/HSD,minimum,time,2.57993117e8
MutableArithmetics for IPM/HSD,minimum,gctime,2.9403466e7
MutableArithmetics for IPM/HSD,minimum,memory,442052896
MutableArithmetics for IPM/HSD,minimum,allocs,3238720
MutableArithmetics for IPM/HSD,median,time,4.22323273e8
MutableArithmetics for IPM/HSD,median,gctime,1.282365305e8
MutableArithmetics for IPM/HSD,median,memory,442052896
MutableArithmetics for IPM/HSD,median,allocs,3238720
MutableArithmetics for IPM/HSD,maximum,time,4.56330849e8
MutableArithmetics for IPM/HSD,maximum,gctime,1.57061172e8
MutableArithmetics for IPM/HSD,maximum,memory,442052896
MutableArithmetics for IPM/HSD,maximum,allocs,3238720
IPM/HSD: use logical slicing ...,minimum,time,2.40996648e8
IPM/HSD: use logical slicing ...,minimum,gctime,2.5335783e7
IPM/HSD: use logical slicing ...,minimum,memory,386588512
IPM/HSD: use logical slicing ...,minimum,allocs,2833356
IPM/HSD: use logical slicing ...,median,time,3.76039574e8
IPM/HSD: use logical slicing ...,median,gctime,1.06930941e8
IPM/HSD: use logical slicing ...,median,memory,386588512
IPM/HSD: use logical slicing ...,median,allocs,2833356
IPM/HSD: use logical slicing ...,maximum,time,4.00347376e8
IPM/HSD: use logical slicing ...,maximum,gctime,1.27260987e8
IPM/HSD: use logical slicing ...,maximum,memory,386588512
IPM/HSD: use logical slicing ...,maximum,allocs,2833356
```

Fixes ds4dm#122
  • Loading branch information
nsajko committed Apr 20, 2023
1 parent 29568d2 commit be697e6
Show file tree
Hide file tree
Showing 3 changed files with 13 additions and 12 deletions.
1 change: 1 addition & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ Krylov = "0.8, 0.9"
LDLFactorizations = "0.8, 0.9, 0.10"
LinearOperators = "2.0"
MathOptInterface = "1"
MutableArithmetics = "1.2"
QPSReader = "0.2"
TimerOutputs = "0.5.6"
julia = "1.6"
12 changes: 6 additions & 6 deletions src/IPM/HSD/HSD.jl
Original file line number Diff line number Diff line change
Expand Up @@ -112,8 +112,8 @@ function compute_residuals!(hsd::HSD{T}
(
(dat.c, pt.x),
(dat.b, pt.y),
(dat.l .* dat.lflag, pt.zl),
(dat.u .* dat.uflag, pt.zu),
(dat.l[dat.lflag], pt.zl[dat.lflag]),
(dat.u[dat.uflag], pt.zu[dat.uflag]),
),
(
1, -1, -1, 1,
Expand All @@ -133,8 +133,8 @@ function compute_residuals!(hsd::HSD{T}
dot_buf,
(
(dat.b, pt.y),
(dat.l .* dat.lflag, pt.zl),
(dat.u .* dat.uflag, pt.zu),
(dat.l[dat.lflag], pt.zl[dat.lflag]),
(dat.u[dat.uflag], pt.zu[dat.uflag]),
),
(
1, 1, -1,
Expand Down Expand Up @@ -209,8 +209,8 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher
dot_buf,
(
(dat.b, pt.y),
(dat.l .* dat.lflag, pt.zl),
(dat.u .* dat.uflag, pt.zu),
(dat.l[dat.lflag], pt.zl[dat.lflag]),
(dat.u[dat.uflag], pt.zu[dat.uflag]),
),
(
1, 1, -1,
Expand Down
12 changes: 6 additions & 6 deletions src/IPM/HSD/step.jl
Original file line number Diff line number Diff line change
Expand Up @@ -69,8 +69,8 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst
h0 = buffered_dot_weighted_sum!!(
dot_buf,
(
(dat.l .* dat.lflag, (dat.l .* θl) .* dat.lflag),
(dat.u .* dat.uflag, (dat.u .* θu) .* dat.uflag),
(dat.l[dat.lflag], (dat.l .* θl)[dat.lflag]),
(dat.u[dat.uflag], (dat.u .* θu)[dat.uflag]),
((@. (c + (θl * dat.l) * dat.lflag + (θu * dat.u) * dat.uflag)), hx),
(b, hy),
),
Expand Down Expand Up @@ -225,10 +225,10 @@ function solve_newton_system!(Δ::Point{T, Tv},
buffered_dot_weighted_sum!!(
dot_buf,
(
((ξxzl ./ pt.xl) .* dat.lflag, dat.l .* dat.lflag), # l'(Xl)^-1 * ξxzl
((ξxzu ./ pt.xu) .* dat.uflag, dat.u .* dat.uflag),
(((pt.zl ./ pt.xl) .* ξl) .* dat.lflag, dat.l .* dat.lflag),
(((pt.zu ./ pt.xu) .* ξu) .* dat.uflag, dat.u .* dat.uflag),
((ξxzl ./ pt.xl)[dat.lflag], dat.l[dat.lflag]), # l'(Xl)^-1 * ξxzl
((ξxzu ./ pt.xu)[dat.uflag], dat.u[dat.uflag]),
(((pt.zl ./ pt.xl) .* ξl)[dat.lflag], dat.l[dat.lflag]),
(((pt.zu ./ pt.xu) .* ξu)[dat.uflag], dat.u[dat.uflag]),
),
(
-1, 1, -1, -1,
Expand Down

0 comments on commit be697e6

Please sign in to comment.