JuliaSmoothOptimizers / QuadraticModels.jl

Data structures for linear and quadratic optimization problems based on NLPModels.jl
Other
16 stars 10 forks source link

Track allocs #106

Open tmigot opened 2 years ago

tmigot commented 2 years ago

There seems to be two functions of the API allocating

I used the following script to track allocations in QuadraticModels

using Pkg
Pkg.activate(".")

# stdlib
using LinearAlgebra, Printf, SparseArrays, Test

# our packages
using ADNLPModels,
  LinearOperators,
  NLPModels,
  NLPModelsModifiers,
  NLPModelsTest, # main version of NLPModelsTest
  QPSReader,
  QuadraticModels,
  SparseMatricesCOO

function only_nonzeros(table)
  for k in keys(table)
    if table[k] == 0
      pop!(table, k)
    end
  end
  return table
end

# Definition of quadratic problems
qp_problems_Matrix = ["bndqp", "eqconqp"]
qp_problems_COO = ["uncqp", "ineqconqp"]
for qp in [qp_problems_Matrix; qp_problems_COO]
  include(joinpath("problems", "$qp.jl"))
end

for problem in qp_problems_Matrix
  @info "Checking allocs of dense problem $(problem)_QPSData"
  nlp_qps = eval(Symbol(problem * "_QPSData"))()
  print_nlp_allocations(nlp_qps, only_nonzeros(test_allocs_nlpmodels(nlp_qps)))
end

for problem in qp_problems_Matrix
  @info "Checking allocs of dense problem $(problem)_QP_dense"
  nlp_qm_dense = eval(Symbol(problem * "_QP_dense"))()
  print_nlp_allocations(nlp_qm_dense, only_nonzeros(test_allocs_nlpmodels(nlp_qm_dense)))
end

for problem in qp_problems_Matrix
  @info "Checking allocs of dense problem $(problem)_QP_sparse"
  nlp_qm_sparse = eval(Symbol(problem * "_QP_sparse"))()
  print_nlp_allocations(nlp_qm_sparse, only_nonzeros(test_allocs_nlpmodels(nlp_qm_sparse)))
end

for problem in qp_problems_Matrix
  @info "Checking allocs of dense problem $(problem)_QP_symmetric"
  nlp_qm_symmetric = eval(Symbol(problem * "_QP_symmetric"))()
  print_nlp_allocations(nlp_qm_symmetric, only_nonzeros(test_allocs_nlpmodels(nlp_qm_symmetric)))
end

for problem in qp_problems_COO
  @info "Checking allocs of COO problem $(problem)_QPSData"
  nlp_qps = eval(Symbol(problem * "_QPSData"))()
  print_nlp_allocations(nlp_qps, only_nonzeros(test_allocs_nlpmodels(nlp_qps)))
end

for problem in qp_problems_COO
  @info "Checking allocs of COO problem $(problem)_QP"
  nlp_qm_dense = eval(Symbol(problem * "_QP"))()
  print_nlp_allocations(nlp_qm_dense, only_nonzeros(test_allocs_nlpmodels(nlp_qm_dense)))
end

for problem in NLPModelsTest.nlp_problems
  @info "Testing allocs of quadratic approximation of problem $problem"
  nlp = eval(Symbol(problem))()
  x = nlp.meta.x0
  nlp_qm = QuadraticModel(nlp, x)
  print_nlp_allocations(nlp_qm, only_nonzeros(test_allocs_nlpmodels(nlp_qm)))
end

and the results

[ Info: Checking allocs of dense problem bndqp_QPSData
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0

[ Info: Checking allocs of dense problem eqconqp_QPSData
  Problem name: Generic
                        obj: ████████████████████ 496.0
                hess_coord!: ████████████████████ 496.0
            hess_lag_coord!: ████████████████████ 496.0

[ Info: Checking allocs of dense problem bndqp_QP_dense
  Problem name: bndqp_QP
                        obj: ████████████████████ 80.0

[ Info: Checking allocs of dense problem eqconqp_QP_dense
  Problem name: eqconqp_QP
                        obj: ████████████████████ 496.0

[ Info: Checking allocs of dense problem bndqp_QP_sparse
  Problem name: bndqp_QP
                        obj: ████████████████████ 80.0

[ Info: Checking allocs of dense problem eqconqp_QP_sparse
  Problem name: eqconqp_QP
                        obj: ████████████████████ 496.0

[ Info: Checking allocs of dense problem bndqp_QP_symmetric
  Problem name: bndqp_QP
                        obj: ████████████████████ 80.0

[ Info: Checking allocs of dense problem eqconqp_QP_symmetric
  Problem name: eqconqp_QP
                        obj: ████████████████████ 496.0

[ Info: Checking allocs of COO problem uncqp_QPSData
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0

[ Info: Checking allocs of COO problem ineqconqp_QPSData
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0
            hess_lag_coord!: ████████████████████ 80.0

[ Info: Checking allocs of COO problem uncqp_QP
  Problem name: uncqp_QP
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0

[ Info: Checking allocs of COO problem ineqconqp_QP
  Problem name: ineqconqp_QP
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0
            hess_lag_coord!: ████████████████████ 80.0

[ Info: Testing allocs of quadratic approximation of problem BROWNDEN
  Problem name: Generic
                        obj: ██████████████⋅⋅⋅⋅⋅⋅ 96.0
                hess_coord!: ████████████████████ 144.0

[ Info: Testing allocs of quadratic approximation of problem HS5
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0

[ Info: Testing allocs of quadratic approximation of problem HS6
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████⋅⋅⋅⋅ 64.0
            hess_lag_coord!: ████████████████⋅⋅⋅⋅ 64.0

[ Info: Testing allocs of quadratic approximation of problem HS10
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0
            hess_lag_coord!: ████████████████████ 80.0

[ Info: Testing allocs of quadratic approximation of problem HS11
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0
            hess_lag_coord!: ████████████████████ 80.0

[ Info: Testing allocs of quadratic approximation of problem HS13
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0
            hess_lag_coord!: ████████████████████ 80.0

[ Info: Testing allocs of quadratic approximation of problem HS14
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████████ 80.0
            hess_lag_coord!: ████████████████████ 80.0

[ Info: Testing allocs of quadratic approximation of problem LINCON
  Problem name: Generic
                        obj: ████████████████████ 176.0
                hess_coord!: ████████████████████ 176.0
            hess_lag_coord!: ████████████████████ 176.0

[ Info: Testing allocs of quadratic approximation of problem LINSV
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████⋅⋅⋅⋅ 64.0
            hess_lag_coord!: ████████████████⋅⋅⋅⋅ 64.0

[ Info: Testing allocs of quadratic approximation of problem MGH01Feas
  Problem name: Generic
                        obj: ████████████████████ 80.0
                hess_coord!: ████████████████⋅⋅⋅⋅ 64.0
            hess_lag_coord!: ████████████████⋅⋅⋅⋅ 64.0
tmigot commented 2 years ago

After #107 only obj allocate.

geoffroyleconte commented 2 years ago

Fixing obj would require writing a function that computes x' * H * x without allocating (for all H types).

tmigot commented 2 years ago

Isn't there a 3-arguments dot?

help?> dot
search: dot @__dot__ stdout DEPOT_PATH adjoint Adjoint adjoint! fieldcount codepoint fieldoffset bound_constrained isdisjoint lowrankdowndate lowrankdowndate! denominator VoidListofStates add_to_list! PKGMODE_PROJECT modifyproperty!

  dot(x, y)
  x ⋅ y

  Compute the dot product between two vectors. For complex vectors, the first vector is conjugated.

  dot also works on arbitrary iterable objects, including arrays of any dimension, as long as dot is defined on the elements.

  dot is semantically equivalent to sum(dot(vx,vy) for (vx,vy) in zip(x, y)), with the added restriction that the arguments must have equal lengths.

  x ⋅ y (where ⋅ can be typed by tab-completing \cdot in the REPL) is a synonym for dot(x, y).

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> dot([1; 1], [2; 3])
  5

  julia> dot([im; im], [1; 1])
  0 - 2im

  julia> dot(1:5, 2:6)
  70

  julia> x = fill(2., (5,5));

  julia> y = fill(3., (5,5));

  julia> dot(x, y)
  150.0

  ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────  

  dot(x, A, y)

  Compute the generalized dot product dot(x, A*y) between two vectors x and y, without storing the intermediate result of A*y. As for the two-argument dot(_,_), this acts recursively. Moreover, for complex vectors, the first vector is        
  conjugated.

  │ Julia 1.4
  │
  │  Three-argument dot requires at least Julia 1.4.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> dot([1; 1], [1 2; 3 4], [2; 3])
  26

  julia> dot(1:5, reshape(1:25, 5, 5), 2:6)
  4850

  julia> ⋅(1:5, reshape(1:25, 5, 5), 2:6) == dot(1:5, reshape(1:25, 5, 5), 2:6)
  true
geoffroyleconte commented 2 years ago

Oh I didn't know about this function. It would have to be added to SparseMatricesCOO. Also I would keep the existing behaviour in case H is a LinearOperator.

tmigot commented 2 years ago

Ok, so if I understand. We should implement this 3-args dot in SparseMatricesCOO and in LinearOperators?

geoffroyleconte commented 2 years ago

Yes for SparseMatricesCOO. For LinearOperators this is complicated, I would keep the existing code even though it allocates.

dpo commented 2 years ago

How about preallocating storage in the quadratic model to store H * x? It could be reused for the gradient.