From 9df0f49aad284903b1bbcba400f2a0b1d7a6bb67 Mon Sep 17 00:00:00 2001 From: tmigot Date: Mon, 12 Aug 2024 09:29:51 -0400 Subject: [PATCH] review tuto --- docs/src/sparsepattern.md | 80 +++++++++++++++++++++++++++------------ 1 file changed, 56 insertions(+), 24 deletions(-) diff --git a/docs/src/sparsepattern.md b/docs/src/sparsepattern.md index 71c5572b..1d4ec3da 100644 --- a/docs/src/sparsepattern.md +++ b/docs/src/sparsepattern.md @@ -1,6 +1,6 @@ # Improve sparse derivatives -In this tutorial, we show a simple trick to dramatically improve the computation of sparse Jacobian and Hessian matrices. +In this tutorial, we show a simple trick to potentially improve the computation of sparse Jacobian and Hessian matrices. Our test problem is an academic investment control problem: @@ -11,23 +11,42 @@ Our test problem is an academic investment control problem: \end{aligned} ``` -Using a simple quadrature formula for the objective functional and a forward finite difference for the differential equation, one can obtain a finite-dimensional continuous optimisation problem. -One is implementation is available in the package [`OptimizationProblems.jl`](https://github.com/JuliaSmoothOptimizers/OptimizationProblems.jl). +Using a simple quadrature formula for the objective functional and a forward finite difference for the differential equation, one can obtain a finite-dimensional continuous optimization problem. +One implementation is available in the package [`OptimizationProblems.jl`](https://github.com/JuliaSmoothOptimizers/OptimizationProblems.jl). ```@example ex1 using ADNLPModels -using OptimizationProblems -using Symbolics using SparseArrays +T = Float64 n = 1000000 +N = div(n, 2) +h = 1 // N +x0 = 1 +gamma = 3 +function f(y; N = N, h = h) + @views x, u = y[1:N], y[(N + 1):end] + return 1 // 2 * h * sum((u[k] - 1) * x[k] + (u[k + 1] - 1) * x[k + 1] for k = 1:(N - 1)) +end +function c!(cx, y; N = N, h = h, gamma = gamma) + @views x, u = y[1:N], y[(N + 1):end] + for k = 1:(N - 1) + cx[k] = x[k + 1] - x[k] - 1 // 2 * h * gamma * (u[k] * x[k] + u[k + 1] * x[k + 1]) + end + return cx +end +lvar = vcat(-T(Inf) * ones(T, N), zeros(T, N)) +uvar = vcat(T(Inf) * ones(T, N), ones(T, N)) +xi = vcat(ones(T, N), zeros(T, N)) +lcon = ucon = vcat(one(T), zeros(T, N - 1)) + @elapsed begin - nlp = OptimizationProblems.ADNLPProblems.controlinvestment(n = n, hessian_backend = ADNLPModels.EmptyADbackend) + nlp = ADNLPModel!(f, xi, lvar, uvar, [1], [1], T[1], c!, lcon, ucon; hessian_backend = ADNLPModels.EmptyADbackend) end ``` -After adding the package `Symbolics.jl`, the `ADNLPModel` will automatically try to prepare AD-backend to compute sparse Jacobian and Hessian. +`ADNLPModel` will automatically try to prepare AD-backend to compute sparse Jacobian and Hessian. We disabled the Hessian computation here to focus the measurement on the Jacobian computation. The keyword argument `show_time = true` can also be passed to the problem's constructor to get more detailed information about the time used to prepare the AD backend. @@ -38,37 +57,50 @@ jac_nln(nlp, x) ``` However, it can be rather costly to determine for a given function the sparsity pattern of the Jacobian and the Lagrangian Hessian matrices. -The good news is that it can be quite easy to have a good approximation of this sparsity pattern dealing with problems like our optimal control investment problem, and problem with differential equations in the constraints in general. +The good news is that it can be quite easy to have a good approximation of this pattern dealing with problems like our optimal control investment problem, and problem with differential equations in the constraints in general. -The following example specialize the function `compute_jacobian_sparsity` to manually provide the sparsity pattern. +The following example instantiates the Jacobian backend while manually providing the sparsity pattern. ```@example ex2 using ADNLPModels -using OptimizationProblems -using Symbolics using SparseArrays +T = Float64 n = 1000000 N = div(n, 2) - -function ADNLPModels.compute_jacobian_sparsity(c!, cx, x0; n = n, N = N) - # S = Symbolics.jacobian_sparsity(c!, cx, x0) - S = spzeros(Bool, N - 1, n) - for i =1:(N - 1) - S[i, i] = true - S[i, i + 1] = true - S[i, N + i] = true - S[i, N + i + 1] = true +h = 1 // N +x0 = 1 +gamma = 3 +function f(y; N = N, h = h) + @views x, u = y[1:N], y[(N + 1):end] + return 1 // 2 * h * sum((u[k] - 1) * x[k] + (u[k + 1] - 1) * x[k + 1] for k = 1:(N - 1)) +end +function c!(cx, y; N = N, h = h, gamma = gamma) + @views x, u = y[1:N], y[(N + 1):end] + for k = 1:(N - 1) + cx[k] = x[k + 1] - x[k] - 1 // 2 * h * gamma * (u[k] * x[k] + u[k + 1] * x[k + 1]) end - return S + return cx end +lvar = vcat(-T(Inf) * ones(T, N), zeros(T, N)) +uvar = vcat(T(Inf) * ones(T, N), ones(T, N)) +xi = vcat(ones(T, N), zeros(T, N)) +lcon = ucon = vcat(one(T), zeros(T, N - 1)) @elapsed begin - nlp = OptimizationProblems.ADNLPProblems.controlinvestment(n = n, hessian_backend = ADNLPModels.EmptyADbackend) + J = spzeros(Bool, N - 1, n) + for i =1:(N - 1) + J[i, i] = true + J[i, i + 1] = true + J[i, N + i] = true + J[i, N + i + 1] = true + end + jac_back = ADNLPModels.SparseADJacobian(n, f, N - 1, c!, J) + nlp = ADNLPModel!(f, xi, lvar, uvar, [1], [1], T[1], c!, lcon, ucon; hessian_backend = ADNLPModels.EmptyADbackend, jacobian_backend = jac_back) end ``` -A similar Jacobian matrix is obtained at a lower price. +A similar Jacobian matrix is obtained potentially at a lower price. ```@example ex2 using NLPModels @@ -76,4 +108,4 @@ x = sqrt(2) * ones(n) jac_nln(nlp, x) ``` -The function `compute_hessian_sparsity(f, nvar, c!, ncon)` does the same for the Lagrangian Hessian. +The same can be done for the Lagrangian Hessian.