Certifying a local minimum

In a quantum physics problem, I need to find minimum of the univariate polynomial objective function obj(x) defined via the following code:

using DynamicPolynomials
using TSSOS
using SpecialFunctions
using LinearAlgebra

@polyvar x

H0 = Matrix(Diagonal([1, 1.1, 1.2]))

V = [
    0 0 1;
    0 0 1;
    1 1 0
]

O = Diagonal([0.1 1 0])

function exp_chebyshev(Δt::Real, Ω::AbstractMatrix, order::Integer)

    Tₙ₋₁ = I
    Tₙ  = Ω

    # The first two terms of Chebyshev series for exp
    series = besselj(0, Δt) * Tₙ₋₁ + 2 * besselj(1, Δt) * Tₙ

    for n=2:order
        Tₙ₊₁  = 2 * Ω * Tₙ + Tₙ₋₁

        series .+= 2 * besselj(n, Δt) * Tₙ₊₁

        (Tₙ, Tₙ₋₁) = (Tₙ₊₁, Tₙ) 
    end

    series
end

function real_poly(p::Polynomial)
    #=
    Real part of the polynomial
    =#
    sum(
        real(c) * m for (c, m) in zip(coefficients(p), monomials(p))# if ~isapproxzero(abs(c))
    )
end

𝓤 = exp_chebyshev(0.5, -im * (H0 + x * V), 10);

ψ = 𝓤 * [1; 0; 0]

obj = (1. - real_poly(ψ' * O * ψ)) ^ 2;

Theoretically, we know that obj(x) has a local minimum at x=0 and its global minimum of 0 is located at x ≠ 0. This is confirmed by simply plotting the polynomial

using Plots

ξ = range(-10, 10, 1000)
plot(ξ, obj.(ξ), label=nothing)
ylims!(0, 2)
xlabel!("x")
ylabel!("obj(x)")

obj(x)

If the optimization is run,

opt,sol,data = tssos_first(obj, variables(obj), solution = true);

previous_sol = sol
previous_opt = opt

while ~isnothing(sol)
    previous_sol = sol
    previous_opt = opt

    opt,sol,data = tssos_higher!(data; QUIET = true, solution = true)
end

previous_sol

a local minima is certified

************************TSSOS************************
TSSOS is launching...
Starting to compute the block structure...
------------------------------------------------------
The sizes of PSD blocks:
[11, 10]
[1, 1]
------------------------------------------------------
Obtained the block structure. The maximal size of blocks is 11.
Assembling the SDP...
There are 22 affine constraints.
Solving the SDP...
Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 22              
  Cones                  : 0               
  Scalar variables       : 1               
  Matrix variables       : 3               
  Integer variables      : 0               

Optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Freed constraints in eliminator : 0
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Problem
  Name                   :                 
  Objective sense        : max             
  Type                   : CONIC (conic optimization problem)
  Constraints            : 22              
  Cones                  : 0               
  Scalar variables       : 1               
  Matrix variables       : 3               
  Integer variables      : 0               

Optimizer  - threads                : 56              
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 22
Optimizer  - Cones                  : 2
Optimizer  - Scalar variables       : 5                 conic                  : 5               
Optimizer  - Semi-definite variables: 2                 scalarized             : 121             
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 234               after factor           : 234             
Factor     - dense dim.             : 0                 flops                  : 2.23e+04        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   1.7e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.00  
1   4.0e-01  2.4e-01  1.3e-01  1.17e+00   5.286737457e-01   4.802695354e-01   2.4e-01  0.00  
2   5.7e-02  3.4e-02  4.6e-03  1.04e+00   7.693185649e-01   7.821035527e-01   3.4e-02  0.01  
3   1.3e-02  7.9e-03  3.9e-04  1.41e+00   8.099679449e-01   8.131886037e-01   7.9e-03  0.01  
4   3.7e-03  2.2e-03  5.4e-05  1.54e+00   8.090341232e-01   8.096213331e-01   2.2e-03  0.01  
5   4.0e-04  2.4e-04  1.7e-06  1.15e+00   8.099467443e-01   8.100147297e-01   2.4e-04  0.01  
6   4.3e-05  2.6e-05  6.0e-08  1.08e+00   8.099953738e-01   8.100026900e-01   2.6e-05  0.01  
7   5.9e-06  3.5e-06  3.1e-09  9.90e-01   8.099992414e-01   8.100001773e-01   3.5e-06  0.01  
8   7.3e-07  4.4e-07  1.3e-10  1.02e+00   8.099999289e-01   8.100000483e-01   4.4e-07  0.01  
9   7.7e-08  4.6e-08  4.1e-12  1.06e+00   8.099999977e-01   8.100000110e-01   4.6e-08  0.01  
10  1.1e-08  6.7e-09  1.8e-13  1.09e+00   8.099999976e-01   8.099999997e-01   6.2e-09  0.01  
Optimizer terminated. Time: 0.01    

SDP solving time: 0.012239247 seconds.
optimum = 0.809999997628379
Global optimality certified!
No higher TSSOS hierarchy!

1-element Vector{Float64}:
 0.0

I played with this problem a lot with no success. Any help would be very much appreciated!

Hi Dibondar, I noticed that the coefficients of

obj = 3.360124005463267e-36x⁴⁰ - 2.063135343212682e-33x³⁸ + 7.737864395102542e-31x³⁶ - 2.0411920994033823e-28x³⁴ + 4.1216331637551674e-26x³² - 6.613486025637102e-24x³⁰ + 8.598482830688325e-22x²⁸ - 9.159301080688455e-20x²⁶ + 8.040326833190711e-18x²⁴ - 5.826362292706338e-16x²² + 3.478249843620919e-14x²⁰ - 1.7004343112432733e-12x¹⁸ + 6.729812063761517e-11x¹⁶ - 2.1131330221589752e-9x¹⁴ + 5.0499295994001044e-8x¹² - 8.198400590466459e-7x¹⁰ + 5.0724934755230575e-6x⁸ + 0.00013545758728621945x⁶ - 0.0040597146294860855x⁴ + 0.04496251247664391x² + 0.809999999999998

stretch across 36 orders of magnitude which is pathological for the SDP solver. So you should scale the problem to avoid this. For instance, let x -> 10x so that

obj = 3.360124005463267e4x^40 - 2.063135343212682e5x^38 + 7.737864395102542e5x^36 - 2.0411920994033823e6x^34 + 4.1216331637551674e6x^32 - 6.613486025637102e6x^30 + 8.598482830688325e6x^28 - 9.159301080688455e6x^26 + 8.040326833190711e6x^24 - 5.826362292706338e6x^22 + 3.478249843620919e6x^20 - 1.7004343112432733e6x^18 + 6.729812063761517e5x^16 - 2.1131330221589752e5x^14 + 5.0499295994001044e4x^12 - 8.198400590466459e3x^10 + 5.0724934755230575e2x^8 + 135.45758728621945x^6 - 40.597146294860855x^4 + 4.496251247664391x^2 + 0.809999999999998

Then it gives the correct global minimum 0 located at x = (10×) 0.8101203314765792.

wangjie212 / TSSOS

Certifying a local minimum #6