Bose–Einstein Condensates

Bose–Einstein Condensates#

Problem Description#

In Bose–Einstein condensates (BEC), under some proper discretization, such as finite difference, sine pseudospectral and Fourier pseudospectral methods, we obtain its discrete version as

\[\begin{split} \begin{aligned} \min_{X \in \mathbb{C}^{n}}\quad & \frac{1}{2} x^H A x + \frac{\alpha}{2} \sum_{i = 1}^n |x_i|^4 \\ \text{s. t.} \quad & ||x||_2 = 1, \end{aligned} \end{split}\]

where \(A\) is a Hermitian matirx, and \(\alpha\) is a parameter. Consider a simplified cases where \(x\) and \(A\) are both real, we each the following optimization problem over the sphere,

\[\begin{split} \begin{aligned} \min_{X \in \mathbb{R}^{n}}\quad & \frac{1}{2} x^H A x + \frac{\alpha}{2} \sum_{i = 1}^n |x_i|^4 \\ \text{s. t.} \quad & ||x||_2 = 1, \end{aligned} \end{split}\]

In this part, we aim to show how to solve these problems with cdopt package.

Importing modules#

We first import all the necessary modules for this optimization problem.

import cdopt 
import numpy as np
import scipy as sp
from scipy.sparse import csr_matrix, diags
from scipy.sparse.linalg import spsolve
import time
import autograd 
import autograd.numpy as anp

Generating datas#

We generate necessary data. Notice that \(L\) is a sparse matrix, we apply the functions from scipy.sparse to accelerate the computation.

n = 1000
alpha = 1
A = np.ones((n,n))

Set functions and problems#

Then we set the objective function and the Stiefel manifold. Notice that all the existing AD packages has limited support for operations on sparse matrices, we manually define the gradient and Hessian-vector product of the objective function.

def obj_fun(X):
    return 0.5 * anp.sum( X* (A @ X) ) + alpha * anp.sum( X ** 4 )


M = cdopt.manifold_np.sphere_np((n,1))   # The sphere

Describe the optimization problem#

The optimization problem can be described by the manifold and the drivatives of objective function.

problem_test = cdopt.core.problem(M, obj_fun, beta = 'auto')  # describe the optimization problem

Apply optimization solvers#

After describe the optimization problem, we can directly function value, gradient and Hessian-vector product from the cdopt.core.Problem class.

# the vectorized function value, gradient and Hessian-vector product of the constraint dissolving function. Their inputs are numpy 1D array, and their outputs are float or numpy 1D array.
cdf_fun_np = problem_test.cdf_fun_vec_np   
cdf_grad_np = problem_test.cdf_grad_vec_np 
cdf_hvp_np = problem_test.cdf_hvp_vec_np


## Apply limit memory BFGS solver from scipy.minimize 
from scipy.optimize import fmin_bfgs, fmin_cg, fmin_l_bfgs_b, fmin_ncg
Xinit = problem_test.Xinit_vec_np # set initial point

# optimize by L-BFGS method
t_start = time.time()
out_msg = sp.optimize.minimize(cdf_fun_np, Xinit,method='L-BFGS-B',jac = cdf_grad_np, options={'disp': None, 'maxcor': 10, 'ftol': 0, 'gtol': 1e-06, 'eps': 0e-08,})
t_end = time.time() - t_start

# Statistics
feas = M.Feas_eval(M.v2m(M.array2tensor(out_msg.x)))   # Feasibility
stationarity = np.linalg.norm(out_msg['jac'],2)   # stationarity

result_lbfgs = [out_msg['fun'], out_msg['nit'], out_msg['nfev'],stationarity,feas, t_end]

# print results
print('Solver   fval         iter   f_eval   stationarity   feaibility     CPU time')
print('& L-BFGS & {:.2e}  & {:}  & {:}    & {:.2e}     & {:.2e}     & {:.2f} \\\\'.format(*result_lbfgs))

Solver   fval         iter   f_eval   stationarity   feaibility     CPU time
& L-BFGS & 1.00e-03  & 384  & 475    & 6.48e-06     & 1.87e-06     & 0.71 \\

# optimize by Newton GLTR trust-region method
t_start = time.time()
out_msg = sp.optimize.minimize(cdf_fun_np, Xinit, method='trust-krylov',jac = cdf_grad_np, hessp = cdf_hvp_np, options={'gtol': 1e-06, 'disp': False})
t_end = time.time() - t_start

# Statistics
feas = M.Feas_eval(M.v2m(M.array2tensor(out_msg.x)))   # Feasibility
stationarity = np.linalg.norm(out_msg['jac'],2)   # stationarity

result_cg = [out_msg['fun'], out_msg['nit'], out_msg['nfev'],stationarity,feas, t_end]

# print results
print('Solver   fval         iter   f_eval   stationarity   feaibility     CPU time')
print('& TR-KRY & {:.2e}  & {:}  & {:}    & {:.2e}     & {:.2e}     & {:.2f} \\\\'.format(*result_cg))

Solver   fval         iter   f_eval   stationarity   feaibility     CPU time
& TR-KRY & 1.00e-03  & 17  & 18    & 5.03e-07     & 4.97e-07     & 0.39 \\

# optimize by Newton conjugate gradient trust-region method  
t_start = time.time()
out_msg = sp.optimize.minimize(cdf_fun_np, Xinit, method='trust-ncg',jac = cdf_grad_np, hessp = cdf_hvp_np, options={'gtol': 1e-06, 'disp': False})
t_end = time.time() - t_start

# Statistics
feas = M.Feas_eval(M.v2m(M.array2tensor(out_msg.x)))   # Feasibility
stationarity = np.linalg.norm(out_msg['jac'],2)   # stationarity

result_cg = [out_msg['fun'], out_msg['nit'], out_msg['nfev'],stationarity,feas, t_end]

# print results
print('Solver   fval         iter   f_eval   stationarity   feaibility     CPU time')
print('& TR-NCG & {:.2e}  & {:}  & {:}    & {:.2e}     & {:.2e}     & {:.2f} \\\\'.format(*result_cg))

Solver   fval         iter   f_eval   stationarity   feaibility     CPU time
& TR-NCG & 1.00e-03  & 28  & 29    & 1.16e-08     & 1.23e-08     & 0.39 \\

Reference#

  1. Hu J, Jiang B, Liu X, et al. A note on semidefinite programming relaxations for polynomial optimization over a single sphere[J]. Science China Mathematics, 2016, 59(8): 1543-1560.

Contents