# Iterative algorithm with structured diagonal Hessian approximation for
solving nonlinear least squares problems

Aliyu Muhammed Awwal1,3 ,  Poom Kumam1,2∗ and  Hassan Mohammad4 1KMUTTFixed
Point Research Laboratory, Department of Mathematics, Room SCL 802 Fixed Point
Laboratory, Science Laboratory Building, Faculty of Science, King Mongkut’s
University of Technology Thonburi (KMUTT), 126 Pracha-Uthit Road, Bang Mod,
Thrung Khru, Bangkok 10140, Thailand  2KMUTT-Fixed Point Theory and
Applications Research Group, Theoretical and Computational Science Center
(TaCS), Science Laboratory Building, Faculty of Science, King Mongkut’s
University of Technology Thonburi (KMUTT), 126 Pracha-Uthit Road, Bang Mod,
Thrung Khru, Bangkok 10140, Thailand  3Department of Mathematics, Faculty of
Science, Gombe State University, Gombe, Nigeria  4Department of Mathematical
Sciences, Faculty of Physical Sciences, Bayero University, Kano. Kano, Nigeria
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

Nonlinear least squares problems are special class of unconstrained
optimization problems in which their gradient and Hessian have special
structures. In this paper, we exploit these structures and proposed a matrix
free algorithm with diagonal Hessian approximation for solving nonlinear least
squares problems. We devise appropriate safeguarding strategies to ensure the
Hessian matrix is positive definite throughout the iteration process. The
proposed algorithm generates descent direction and is globally convergent.
Preliminary numerical experiments shows that the proposed method is
competitive with a recent developed similar methods.

###### Key words and phrases:

Nonlinear least-squares problems; Large-scale problems; Jacobian-free
strategy; Global convergence

###### 2010 Mathematics Subject Classification:

90C30; 65K05; 49M37

##  1\. Introduction

In this paper, we consider nonlinear least squares problems, the special class
of unconstrained optimization problems, of the form

(1.1) $\min_{x\in\mathbb{R}^{n}}f(x),~{}~{}~{}f(x)=\frac{1}{2}\|F(x)\|^{2},$

where $F(x)=(F_{1}(x),\cdots,F_{m}(x))^{T}$ and each residual
$F_{i}:\mathbb{R}^{n}\to\mathbb{R},$ $i=1,\cdots,m$ (usually $m\geq n$), is
twice continuously differentiable function. Let $J(x)\in\mathbb{R}^{m\times
n}$ denotes Jacobian of the residual function $F(x)$ and $g(x),$ denote the
gradient of the objective function $f$, $\nabla f(x_{k})$, and $H(x)$ denote
the Hessian of the objective function $\nabla^{2}f(x)$. The gradient and
Hessian of problem (1.1) have special structures and are respectively given by

(1.2) $g(x):=\sum_{i=1}^{m}F_{i}(x)\nabla F_{i}(x)=J(x)^{T}F(x)$ (1.3)
$H(x):=\sum_{i=1}^{m}\nabla F_{i}(x)\nabla
F_{i}(x)^{T}+\sum_{i=1}^{m}F_{i}(x)\nabla^{2}F_{i}(x)=J(x)^{T}J(x)+C(x),$

where $F_{i}$ is the $i$th component of $F$, $\nabla^{2}F_{i}(x)$ is its
Hessian, and $C(x)$ is a square matrix representing the second term of the
Hessian.

Nonlinear least squares problems have been studied extensively, and many
iterative algorithms for solving them have been proposed. These generally fall
into two categories, namely, general unconstrained optimization algorithms
that includes Newton’s method and quasi-Newton methods; and special methods,
that take the special structure of the problem into account, which constitute
Gauss-Newton method, Levenberg-Marquardt method and Structured quasi-Newton
methods (see [22, 25, 26, 23, 6]). For a brief survey of methods for
addressing nonlinear least squares problems, interested reader may refer to
the recent articles by Mohammad et al. [17] and Yuan [24].

The study of efficient algorithm for nonlinear least squares (NLS) problems is
important because of its numerous areas of applications such as data fitting,
optimal control, parameter estimation, experimental design, data assimilation,
and imaging problems (see [7, 10, 13, 2, 1, 21]). For instance, it is often
common to measure the discrepancy between a proposed parametrized model and
the observed behavior of a given system. To select values for the parameters
that best match the model to the data, it is usual to minimize the sum of the
squares of the residuals $F_{i}^{{}^{\prime}s}$. The special structure of
problem (1.1) can always be explored to devise efficient algorithms for
obtaining its solution. For example Kobayashi et al. [11] exploits the
structure of the nonlinear least squares by introducing a class of matrix-free
structured methods that falls into the category of conjugate gradient
algorithms with modified secant condition for solving nonlinear least squares
problems. Motivated by their idea, Dehghani and Mahdavi-Amiri [3] proposed a
modified secant relation specifically to get more information of the Hessian
of the nonlinear least squares objective function. Furthermore, they proposed
another class of conjugate gradient methods for addressing nonlinear least
squares problems. In another attempt, but different approach, Mohammad and
Waziri [16] proposed two structured Barzilai-Borwein step sizes for solving
nonlinear least squares.

Recently, Mohammad and Sandra [15] proposed a diagonal Hessian approximation
method for nonlinear least squares problems in which the diagonal
approximation of the Hessian of the objective function is obtained using a
structured secant condition that have some information of the exact Hessian.
However, as a final remarks, the authors comment on the need for further
research that investigate a better approximation of the Hessian matrix that
involved its special structure.

We feel that approximating the first and second terms of the Hessian matrix
(1.3) will lead to substantial lost of information about the Hessian. In this
paper, we proposed a diagonal Hessian with better approximation by exploiting
the special structure of problem (1.1) and obtained a matrix-free algorithm.
The main difference between our method and the method in [15], is that, in
building our diagonal matrix, we take the whole information of the first term
of (1.3) and approximate its second term. By this, our diagonal matrix
contains more information than the one proposed in [15]. Our proposed method
generates descent directions and is globally convergent.

The remaining of this paper is organized as follows. In Section 2, we present
the proposed method and its algorithm. The convergence analysis is discussed
in Section 3 and in Section 4, we give numerical experiments. Throughout this
article, we use the following notations for the objective
function$f(x_{k})=f_{k},\text{ for the residual }~{}F(x_{k})=F_{k},\text{ for
any matrix}~{}A(x_{k})=A_{k},~{}\text{ and }\|\cdot\|$ for the Euclidean norm
of vectors and the induced 2-norm of matrices.

##  2\. Proposed method

An important concept of a structured quasi-Newton method for nonlinear least
squares is the structure principle [4]. Now, we provide the necessary elements
to develop our proposed diagonal Hessian approximation, taking the special
structures of the gradient and Hessian of the problem (1.1) into account.
Consider the second term of the Hessian matrix (1.3). Suppose that at certain
iteration $k-1,~{}~{}k\geq 1$ the second term of the Hessian matrix (1.3) is

(2.1) $C(x_{k-1})=\sum_{i=1}^{m}F_{i}(x_{k-1})\nabla^{2}F_{i}(x_{k-1}),$

so that the updating matrix $C(x_{k})$ which satisfies the secant equation

$C(x_{k})s_{k-1}=y_{k-1},$

can be obtained as follows.
The Taylor’s expansion of $\nabla F_{i}(x_{k-1})$ can be written as

(2.2) $\nabla F_{i}(x_{k-1})=\nabla
F_{i}(x_{k})+\nabla^{2}F_{i}(x_{k})^{T}(x_{k-1}-x_{k})+\textit{{o}}(\|x_{k-1}-x_{k}\|),$

where $\textit{{o}}:\mathbb{R}_{+}\to\mathbb{R}^{n}$ such that for each
$i=1,\cdots,n,$ $\lim\limits_{\xi\to 0}\frac{o^{i}(\xi)}{\xi}=0.$
Let $s_{k-1}=x_{k}-x_{k-1}.$ Multiplying (2.2) by $F_{i}(x_{k})$ and
rearranging, we have

(2.3) $F_{i}(x_{k})\nabla^{2}F_{i}(x_{k})^{T}s_{k-1}=F_{i}(x_{k})\nabla
F_{i}(x_{k})-F_{i}(x_{k})\nabla
F_{i}(x_{k-1})+F_{i}(x_{k})\textit{{o}}(\|s_{k-1}\|),$

By summing both sides of (2.3) for $i=1,\cdots,m$ we obtain

(2.4)
$C(x_{k})s_{k-1}=J(x_{k})^{T}F(x_{k})-J(x_{k-1})^{T}F(x_{k})+(F_{k}^{T}\textbf{1}_{m})\textit{{o}}(\|s_{k-1}\|).$

Putting (2.4) into (1.3) implies

(2.5)
$H_{k}s_{k-1}=J_{k}^{T}J_{k}s_{k-1}+(J_{k}-J_{k-1})^{T}F(x_{k})+(F_{k}^{T}\textbf{1}_{m})\textit{{o}}(\|s_{k-1}\|).$

Let $D_{k}\approx H_{k}$ such that $D_{k}$ is a diagonal matrix approximately
satisfying the secant equation

(2.6) $D_{k}s_{k-1}\approx y_{k-1},$

where $y_{k-1}=J_{k}^{T}J_{k}s_{k-1}+(J_{k}-J_{k-1})^{T}F(x_{k}).$ For
convenience, we denote the first and second terms of $y_{k-1}$ as

(2.7)
$\hat{y}_{k-1}=J_{k}^{T}J_{k}s_{k-1}~{}~{}\text{and}~{}~{}\overline{y}_{k-1}=(J_{k}-J_{k-1})^{T}F(x_{k}).$

The following Lemma comes from [15] and will be useful in defining the entries
of the diagonal matrix $D_{k}$ in view of the secant equation (2.6).

######  Lemma 2.1.

Let Let $D=\text{diag}(d)$ be a diagonal matrix in $\mathbb{R}^{n\times n},$
and let $c$ and $s$ be vectors in $\mathbb{R}^{n}.$ Then, the solution of the
constrained linear least-squares problem with simple bounds:

(2.8) $\min_{d\in\mathbb{R}^{n}}\frac{1}{2}\|diag(d)s-c\|^{2}$
$subject~{}~{}to~{}-d\leq 0$

is given by

(2.9)
${d^{i}}=\begin{cases}\frac{c^{i}}{s^{i}},&\text{if}~{}~{}\frac{c^{i}}{s^{i}}>0,\\\
0,&\text{if}~{}~{}\frac{c^{i}}{s^{i}}\leq
0,~{}\text{or}~{}s^{i}=0,\end{cases}i=1,\cdots,n.$

Let the diagonal matrix $D_{k}$ be decompose into two diagonal matrices i.e.
$D_{k}=diag(a_{k})+diag(b_{k})$, where $a_{k},~{}b_{k}$ are vectors
representing the first and second terms of $y_{k-1}$ respectively. For the
diagonal matrix to be positive definite, all the diagonal entries $a_{k}^{i}$
and $b_{k}^{i}$ for $i=1,\cdots,n,$ must be strictly positive. Now, applying
Lemma 2.1 to the secant equation (2.6) for $i=1,\cdots,n,$ Equation (2.9)
becomes

(2.10)
${a_{k}^{i}}=\begin{cases}\frac{\hat{y}_{k-1}^{i}}{s_{k-1}^{i}},&\text{if}~{}~{}\frac{\hat{y}_{k-1}^{i}}{s_{k-1}^{i}}>0,\\\
0,&\text{if}~{}~{}\frac{\hat{y}_{k-1}^{i}}{s_{k-1}^{i}}\leq
0,~{}\text{or}~{}s^{i}=0,\end{cases}i=1,\cdots,n.$

and

(2.11)
${b_{k}^{i}}=\begin{cases}\frac{\overline{y}_{k-1}^{i}}{s_{k-1}^{i}},&\text{if}~{}~{}\frac{\overline{y}_{k-1}^{i}}{s_{k-1}^{i}}>0,\\\
0,&\text{if}~{}~{}\frac{\overline{y}_{k-1}^{i}}{s_{k-1}^{i}}\leq
0,~{}\text{or}~{}s^{i}=0,\end{cases}i=1,\cdots,n.$

From the above Lemma 2.1, the diagonal matrix has nonnegative diagonal entries
which means it is positive semidefinite matrix. In order to ensure the
diagonal matrix is positive definite, in the next subsection we provide a
safeguard strategy similar to the one given in [15], that guarantee each
diagonal entries $a_{k}^{i}$ and $b_{k}^{i}$ for $i=1,\cdots n,$ is strictly
positive.

###  2.1. Safeguarding strategy

We consider the situation in which $\hat{y}_{k-1}^{i},s_{k-1}^{i}$ and
$\overline{y}_{k-1}^{i},s_{k-1}^{i}$ have different signs with
$s_{k-1}^{i}\neq 0.$ Let $\gamma\in(0,1)$ be a shrinking parameter and
$\rho>0$ be a tolerance for ensuring strictly values.

  * Case (a)

    * Suppose $s_{k-1}^{i}>0$.

    * (ai)

If $\hat{y}_{k-1}^{i}\leq 0,$ i.e., $(J_{k}^{T}J_{k}s_{k-1})^{i}\leq 0.$ Then
redefine $\hat{y}_{k-1}^{i}$ as

(2.12) $\hat{y}_{k-1}^{i}=\gamma\max\\{|(J_{k}^{T}J_{k}s_{k-1})^{i}|,\rho\\}$

so that $a_{k}^{i}=\frac{\hat{y}_{k-1}^{i}}{s_{k-1}^{i}}>0.$

    * (aii)

If $\overline{y}_{k-1}^{i}\leq 0,$ i.e.,
$(J_{k}^{T}F_{k})^{i}\leq(J_{k-1}^{T}F_{k})^{i}.$ Then redefine
$\overline{y}_{k-1}^{i}$ as

(2.13)
$\overline{y}_{k-1}^{i}=\gamma\left\\{\max\left\\{\max|(J_{k}^{T}F_{k})^{i}|,|(J_{k-1}^{T}F_{k})^{i}|\right\\},\rho\right\\}$

so that $b_{k}^{i}=\frac{\overline{y}_{k-1}^{i}}{s_{k-1}^{i}}>0.$

  * Case(b)

    * Suppose $s_{k-1}^{i}<0$.

    * (bi)

If $\hat{y}_{k-1}^{i}\geq 0,$ i.e., $(J_{k}^{T}J_{k}s_{k-1})^{i}\geq 0.$ Then
redefine $\hat{y}_{k-1}^{i}$ as

(2.14) $\hat{y}_{k-1}^{i}=-\gamma\max\\{(J_{k}^{T}J_{k}s_{k-1})^{i},\rho\\}$

so that $a_{k}^{i}=\frac{\hat{y}_{k-1}^{i}}{s_{k-1}^{i}}>0.$

    * (bii)

If $\overline{y}_{k-1}^{i}\geq 0,$ i.e.,
$(J_{k}^{T}F_{k})^{i}\geq(J_{k-1}^{T}F_{k})^{i}.$ Then redefine
$\overline{y}_{k-1}^{i}$ as

(2.15)
$\overline{y}_{k-1}^{i}=-\gamma\left\\{\max\left\\{\max|(J_{k}^{T}F_{k})^{i}|,|(J_{k-1}^{T}F_{k})^{i}|\right\\},\rho\right\\}$

so that $b_{k}^{i}=\frac{\overline{y}_{k-1}^{i}}{s_{k-1}^{i}}>0.$

In a situation where $s_{k-1}^{i}=0,$ then $a_{k-1}^{i}$ and $b_{k-1}^{i}$
will assume any suitable nonnegative safeguarding value.

###  2.2. Algorithm

In this subsection, we present the proposed algorithm. Let $d_{k}$ and $g_{k}$
denote the search direction and the gradient of the objective function (1.1)
respectively. The search direction $d_{k}$, is obtained by solving the linear
systems

$H_{k}d_{k}=-g_{k},$

where

(2.16) $H_{k}=\begin{cases}I,&\text{if}~{}~{}k=0\\\
diag(a_{k}+b_{k}),&\text{if}~{}~{}k\geq 1\end{cases}$

is a diagonal matrix whose entries are computed by

(2.17)
$a_{k}^{i}+b_{k}^{i}=\begin{cases}\frac{\hat{y}_{k-1}^{i}+\overline{y}_{k-1}^{i}}{s_{k-1}^{i}}&\text{if}~{}~{}s_{k-1}^{i}\neq
0\\\ 1&\text{if}~{}~{}s_{k-1}^{i}=0,\end{cases}i=1,\cdots,n.$

The vectors $\hat{y}_{k-1}$ and $\overline{y}_{k-1}$ are defined by (2.7) with
some of their components possibly redefined by (2.12)$-$(2.15). Furthermore,
we safeguard the diagonal entries of the diagonal Hessian $H_{k}$ from
assuming extremely small and extremely large values by means of projecting
them into a given scalar interval $[l,u],$ such that $0<l\leq 1\leq
u<<+\infty.$ Hence, the $i$th diagonal entry of our Hessian matrix $H_{k}$ in
which $y_{k-1}^{i}=\hat{y}_{k-1}^{i}+\overline{y}_{k-1}^{i}$ for each
$i=1,2,...,n$ is given by

(2.18)
${h_{k}^{i}}=\begin{cases}\frac{y_{k-1}^{i}}{s_{k-1}^{i}},&\text{if}~{}~{}l\leq\frac{y_{k-1}^{i}}{s_{k-1}^{i}}\leq
u,\\\ l,&\text{if}~{}~{}\frac{y_{k-1}^{i}}{s_{k-1}^{i}}<l,\\\
u,&\text{if}~{}~{}\frac{y_{k-1}^{i}}{s_{k-1}^{i}}>u,\\\
1,&\text{if}~{}~{}s_{k-1}^{i}=0\end{cases}$

Given a starting point $x_{0},$ we compute the next iterate via

(2.19) $x_{k+1}=x_{k}+\alpha_{k}d_{k},~{}~{}k=0,1,2,\cdots.$

Here the search direction is given by

(2.20) $d_{k}=H_{k}^{-1}g_{k},$

$H_{k}$ is a structured diagonal Hessian with diagonal entries $h_{k}^{i}$,
the gradient $g_{k}=J_{k}^{T}F_{k}$ and $J_{k}$ and $F_{k}$ are the Jacobian
matrix and function evaluation at $x_{k}$ respectively.
We adopt the non-monotone line search proposed by Zhang and Hager [27] to
determine the step length $\alpha_{k}$. Let the search direction $d_{k}$
defined by (2.20) be a descent direction, then the step length $\alpha_{k}>0$
in (2.19) should satisfy the following non-monotone Armijo-type line search
technique

(2.21) $f(x_{k}+\alpha_{k}d_{k})\leq P_{k}+\theta\alpha_{k}g_{k}^{T}d_{k},$

where

(2.22) $\begin{cases}P_{0}=f(x_{0})\\\
P_{k+1}=\frac{\eta_{k}Q_{k}P_{k}+f(x_{k+1})}{Q_{k+1}}\\\ Q_{0}=1,\\\
Q_{k+1}=\eta_{k}Q_{k}+1,\end{cases}$

and $\theta\in(0,1),$ $\eta_{k}\in[0,1].$
Next we give the following remarks:

  * Remark A
  * (i)

Note that the $P_{k+1}$ in the above line search technique is a convex
combination of $P_{k}$ and $f(x_{k+1}).$ Since $P_{0}=f(x_{0})$, it follows
that the sequence $\\{P_{k}\\}$ is a convex combination of the function values
$f(x_{i}),$ for $i=0,1,2,\cdots,k.$

  * (ii)

The parameter $\eta_{k}$ controls the degree of monotonicity. If for each $k,$
$\eta_{k}=0,$ then the line search (2.21) is the usual monotone (Armijo-type);
otherwise, it is non-monotone.

  * (iii)

If for each $k,$ $\eta_{k}=1,$ then $P_{k}=\psi_{k}$ where

(2.23) $\psi_{k}=\frac{1}{k+1}\sum_{i=1}^{k}f(x_{i})$

We now formally state the steps of our proposed iterative algorithm with
structured diagonal Hessian.

  * Algorithm 1: Algorithm with Structured Diagonal Hessian (ASDH)

  * Step 0.

Given $x_{0}\in\mathbb{R}^{n}$, $\gamma,\theta\in(0,1)$,
$0\leq\eta_{\min}\leq\eta_{\max}\leq 1,$ $0<l\leq 1\leq u,$ $\rho,\epsilon>0$,
and $k_{\max}\in\mathbb{N}.$

  * Step 1.

Set $k=0,$ $H_{k}=I,$ $Q_{k}=1$. Compute $F_{k}$ and $g_{k},$ and set
$P_{k}=f_{k}.$

  * Step 2.

If $\|g_{k}\|\leq\epsilon$ and $k\geq k_{\max},$ stop. Else compute $d_{k}$
using (2.20).

  * Step 3.

Perform nonmonotone line search

    * Step 3.1.

Set $\alpha=1,$

    * Step 3.2.

if the following inequality

$f(x_{k}+\alpha d_{k})\leq C_{k}+\alpha\theta g_{k}^{T}d_{k}.$

holds, then proceed to Step 4.
Else, set $\alpha_{k}=\alpha/2$ and repeat Step 3.2

  * Step 4.

Set $\alpha_{k}=\alpha$ and compute the next iterate using (2.19).

  * Step 5.

Set $s_{k}=\alpha_{k}d_{k}$ and compute $\hat{y}_{k}$ and $\overline{y}_{k}$
using (2.7).

  * Step 6.

Safeguard $\hat{y}_{k}$ and $\overline{y}_{k}$ using (2.12)$-$(2.15) in case
any of them has different sign with $s_{k}^{i}$.

  * Step 7.

Update the diagonal Hessian $H_{k+1}$ using (2.18).

  * Step 8.

Choose $\eta_{k}\in[\eta_{\min},\eta_{\max}]$ and compute $Q_{k+1}$ and
$P_{k+1}$ using (2.22).

  * Step 9.

Set $k=k+1$ and go to step 2.

  * Remark B
  * (i)

Though the vectors $g_{k},$ $\hat{y}_{k}$ and $\overline{y}_{k}$ are in the
form of matrix-vector products, these products were obtained by writing a
MATLAB code that computes them directly without forming or storing the
Jacobian matrix.

  * (ii)

Since the Hessian $H_{k}$ is a diagonal matrix, its inverse $H_{k}^{-1}$ is
obtained by taking the reciprocal of each $h_{i}$ $(i=1,2,\cdots,n)$ using
Equation (2.18) taking into account our safeguarding rule. Therefore, the
product $H_{k}^{-1}g_{k}$ is simply component-wise vector multiplications.

From the discussions in Remark B above, we can see that our proposed method is
matrix-free and therefore suitable for large-scale problems.

##  3\. Convergence Analysis

In this section, we discuss the global convergence of our proposed method. We
begin by stating the following assumptions which will be useful in our
analysis.

  * A1.

The level set $\mathcal{D}=\\{x\in\mathbb{R}^{n}~{}|~{}f(x)\leq f(x_{0})\\}$
is bounded, i.e. there exists a positive constant $\omega$ such that
$\|x\|\leq\omega$ for all $x\in\mathcal{D}.$

  * A2.

There exist constants $L_{1}$ and $L_{2}$ such that for all
$x,y\in\mathcal{D},$ we have

(3.1) $\|J(x)-J(y)\|\leq L_{1}\|x-y\|$

and

(3.2) $\|F(x)-F(y)\|\leq L_{2}\|x-y\|$

From (3.1) and (3.2), we can obtain the followings

$\|g(x)-g(y)\|\leq l\|x-y\|,$
$\|F(x)\|\leq\omega_{1},~{}~{}\|J(x)\|\leq\omega_{2},~{}~{}\|g(x)\|\leq\gamma_{3},$

where $l,$ $\omega_{1},$ $\omega_{2}$ and $\omega_{3}$ are positive constants.

######  Lemma 3.1.

Let the sequence of search directions $\\{d_{k}\\}$ be generated by the ASDH
algorithm, then there exist $m_{1},m_{2}$ positive constants such that for all
$k=0,1,2,\cdots,$ the following relations hold

(3.3) $g_{k}^{T}d_{k}\leq-m_{1}\|g_{k}\|^{2},$

and

(3.4) $\|d_{k}\|\leq m_{2}\|g_{k}\|.$

###### Proof.

From the definition of $H_{k}^{-1},$ the diagonal entries defined by (2.18) is
bounded for each $i$ and for all $k,$ i.e. $l\leq h_{k}^{i}\leq u.$
Now, by the Equation (2.20) we have

$\begin{split}g_{k}^{T}d_{k}&=-g_{k}^{T}H_{k}^{-1}g_{k}\\\
&=-\sum_{i=1}^{n}(g_{k}^{i})^{2}/h_{k}^{i}\\\
&\leq-(1/u)\sum_{i=1}^{n}(g_{k}^{i})^{2}\\\ &=-(1/u)\|g_{k}\|^{2}.\end{split}$

If we let $m_{1}=1/u,$ then (3.3) holds.
In a similar way, since the diagonal matrix is always symmetric, we get

$\begin{split}\|d_{k}\|^{2}&=g_{k}^{T}H_{k}^{-2}g_{k}\\\
&=\sum_{i=1}^{n}(g_{k}^{i}/h_{k}^{i})^{2}\\\
&\leq(1/l^{2})\sum_{i=1}^{n}(g_{k}^{i})^{2}\\\
&=(1/l^{2})\|g_{k}\|^{2}.\end{split}$

If we let $m_{2}=1/l,$ we obtain (3.4) and the proof is complete. ∎

Equations (3.3), (3.4) and the Proposition 1 in [15] imply tha the ASDH
Algorithm is well-defined. Next, we state the following result which comes
from Lemma 1.1 in [27].

######  Lemma 3.2.

The iterates generated by the ASDH Algorithm satisfy $f(x_{k})\leq
P_{k}\leq\psi_{k},$ for all $k\geq 0,$ where $P_{k}$ and $\psi_{k}$ are
defined by (2.22) and (2.23) respectively.

######  Lemma 3.3.

The sequence of iterates $\\{x_{k}\\}$ generated by the ASDH Algorithm satisfy
$f(x_{k})\leq f(x_{0}),$ for each $k\geq 0.$

###### Proof.

By substituting $\eta_{k}Q_{k}=Q_{k+1}-1,$ in $P_{k+1}$ defined in Equation
(2.22) we have

$\begin{split}P_{k+1}&=\frac{(Q_{k+1}-1)P_{k}+f(x_{k+1})}{Q_{k+1}},\\\
&\leq\frac{Q_{k+1}P_{k}-P_{k}+P_{k}+\alpha_{k}\theta
m_{1}g_{k}^{T}d_{k}}{Q_{k+1}},\\\ &\leq\frac{Q_{k+1}P_{k}-\alpha_{k}\theta
m_{1}\|g_{k}\|^{2}}{Q_{k+1}},\\\ &=P_{k}-\frac{\alpha_{k}\theta
m_{1}\|g_{k}\|^{2}}{Q_{k+1}}\\\ &\leq P_{k}.\end{split}$

The first two inequalities respectively come from (2.21) and (3.3). Now, since
$f_{k}\leq P_{k},$ (Lemma 3.2), we obtain

$f(x_{k+1})\leq P_{k+1}\leq P_{k}\leq P_{k-1}\leq\cdots\leq P_{0}=f(x_{0}),$

so that $\\{x_{k}\\}\subset\mathcal{D}$ and the proof is complete. ∎

######  Theorem 3.4.

Let $f(x)$ be defined by problem (1.1) and suppose assumptions A1 and A2 hold.
Then, the sequence of iterates $\\{x_{k}\\}$ generated by the ASDH Algorithm
is contained in the level set $\mathcal{D}$ and

(3.5) $\liminf_{k\to\infty}\|g_{k}\|=0.$

Moreover, if $\eta_{\max}<1$ then

(3.6) $\lim_{k\to\infty}\|g_{k}\|=0.$

###### Proof.

The proof of this theorem follows directly from [27]. ∎

##  4\. Numerical Experiments

In this section, we turn our attention to numerical experiments to assess the
performance of our proposed ASDH method compared to SDHAM method proposed in
[15]. The SDHAM method is also a matrix-free algorithm for nonlinear least
squares problems which exploits the special structure of the Hessian of the
objective function; and generates its search direction using diagonal Hessian
approximation similar to our proposed ASDH algorithm.

In our experiment, we solved 30 test problems of which 22 are large scale and
8 are small scale (see Table 1). We vary the dimensions of the large scale
problems as 1000, 5000; 10,000. All the test problems considered are properly
cited. The parameters used in the experiment are as follows

  * •

ASDH algorithm: $\gamma=0.2,$
$\eta_{k}=0.75e^{(-(k/45)^{2})}+0.1~{}\text{with}~{}\eta_{\min}=0.1,$
$\eta_{\max}=0.85,$ $l=10^{-30},$ $u=10^{30};$ and $\rho=0.0001$

  * •

SDHAM algorithm: All parameters are as presented in [15].

All codes were written in MATLAB R2017a and run on a PC with intel COREi5
processor with 4GB of RAM and CPU 2.3GHZ speed. The iteration is terminated
whenever the inequality $\|g_{k}\|\leq 10^{-4}$ is satisfied. Failure, denoted
by F, is recorded when the number of iterations exceeds 1,000 and the stopping
criterion mentioned above has not been satisfied.

In Tables 2$-$5, we report the results of the following information: the
number of iterations (NITER) needed by each solver to converge to an
approximate solution, the number of function evaluation (NFVAL), the number of
matrix-vector product (NMVP) the CPU time in seconds (TIME), and the objective
function $f$ value at the minimizer (FVALUE). The test problems are denoted by
P$i,$ $i=1,2,\cdots,30.$

In the Figures 1$-$4, we adopt the popular performance profile by Dolan and
Moré [5] to compare the performance of the ASDH method with that of SDHAM
method based on the number of iterations, number of functions evaluation,
number of matrix-vector product and CPU time. Though the two methods are
competitive, it can be seen from the Figures 1$-$4 that all the curves with
respect to our proposed ASDH method stay longer on the vertical axis which
means it solves more problems with less NITER, NFVAL, NMVP and TIME compared
to the SDHAM method. Specifically, it can observed from Figures 1$-$3 that our
method solves about $80\%$ of the test problems with least NITER, NFVAL, NMVP.
Also, from Figure 4 we can see that our method solves about $70\%$ of the test
problems with least CPU TIME.
Moreover, from the information reported in Table 2$-$5, it can be seen that
ASDH method solves all the test problems without any failure while the SDHAM
recorded 3 failures. In the overall experiments, ASDH method needs less number
of iterations, number of functions evaluation, number of matrix-vector product
and CPU time to obtain the minimizer of most of the test problems compared to
the SDHAM methods.

Table 1. List of test problems with references and starting points Problems | Function name | Starting point
---|---|---
Large scale
P1 | Penalty function I [12] |  $(1/3,1/3,\cdots,1/3)^{T}$ [t]
P2 | Trigonometric function [19] | $(1/n,\cdots,1/n)^{T}$
P3 | Discrete boundary value [19] | $(\frac{1}{n+1}(\frac{1}{n+1}-1),\cdots,\frac{1}{n+1}(\frac{n}{n+1}-1))^{T}$
P4 | Linear function full rank [19] | $(1,1,\cdots,1)^{T}$
P5 | Problem 202 [14] | $(2,2,\cdots,2)^{T}$
P6 | Problem 206 [14] | $(1/n,\cdots,1/n)^{T}$
P7 | Problem 212 [14] | $(0.5,\cdots,0.5)^{T}$
P8 | Strictly convex function I [20] | $(1/n,2/n,\cdots,1)^{T}$
P9 | Strictly convex function II [20] | $(1,1,\cdots,1)^{T}$
P10 | Brown almost linear [19] | $(0.5,\cdots,0.5)^{T}$
P11 | Exponential function I [12] | $(\frac{n}{n-1},\cdots,\frac{n}{n-1})^{T}$
P12 | Singular function [12] | $(1,1,\cdots,1)^{T}$
P13 | Logarithmic function [12] | $(1,1,\cdots,1)^{T}$
P14 | Extended Freudenstein and Roth [12] | $(6,3,6,3,\cdots,6,3)^{T}$
P15 | Extended Powell singular [12] | $(1.5E-4,\cdots,1.5E-4)^{T}$
P16 | Function21 [12] | $(-1,-1,\cdots,-1)^{T}$
P17 | Broyden tridiagonal function [19] | $(-1,-1,\cdots,-1)^{T}$
P18 | Generalized Broyden tridiagonal [14] | $(-1,-1,\cdots,-1)^{T}$
P19 | Extended Rosenbrock [19] | $(-1,1,-1,1,\cdots,-1,1)^{T}$
P20 | Extended Himmelblau [18] | $(1,1/n,1,1/n,\cdots,1,1/n)^{T}$
P21 | Function 27 [12] | $(100,1/n^{2},1/n^{2},\cdots,1/n^{2})^{T}$
P22 | Trigonometric logarithmic function [15] |  $(1,1,\cdots,1)^{T}$ [b]
Small scale |  |
P23 | Bard function [19] |  $(-1000,-1000,-1000)^{T}$ [t]
P24 | Brown badly scaled [19] | $(1,1)^{T}$
P25 | Jennrich and Sampson [19] | $(0.2,0.2)^{T}$
P26 | Box 3D function [19] | $(0,0.1)^{T}$
P27 | Rank deficient Jacobian [8] | $(-1,1)^{T}$
P28 | Rosenbrock function [19] | $(-1,1)^{T}$
P29 | Parameterized problem [9] | $(10,10)^{T}$
P30 | Freudenstein and Roth function [19] |  $(0.5,-2)^{T}$ [b]
Table 2. Numerical results of our ASDH and SDHAM methods for large scale
problems $1-22$ with dimension $n=m=1,000$

| ASDH | SDHAM
---|---|---
Problem | NITER | NFEVAL | NMVP | TIME | FVALUE | ITER | NFEVAL | NMVP | TIME | FVALUE
P1 | 6 | 7 | 19 | 0.022139 | 2.36E-06 | 5 | 6 | 16 | 0.059823 | 1.23E-05
P2 | 193 | 354 | 580 | 0.26363 | 1.34E-07 | 5 | 12 | 16 | 0.011127 | 1.89E-15
P3 | 18 | 20 | 55 | 0.063171 | 2.52E-08 | 20 | 23 | 61 | 0.074233 | 2.33E-08
P4 | 2 | 3 | 7 | 0.02106 | 0.502 | 2 | 3 | 7 | 0.020876 | 0.502
P5 | 5 | 6 | 16 | 0.013713 | 7.22E-13 | 5 | 6 | 16 | 0.017462 | 1.39E-10
P6 | 43 | 49 | 130 | 0.099136 | 3E-08 | 44 | 57 | 133 | 0.03711 | 2.96E-08
P7 | 5 | 6 | 16 | 0.012173 | 9.06E-12 | 7 | 8 | 22 | 0.015851 | 3.32E-12
P8 | 5 | 6 | 16 | 0.007124 | 500 | 5 | 6 | 16 | 0.01575 | 500
P9 | 8 | 23 | 25 | 0.039868 | 1669168 | 8 | 23 | 25 | 0.012976 | 1669168
P10 | 2 | 23 | 7 | 0.005995 | 1.99E-13 | F | F | F | F | F
P11 | 4 | 5 | 13 | 0.003988 | 9.44E-08 | 5 | 6 | 16 | 0.008688 | 4.53E-08
P12 | 22 | 41 | 67 | 0.061689 | 55.0298 | 25 | 44 | 76 | 0.081369 | 55.0298
P13 | 6 | 7 | 19 | 0.025038 | 8.45E-17 | 6 | 7 | 19 | 0.012527 | 1.57E-12
P14 | 19 | 27 | 58 | 0.038295 | 8.12E-12 | 21 | 30 | 64 | 0.018482 | 6.24E-14
P15 | 1 | 7 | 4 | 0.018509 | 1.21E-12 | 1 | 7 | 4 | 0.009089 | 1.21E-12
P16 | 85 | 116 | 256 | 0.13824 | 6.2E-10 | 75 | 101 | 226 | 0.21378 | 2.51E-10
P17 | 247 | 507 | 742 | 0.26142 | 0.35626 | 19 | 34 | 58 | 0.10981 | 1.63E-11
P18 | 72 | 125 | 217 | 0.1084 | 0.76757 | 103 | 148 | 310 | 0.13672 | 0.026797
P19 | 1 | 2 | 4 | 0.009725 | 0 | 1 | 2 | 4 | 0.017676 | 0
P20 | 17 | 22 | 52 | 0.018535 | 7.54E-13 | 30 | 34 | 91 | 0.088546 | 5.78E-11
P21 | 17 | 32 | 52 | 0.022589 | 6.37E-07 | 18 | 33 | 55 | 0.064815 | 2.16E-07
P22 | 6 | 7 | 19 | 0.02308 | 8.16E-17 | 6 | 9 | 19 | 0.29798 | 1.03E-14

Table 3. Numerical results of our ASDH and SDHAM methods for large scale
problems $1-22$ with dimension $n=m=5,000$

| ASDH | SDHAM
---|---|---
Problem | NITER | NFEVAL | NMVP | TIME | FVALUE | ITER | NFEVAL | NMVP | TIME | FVALUE
P1 | 15 | 16 | 46 | 0.066157 | 7.93E-05 | 15 | 16 | 46 | 0.17645 | 7.92E-05 [t]
P2 | 120 | 273 | 361 | 0.76396 | 7.91E-08 | F | F | F | F | F
P3 | 3 | 5 | 10 | 0.044907 | 1.44E-09 | 3 | 5 | 10 | 0.16634 | 1.44E-09
P4 | 2 | 3 | 7 | 0.027567 | 0.5004 | 2 | 3 | 7 | 0.09278 | 0.5004
P5 | 5 | 6 | 16 | 0.029966 | 7.23E-13 | 5 | 6 | 16 | 0.052194 | 6.94E-10
P6 | 7 | 9 | 22 | 0.060726 | 2.32E-09 | 6 | 8 | 19 | 0.11743 | 2.83E-09
P7 | 5 | 6 | 16 | 0.018999 | 9.11E-12 | 7 | 8 | 22 | 0.1734 | 3.32E-12
P8 | 5 | 6 | 16 | 0.075254 | 2500 | 5 | 6 | 16 | 0.047758 | 2500
P9 | 8 | 27 | 25 | 0.071887 | 2.08E+08 | 9 | 28 | 28 | 0.13815 | 2.08E+08
P10 | 3 | 29 | 10 | 0.027006 | 5.17E-16 | 2 | 28 | 7 | 0.072295 | NaN
P11 | 3 | 4 | 10 | 0.047036 | 9.54E-08 | 4 | 5 | 13 | 0.046785 | 3.62E-08
P12 | 27 | 50 | 82 | 0.25832 | 1386.256 | 29 | 52 | 88 | 0.26293 | 1386.255
P13 | 6 | 7 | 19 | 0.040228 | 3.93E-16 | 6 | 7 | 19 | 0.062598 | 7.6E-12
P14 | 21 | 29 | 64 | 0.062354 | 1.41E-13 | 21 | 30 | 64 | 0.1584 | 3.12E-13
P15 | 1 | 7 | 4 | 0.024886 | 6.03E-12 | 1 | 7 | 4 | 0.12042 | 6.03E-12
P16 | 69 | 96 | 208 | 0.36005 | 5.97E-10 | 83 | 113 | 250 | 0.62034 | 4.57E-10
P17 | 140 | 305 | 421 | 0.32061 | 2.26E-10 | 19 | 34 | 58 | 0.20367 | 1.63E-11
P18 | 58 | 98 | 175 | 0.21311 | 2.26E-10 | 97 | 132 | 292 | 0.3257 | 0.031054
P19 | 1 | 2 | 4 | 0.005155 | 0 | 1 | 2 | 4 | 0.088945 | 0
P20 | 40 | 50 | 121 | 0.21105 | 1.15E-10 | 30 | 35 | 91 | 0.1997 | 3.18E-11
P21 | 17 | 32 | 52 | 0.10225 | 6.37E-07 | 18 | 33 | 55 | 0.10118 | 2.16E-07
P22 | 6 | 7 | 19 | 0.048698 | 3.9E-16 | 6 | 9 | 19 | 0.56375 | 3.33E-12 [b]

Table 4. Numerical results of our ASDH and SDHAM methods for large scale
problems $1-22$ with dimension $n=m=10,000$

| ASDH | SDHAM
---|---|---
Problem | NITER | NFEVAL | NMVP | TIME | FVALUE | ITER | NFEVAL | NMVP | TIME | FVALUE
P1 | 23 | 24 | 70 | 0.49096 | 0.000142 | 23 | 24 | 70 | 0.13661 | 0.000142 [t]
P2 | 81 | 192 | 244 | 0.67185 | 7.36E-08 | F | F | F | F | F
P3 | 2 | 4 | 7 | 0.095771 | 1.58E-09 | 2 | 4 | 7 | 0.044607 | 1.58E-09
P4 | 2 | 3 | 7 | 0.055015 | 0.5002 | 2 | 3 | 7 | 0.030105 | 0.5002
P5 | 5 | 6 | 16 | 0.06105 | 7.25E-13 | 5 | 6 | 16 | 0.040306 | 1.39E-09
P6 | 3 | 5 | 10 | 0.020898 | 7.16E-10 | 3 | 5 | 10 | 0.025454 | 7.16E-10
P7 | 5 | 6 | 16 | 0.040219 | 9.17E-12 | 7 | 8 | 22 | 0.056896 | 3.32E-12
P8 | 5 | 6 | 16 | 0.1955 | 5000 | 5 | 6 | 16 | 0.051594 | 5000
P9 | 9 | 30 | 28 | 0.099793 | 1.67E+09 | 10 | 31 | 31 | 0.15435 | 1.67E+09
P10 | 3 | 4 | 10 | 0.019368 | 4.77E-08 | 3 | 4 | 10 | 0.021593 | 7.24E-08
P11 | 199 | 3849 | 598 | 4.6179 | 2.49E-15 | 2 | 19 | 7 | 0.025921 | 1.13E-13
P12 | 28 | 53 | 85 | 0.27666 | 5550.287 | 30 | 55 | 91 | 0.51292 | 5550.286
P13 | 6 | 7 | 19 | 0.066932 | 7.79E-16 | 6 | 7 | 19 | 0.03482 | 1.51E-11
P14 | 21 | 29 | 64 | 0.11086 | 2.83E-13 | 21 | 30 | 64 | 0.16018 | 6.24E-13
P15 | 1 | 7 | 4 | 0.011959 | 1.21E-11 | 1 | 7 | 4 | 0.014576 | 1.21E-11
P16 | 89 | 122 | 268 | 0.90064 | 5.94E-10 | 86 | 112 | 259 | 0.84402 | 8.47E-10
P17 | 286 | 620 | 859 | 1.5625 | 0.35626 | 19 | 34 | 58 | 0.20675 | 1.63E-11
P18 | 56 | 80 | 169 | 0.28569 | 3.2E-10 | 238 | 362 | 715 | 1.3739 | 0.010024
P19 | 1 | 2 | 4 | 0.016878 | 0 | 1 | 2 | 4 | 0.018502 | 0
P20 | 24 | 29 | 73 | 0.14056 | 6.17E-11 | 26 | 31 | 79 | 0.23926 | 4.71E-11
P21 | 17 | 32 | 52 | 0.11498 | 6.37E-07 | 18 | 33 | 55 | 0.17027 | 2.16E-07
P22 | 6 | 7 | 19 | 0.031825 | 7.77E-16 | 6 | 9 | 19 | 0.12546 | 8.61E-12 [b]

Table 5. Numerical results of our ASDH and SDHAM methods for small scale
problems $23-30$

|  | ASDH | SDHAM
---|---|---|---
Problem | DIM | ITER | NFEVAL | NMVP | TIME | FVALUE | ITER | NFEVAL | NMVP | TIME | FVALUE
P23 | n=3, m=15 | 1 | 4 | 4 | 0.026898 | 8.7448 | 1 | 4 | 4 | 0.16203 | 8.7448 [t]
P24 | n=2, m=3 | 9 | 12 | 28 | 0.010317 | 7.16E-14 | 10 | 12 | 31 | 0.060814 | 1.05E-27
P25 | n=2, m=20 | 4 | 20 | 13 | 0.00827 | 62.1811 | 5 | 21 | 16 | 0.046634 | 62.1811
P26 | n=3, m=10 | 47 | 64 | 142 | 0.063188 | 4.06E-09 | 38 | 50 | 115 | 0.12071 | 1.83E-10
P27 | n=2, m=3 | 6 | 9 | 19 | 0.008934 | 1.2905 | 6 | 9 | 19 | 0.077812 | 1.2905
P28 | n=m=2 | 1 | 2 | 4 | 0.005085 | 0 | 1 | 2 | 4 | 0.040352 | 0
P29 | n=2, m=3 | 6 | 22 | 19 | 0.018706 | 0.49999 | 6 | 22 | 19 | 0.039052 | 0.49999
P30 | n=m=2 | 81 | 113 | 244 | 0.036083 | 24.4921 | 93 | 133 | 280 | 0.10623 | 24.4921 [b]

Figure 1. Performance profile with respect to number of iterations  Figure 2.
Performance profile with respect to number of function evaluations  Figure 3.
Performance profile with respect to number of matrix-vector product Figure 4.
Performance profile with respect to CPU time

##  5\. Conclusions

We have proposed an iterative algorithm with structured diagonal Hessian
approximation for solving nonlinear least square problems. The proposed
algorithm neither forms nor stores matrices which make it suitable for large
scale problems. We have devised appropriate safeguards to ensure the search
directions generated by our proposed algorithm are descent. The proposed ASDH
method was developed using a diagonal Hessian approximation that contains more
information of the objective function than the one proposed by Mohammad and
Santos [15]. We presented some preliminary numerical experiments to show the
robustness of our ASDH method.

Acknowledgements: The first and second authors acknowledge the financial
support provided by King Mongkut’s University of Technology Thonburi through
the ”KMUTT 55th Anniversary Commemorative Fund”. The first author was
supported by the Petchra Pra Jom Klao Doctoral Scholarship Academic for Ph.D.
Program at KMUTT. This project was supported by the Theoretical and
Computational Science (TaCS) Center under Computational and Applied Science
for Smart Innovation Research Cluster (CLASSIC), Faculty of Science, KMUTT.

## References

  * [1] Tilman Barz, Stefan Körkel, Günter Wozny, et al.  Nonlinear ill-posed problem analysis in model-based parameter estimation and experimental design.  Computers & Chemical Engineering, 77:24–42, 2015.
  * [2] Anastasia Cornelio.  Regularized nonlinear least squares methods for hit position reconstruction in small gamma cameras.  Applied Mathematics and Computation, 217(12):5589–5595, 2011.
  * [3] R Dehghani and N Mahdavi-Amiri.  Scaled nonlinear conjugate gradient methods for nonlinear least squares problems.  Numerical Algorithms, pages 1–20, 2018.
  * [4] JE Dennis, Héctor J Martinez, and Richard A Tapia.  Convergence theory for the structured bfgs secant method with an application to nonlinear least squares.  Journal of Optimization Theory and Applications, 61(2):161–178, 1989.
  * [5] Elizabeth D Dolan and Jorge J Moré.  Benchmarking optimization software with performance profiles.  Mathematical programming, 91(2):201–213, 2002.
  * [6] R Fletcher and C Xu.  Hybrid methods for nonlinear least squares.  IMA Journal of Numerical Analysis, 7(3):371–389, 1987.
  * [7] Gene Golub and Victor Pereyra.  Separable nonlinear least squares: the variable projection method and its applications.  Inverse problems, 19(2):R1, 2003.
  * [8] Douglas S Gonçalves and Sandra A Santos.  Local analysis of a spectral correction for the gauss-newton model applied to quadratic residual problems.  Numerical Algorithms, 73(2):407–431, 2016.
  * [9] Jürgen Huschens.  On the use of product structure in secant methods for nonlinear least squares problems.  SIAM Journal on Optimization, 4(1):108–129, 1994.
  * [10] Seung-Jean Kim, Kwangmoo Koh, Michael Lustig, Stephen Boyd, and Dimitry Gorinevsky.  An interior-point method for large-scale $l_{1}-$regularized least squares.  IEEE journal of selected topics in signal processing, 1(4):606–617, 2007.
  * [11] Michiya Kobayashi, Yasushi Narushima, and Hiroshi Yabe.  Nonlinear conjugate gradient methods with structured secant condition for nonlinear least squares problems.  Journal of computational and applied mathematics, 234(2):375–397, 2010.
  * [12] William La Cruz, José Mario Martínez, and Marcos Raydan.  Spectral residual method without gradient information for solving large-scale nonlinear systems: theory and experiments.  $http://kuainasi.ciens.ucv.ve/mraydan/download_{p}apers/TechRep.pdf$, 2004.
  * [13] Junhong Li, Feng Ding, and Guowei Yang.  Maximum likelihood least squares identification method for input nonlinear finite impulse response moving average systems.  Mathematical and Computer Modelling, 55(3-4):442–450, 2012.
  * [14] Ladislav Lukšan and Jan Vlcek.  Test problems for unconstrained optimization.  Academy of Sciences of the Czech Republic, Institute of Computer Science, Technical Report, (897), 2003.
  * [15] Hassan Mohammad and Sandra A Santos.  A structured diagonal hessian approximation method with evaluation complexity analysis for nonlinear least squares.  Computational and Applied Mathematics, 37(5):6619–6653, 2018.
  * [16] Hassan Mohammad and Mohammed Yusuf Waziri.  Structured two-point stepsize gradient methods for nonlinear least squares.  Journal of Optimization Theory and Applications, 181(1):298–317, 2019.
  * [17] Hassan Mohammad, Mohammed Yusuf Waziri, and Sandra Augusta Santos.  A brief survey of methods for solving nonlinear least-squares problems.  Numerical Algebra, Control & Optimization, 9(1):1–13, 2019.
  * [18] Jamil Momin and Yang Xin-She.  A literature survey of benchmark functions for global optimization problems.  Journal of Mathematical Modelling and Numerical Optimisation, 4(2):150–194, 2013.
  * [19] Jorge J Moré, Burton S Garbow, and Kenneth E Hillstrom.  Testing unconstrained optimization software.  Technical report, Argonne National Lab., IL (USA), 1978.
  * [20] Marcos Raydan.  The barzilai and borwein gradient method for the large scale unconstrained minimization problem.  SIAM Journal on Optimization, 7(1):26–33, 1997.
  * [21] Li Min Tang.  A regularization homotopy iterative method for il–posed nonlinear least squares problem and its application.  In Applied Mechanics and Materials, volume 90, pages 3268–3273. Trans Tech Publ, 2011.
  * [22] CX Xu.  Hybrid method for nonlinear least-square problems without calculating derivatives.  Journal of optimization Theory and Applications, 65(3):555–574, 1990.
  * [23] Ya-Xiang Yuan.  Subspace methods for large scale nonlinear equations and nonlinear least squares.  Optimization and Engineering, 10(2):207–218, 2009.
  * [24] Ya-Xiang Yuan.  Recent advances in numerical methods for nonlinear equations and nonlinear least squares.  Numerical algebra, control and optimization, 1(1):15–34, 2011.
  * [25] Hongchao Zhang and Andrew R Conn.  On the local convergence of a derivative-free algorithm for least-squares minimization.  Computational optimization and applications, 51(2):481–507, 2012\.
  * [26] Hongchao Zhang, Andrew R Conn, and Katya Scheinberg.  A derivative-free algorithm for least-squares minimization.  SIAM Journal on Optimization, 20(6):3555–3576, 2010.
  * [27] Hongchao Zhang and William W Hager.  A nonmonotone line search technique and its application to unconstrained optimization.  SIAM journal on Optimization, 14(4):1043–1056, 2004.


# Casimir effect with one large extra dimension

Andrea Erdas<EMAIL_ADDRESS>Department of Physics, Loyola University
Maryland, 4501 North Charles Street, Baltimore, Maryland 21210, USA

###### Abstract

In this work I study the Casimir effect of a massive complex scalar field in
the presence of one large compactified extra dimension. I investigate the case
of a scalar field confined between two parallel plates in the macroscopic
three dimensions, and examine the cases of Dirichlet and mixed (Dirichlet-
Neumann) boundary conditions on the plates. The case of Neumann boundary
conditions is uninteresting, since it yields the same result as the case of
Dirichlet boundary conditions. The scalar field also permeates a fourth
compactified dimension of a size that could be comparable to the distance
between the plates. This investigation is carried out using the
$\zeta$-function regularization technique that allows me to obtain exact
expressions for the Casimir energy and pressure. I discover that, when the
compactified length of the extra dimension is similar to the plate distance,
or slightly larger, the Casimir energy and pressure become significantly
different than their standard three dimensional values, for either Dirichlet
or mixed boundary conditions. Therefore, the Casimir effect of a quantum field
that permeates a compactified fourth dimension, could be used as an effective
tool to explore the existence of large compactified extra dimensions.

###### pacs:

03.70.+k, 11.10.-z, 11.25.Mj, 12.20.Ds.

##  I Introduction

Hendrik Casimir first recognized that an attractive force exists between two
uncharged, conducting parallel plates in vacuum, entirely due to quantum
effects Casimir:1948dh . Sparnaay was first to attempt measuring this force
Sparnaay:1958wg , ten years after Casimir’s original paper. Many increasingly
accurate experiments followed throughout the decades Bordag:2001qi ;
Bordag:2009zz  that fully confirmed the Casimir effect, caused by quantum
fluctuations of the electromagnetic field in a vacuum constrained by two
conducting parallel plates. Fluctuations of quantum fields of different types
also produce Casimir forces, that depend strongly on the boundary conditions
at the plates of the field considered. For a scalar field, Dirichlet or
Neumann boundary conditions cause attraction between the plates, while mixed
(Dirichlet-Neumann) boundary conditions, cause a repulsive force Boyer:1974 .
Casimir forces depend strongly also on the geometrical shape of the conducting
plates, as discovered by Boyer Boyer:1968uf , when he realized that spherical
plates are subject to a repulsive force, caused by either the electromagnetic
field, or a scalar field, satisfying Dirichlet boundary conditions at the
spherical shell.

The existence of compactified extra dimensions has been proposed by Kaluza and
Klein a century ago, as an attempt to unify gravity and electromagnetism
Kaluza:1921tu ; Klein:1926tv . The idea has been revived twenty yers ago in a
series of famous papers by Arkani-Hamed, Dimopoulos and Dvali Arkani-
Hamed:1998jmv ; Antoniadis:1998ig ; Arkani-Hamed:1998sfv , where large compact
extra dimensions were proposed within a new framework, embedded in string
theory, for solving the hierarchy problem. The extra dimensions size proposed
by these authors should range between roughly a millimeter and a TeV-1. One
year later, Randall and Sundrum Randall:1999ee  found a different solution to
the hierarchy problem involving a five-dimensional AdS spacetime with warped
geometry. All these scenarios modify gravity by producing deviations from
Newton’s law of gravitation at the submillimiter distance. Several experiments
have been carried out in the last twenty years that test gravity at the
submillimiter scale, see Refs. Adelberger:2009zz ; Murata:2014nra  for a
comprehensive list of experiments, placing deviations from Newton’s gravity at
or below $40\mu$m. Cosmological and astrophysical constraints also exist on
the size and number of the compactified extra-dimensions, as well as collider
constraints from LHC.

The Casimir effect can be a powerful experimental tool to search for the
existence of a large compactified extra dimension, as long as an accurate
theoretical investigation is carried out to understand exactly the implication
of such extra dimension on the Casimir effect experiments in the lab. Several
papers have been published on this subject, see Refs. Teo:2008ah ; Teo:2009tm
for a comprehensive list, but the exact and detailed analysis that links
compactification length, plate distance and mass of the field to deviation of
Casimir energy and pressure from their well known three dimensional values, is
missing. In this paper I will undertake that detailed analysis by
investigating a massive scalar field constrained between two parallel plates
in three dimensions and satisfying either Dirichlet or mixed boundary
conditions on the plates. The scalar field also permeates a compactified
fourth dimension sized at the submillimiter scale. I will use a well known
regularization method introduced by Hawking Hawking:1976ja , the zeta function
technique, to regularize the Casimir energy Elizalde:1988rh ; Elizalde:2007du
.

In Sec. II of this work I describe the model of a massive scalar field
permeating a compactified extra-dimension and obtain the exact zeta function
in the two different cases of Dirichlet and mixed boundary conditions at the
parallel plates. Sec. III is dedicated to Dirichlet boundary conditions and I
obtain exact and simple expressions for the Casimir energy in the two cases of
a field with light mass and heavy mass respectively. In Sec. IV I investigate
mixed boundary conditions and find exact and simple expressions of the Casimir
energy for the two cases of light and heavy mass. In Sec. V I use the Casimir
energy obtained in the previous two sections to evaluate the Casimir pressure
for Dirichlet and mixed boundary conditions in the two cases of light and
heavy mass. A discussion of my results and my conclusions are in Sec. VI.

##  II The model and its zeta function

In this paper I investigate a complex scalar field
$\phi(x_{1},x_{2},x_{3},x_{4},t)$ of mass $M$ in ${\bf R}^{3}\times{\bf
S}^{1}$, where ${\bf S}^{1}$ is compactified into a circle of circumference
$b$. Two square parallel plates of area $L^{2}$ are present in ${\bf R}^{3}$
and are perpendicular to the $x_{3}$-axis and located at $x_{3}=0$ and
$x_{3}=a$. The following periodic boundary condition on $\phi$

$\phi(x_{1},x_{2},x_{3},0,t)=\phi(x_{1},x_{2},x_{3},b,t)$ (II.1)

implements the presence of the compactified extra-dimension ${\bf S}^{1}$
Elizalde:1988rh . In addition to the boundary condition of Eq. (II.1), I will
impose either Dirichlet boundary conditions

$\phi(x_{1},x_{2},0,x_{4},t)=\phi(x_{1},x_{2},a,x_{4},t)=0,$ (II.2)

or mixed boundary conditions

$\phi(x_{1},x_{2},0,x_{4},t)={\partial\phi\over\partial
x_{3}}(x_{1},x_{2},a,x_{4},t)=0,$ (II.3)

on the scalar field at the plates. I will not investigate Neumann boundary
conditions since they produce the same results as Dirichlet ones. My goal is
to understand how the presence of one compactified extra dimension modifies
the Casimir effect due to the complex scalar field, and will use the
generalized zeta function technique to investigate this problem.

I use Euclidean time $\tau$ and begin by writing the Klein Gordon operator

$D_{E}=-{\partial^{2}\over\partial\tau^{2}}-\nabla^{2}-{\partial^{2}\over\partial
x_{4}^{2}}+M^{2},$ (II.4)

where the contribution of the compactified extra dimension is shown
explicitly. The eigenvalues of $D_{E}$ are:

$\left\\{k_{0}^{2}+k_{1}^{2}+k_{2}^{2}+k_{3}^{2}+\left({2m\pi\over
b}\right)^{2}+M^{2}\right\\},$ (II.5)

where $k_{0},k_{1},k_{2}\in\Re$, $m=0,\pm 1,\pm 2,\cdots$, and $k_{3}$ takes
the following discrete values when considering Dirichlet boundary conditions

$k_{3}={n\pi\over a},$ (II.6)

with $n=1,2,\cdots$. When mixed boundary conditions are considered, $k_{3}$
takes different discrete values

$k_{3}=\left(n+\textstyle{\frac{1}{2}}\right){\pi\over a},$ (II.7)

with $n=0,1,2,\cdots$.

Using the eigenvalues of Eq. (II.5) and $k_{3}$ of Eq. (II.6), I construct the
zeta function for the case of Dirichlet boundary conditions

$\zeta(s)=\mu^{2s}{L^{2}\over(2\pi)^{3}}\int
d^{3}k\sum_{n=1}^{\infty}\sum_{m=-\infty}^{\infty}\left[k^{2}+\left({n\pi\over
a}\right)^{2}+\left({2m\pi\over b}\right)^{2}+M^{2}\right]^{-s},$ (II.8)

where $d^{3}k=dk_{0}dk_{1}dk_{2}$, $k^{2}=k_{0}^{2}+k_{1}^{2}+k_{2}^{2}$ and,
as it is done routinely when using the zeta function technique, I introduce
the multiplicative parameter $\mu$ with dimension of mass to keep the zeta
function with dimensions of energy for all values of $s$ Hawking:1976ja . For
mixed boundary conditions, I use Eqs. (II.5) and (II.7) and find

$\zeta(s)=\mu^{2s}{L^{2}\over(2\pi)^{3}}\int
d^{3}k\sum_{n=0}^{\infty}\sum_{m=-\infty}^{\infty}\left[k^{2}+\left(n+\textstyle{\frac{1}{2}}\right)^{2}\left({\pi\over
a}\right)^{2}+\left({2m\pi\over b}\right)^{2}+M^{2}\right]^{-s}.$ (II.9)

##  III Dirichlet boundary conditions

When investigating Dirichlet boundary conditions, I use the following identity

$z^{-s}=\frac{1}{\Gamma(s)}\int_{0}^{\infty}dt\,t^{s-1}e^{-zt}$ (III.1)

in Eq. (II.8), then do the straightforward $k$-integration to obtain

$\zeta(s)={\mu^{2s}\over\Gamma(s)}{L^{2}\over(2\sqrt{\pi})^{3}}\sum_{n=1}^{\infty}\sum_{m=-\infty}^{\infty}\int_{0}^{\infty}dt\,\,t^{s-5/2}e^{-({n^{2}\pi^{2}\over
a^{2}}+{4m^{2}\pi^{2}\over b^{2}}+M^{2})t}.$ (III.2)

At this point I do a Poisson resummation of the $n$-sum

$\sum_{n=1}^{\infty}e^{-{n^{2}\pi^{2}\over a^{2}}t}={a\over 2\sqrt{\pi
t}}\sum_{n=-\infty}^{\infty}e^{-{n^{2}a^{2}\over t}}-{\frac{1}{2}},$ (III.3)

drop the $-{1\over 2}$ term since it does not depend on the plate distance
$a$, and obtain

$\zeta(s)={\mu^{2s}\over\Gamma(s)}{L^{2}a\over(4\pi)^{2}}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-M^{2}t}\left(1+2\sum_{n=1}^{\infty}e^{-{n^{2}a^{2}\over
t}}\right)\left(1+2\sum_{m=1}^{\infty}e^{-{4m^{2}\pi^{2}\over
b^{2}}t}\right).$ (III.4)

Terms with only a linear dependence on $a$ can be ignored, since they are
uniform energy density terms that do not depend on the presence of the plates,
and I write

$\zeta(s)=\zeta_{1}(s)+\zeta_{2}(s),$ (III.5)

where

$\zeta_{1}(s)={\mu^{2s}\over\Gamma(s)}{L^{2}a\over
8\pi^{2}}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-M^{2}t}\sum_{n=1}^{\infty}e^{-{n^{2}a^{2}\over
t}},$ (III.6)

depends only on the plate distance $a$, and $\zeta_{2}(s)$ depends also on
$b$, the compactification length of the large extra dimension

$\zeta_{2}(s)={\mu^{2s}\over\Gamma(s)}{L^{2}a\over
4\pi^{2}}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-M^{2}t}\sum_{n=1}^{\infty}e^{-{n^{2}a^{2}\over
t}}\sum_{m=1}^{\infty}e^{-{4m^{2}\pi^{2}\over b^{2}}t}.$ (III.7)

I am going to investigate first the light mass scenario, $M<a^{-1}$, then the
heavy mass one, $M>a^{-1}$.

In the case of low mass, I do a power series expansion of the $e^{-M^{2}t}$
term in $\zeta_{1}(s)$ and retain the first three terms, change the
integration variable, ${n^{2}a^{2}\over t}=t^{\prime}$, to obtain

$\zeta_{1}(s)={(\mu a)^{2s}\over\Gamma(s)}{L^{2}\over
8\pi^{2}a^{3}}\left[\zeta_{R}(4-2s)\Gamma(2-s)-M^{2}a^{2}\zeta_{R}(2-2s)\Gamma(1-s)+{M^{4}a^{4}\over
2}\zeta_{R}(-2s)\Gamma(-s)\right],$ (III.8)

where $\zeta_{R}(s)$ is the Riemann zeta function and $\Gamma(s)$ is the Euler
gamma function. The Casimir energy $E$ is obtained from the zeta function by
taking its derivative at $s=0$

$E=-\zeta^{\prime}(0),$ (III.9)

therefore it is sufficient to evaluate $\zeta(s)$ when $s$ is near zero. For
$s$ close to zero I find

${(\mu a)^{2s}\over\Gamma(s)}\zeta_{R}(4-2s)\Gamma(2-s)\simeq{\pi^{4}\over
90}s+{\cal O}(s^{2}),$ (III.10) ${(\mu
a)^{2s}\over\Gamma(s)}\zeta_{R}(2-2s)\Gamma(1-s)\simeq{\pi^{2}\over 6}s+{\cal
O}(s^{2}),$ (III.11) ${(\mu
a)^{2s}\over\Gamma(s)}\zeta_{R}(-2s)\Gamma(-s)\simeq\left({1\over
2s}+\gamma_{E}+\log{\mu a\over 2\pi}\right)s+{\cal O}(s^{2}),$ (III.12)

where $\gamma_{E}$ is the Euler-Mascheroni constant, and obtain
$E_{1}=-\zeta^{\prime}_{1}(0)$, which is the contribution of $\zeta_{1}(s)$ to
the Casimir energy

$E_{1}=-{L^{2}\over 8\pi^{2}a^{3}}\left[{\pi^{4}\over 90}-{\pi^{2}\over
6}M^{2}a^{2}+{M^{4}a^{4}\over 2}\log(Ma)\right]$ (III.13)

where I neglected uniform energy density terms that do not contribute to the
Casimir energy and made the obvious choice $\mu=M$ indicating on-shell
renormalization. Next, I evaluate $\zeta_{2}(s)$ by changing integration
variable, $t\rightarrow{nab\over 2\pi m}t$, and obtain

$\zeta_{2}(s)=\sum_{m,n=1}^{\infty}\left({\mu^{2}nab\over 2\pi
m}\right)^{s}{L^{2}a\over 4\pi^{2}\Gamma(s)}\left({2\pi m\over
nab}\right)^{2}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-(M^{2}{nab\over 2\pi
m})t}e^{-{2\pi mna\over b}(t+t^{-1})}.$ (III.14)

As long as ${a\over b}>{1\over 2\pi}$, I can evaluate this integral using the
saddle point method and, for $s\ll 1$, obtain

$\zeta_{2}(s)\simeq{L^{2}\over\sqrt{2}(ab)^{3/2}}s\sum_{m,n=1}^{\infty}{m^{3/2}\over
n^{5/2}}e^{-(M^{2}{nab\over 2\pi m})}e^{-{4\pi mna\over b}},$ (III.15)

whose contribution to the Casimir energy, $E_{2}=-\zeta^{\prime}_{2}(0)$, is

$E_{2}=-{L^{2}\over\sqrt{2}(ab)^{3/2}}\sum_{m,n=1}^{\infty}{m^{3/2}\over
n^{5/2}}e^{-(M^{2}{nab\over 2\pi m})}e^{-{4\pi mna\over b}}.$ (III.16)

Below is the Casimir energy, $E=E_{1}+E_{2}$, per unit area, in the case of
light mass

${E\over L^{2}}=-{\pi^{2}\over
720a^{3}}\left[1-{15\over\pi^{2}}y^{2}+{45\over\pi^{4}}y^{4}\log
y+{720\over\sqrt{2}\pi^{2}}x^{3/2}\sum_{m,n=1}^{\infty}{m^{3/2}\over
n^{5/2}}e^{-{ny^{2}\over 2\pi mx}}e^{-{4\pi mnx}}\right],$ (III.17)

where the overall factor $-{\pi^{2}\over 720a^{3}}$ is the three dimensional
Casimir energy per unit area due to a massless scalar, and I introduced the
two dimensionless parameters $x={a\over b}$ and $y=Ma$. In the case of light
mass being examined now, $y<1$. Notice that Eq. (III.17) is valid, and exact,
for $x>{1\over 2\pi}$. In the case of a small or non existent extra dimension,
$x$ will be very large and the last term in the brackets will vanish, yielding
the well known result for the Casimir energy of a complex scalar field with
light mass. Only this last term in the brackets depends on $x$ and it contains
the entire dependence of $E$ on the compactification length $b$ of the extra
dimension. To understand how the presence of the large extra dimension affects
the Casimir energy, it is convenient to introduce the following dimensionless
function

$F_{D}(y,x)=1-{15\over\pi^{2}}y^{2}+{45\over\pi^{4}}y^{4}\log
y+{720\over\sqrt{2}\pi^{2}}x^{3/2}\sum_{m,n=1}^{\infty}{m^{3/2}\over
n^{5/2}}e^{-{ny^{2}\over 2\pi mx}}e^{-{4\pi mnx}},$ (III.18)

which allows me to write the Casimir energy per unit area as

${E\over L^{2}}=-{\pi^{2}\over 720a^{3}}F_{D}(y,x).$ (III.19)

Figure 1, below, displays $F_{D}(y,x)$ as a function of $x$ for three
different values of $y$ and shows that a deviation from the standard three-
dimensional Casimir energy starts to happen when $a\sim 0.65b$ for all three
values of $y$ and rapidly reaches an increase of more than $10\%$ over its
three-dimensional value when the plate distance $a\sim 0.40b$.

Figure 1:  Plot of $F_{D}(y,x)$ as a function of $x={a\over b}$, for $y=0$,
$y=0.1$, and $y=0.2$.

If a large compactified extra dimension is not present, then
$x\rightarrow\infty$ and $F_{D}(y,\infty)$ represents the correction to the
three dimensional Casimir energy due to a non-vanishing, but light, scalar
mass. For a massless scalar field, $y=0$ and $F_{D}(0,\infty)=1$. Figure 2
displays the ratio $F_{D}(y,x)\over F_{D}(y,\infty)$ as a function of $x$ for
the same three values of $y$ used in above.

Figure 2:  Plot of $F_{D}(y,x)\over F_{D}(y,\infty)$ as a function of
$x={a\over b}$, for $y=0$, $y=0.1$, and $y=0.2$.

Notice that the three curves are identical and therefore completely
superimposed onto each other, indicating that the effect of the extra
dimension on the Casimir energy appears for the same plate distance value
regardless of the value of the scalar field mass, as long as it is a light
mass. The value of $F_{D}(y,x)\over F_{D}(y,\infty)$ is $1.1$ for $x\sim 0.4$,
indicating that, when $a\sim 0.4b$, the Casimir energy in the presence of the
extra dimension increases by $10\%$ compared to its value in the absence of
the extra dimension. The increase rises to $20\%$ when $x\sim 0.3$.

Next, I investigate the case of large mass, $M>a^{-1}$. I start with the
expression of $\zeta_{1}(s)$ from Eq. (III.6), change variable of integration,
$t\rightarrow{na\over M}t$, and find

$\zeta_{1}(s)={L^{2}\over
8\pi^{2}\Gamma(s)}\sum_{n=1}^{\infty}\left({n\mu^{2}a\over
M}\right)^{s}{M^{2}\over
n^{2}a}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-nMa(t+t^{-1})}.$ (III.20)

Since $Ma>1$, I integrate using the saddle point method and obtain

$\zeta_{1}(s)={L^{2}\over 8\pi^{3/2}\Gamma(s)}\left({M\over
a}\right)^{3/2}\sum_{n=1}^{\infty}\left({n\mu^{2}a\over
M}\right)^{s}{e^{-2nMa}\over n^{5/2}}.$ (III.21)

I continue and evaluate the expression of $\zeta_{2}(s)$ from Eq. (III.7). I
change variable of integration,
$t\rightarrow{na\over\sqrt{M^{2}+{4\pi^{2}m^{2}\over b^{2}}}}t$, to obtain

$\zeta_{2}(s)={L^{2}\over
4\pi^{2}\Gamma(s)}\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\left({n\mu^{2}a\over\sqrt{M^{2}+{4\pi^{2}m^{2}\over
b^{2}}}}\right)^{s}\left({M^{2}+{4\pi^{2}m^{2}\over b^{2}}\over
n^{2}a}\right)\int_{0}^{\infty}dt\,\,t^{s-3}e^{-na\sqrt{M^{2}+{4\pi^{2}m^{2}\over
b^{2}}}(t+t^{-1})},$ (III.22)

then integrate using the saddle point method and find

$\zeta_{2}(s)={L^{2}\over
4\pi^{3/2}\Gamma(s)}\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\left({n\mu^{2}a\over\sqrt{M^{2}+{4\pi^{2}m^{2}\over
b^{2}}}}\right)^{s}\left({\sqrt{M^{2}+{4\pi^{2}m^{2}\over b^{2}}}\over
a}\right)^{3/2}{e^{-2na\sqrt{M^{2}+{4\pi^{2}m^{2}\over b^{2}}}}\over
n^{5/2}}.$ (III.23)

The Casimir energy $E$ is obtained from $\zeta_{1}(s)$ and $\zeta_{2}(s)$ by
using $E=-\zeta_{1}^{\prime}(0)-\zeta_{2}^{\prime}(0)$, and I find

${E\over L^{2}}=-{1\over 8\pi^{3/2}}\left({M\over
a}\right)^{3/2}\left[\sum_{n=1}^{\infty}{e^{-2nMa}\over
n^{5/2}}+2\sum_{m=1}^{\infty}\sum_{n=1}^{\infty}\left({1+{4\pi^{2}m^{2}\over
M^{2}b^{2}}}\right)^{3/4}{e^{-2na\sqrt{M^{2}+{4\pi^{2}m^{2}\over b^{2}}}}\over
n^{5/2}}\right].$ (III.24)

Since $Ma>1$, only terms with $n=1$ will contribute significantly and, using
the two dimensionless parameters $x=a/b$ and $y=Ma$ introduced previously, I
obtain

${E\over L^{2}}=-{1\over 8\pi^{3/2}}\left({M\over
a}\right)^{3/2}e^{-2Ma}\left[1+2\sum_{m=1}^{\infty}\left({1+{4\pi^{2}m^{2}x^{2}\over
y^{2}}}\right)^{3/4}e^{-2y\left(\sqrt{1+{4\pi^{2}m^{2}x^{2}\over
y^{2}}}-1\right)}\right].$ (III.25)

Notice that the first of the two terms in Eq. (III.25) is the well known
three-dimensional Casimir energy per unit area due to a scalar field with
large mass, while the second term is the correction due to the large
compactified extra dimension. To facilitate understanding the role of the
extra dimension, I introduce the following dimensionless function

$H(y,x)=1+2\sum_{m=1}^{\infty}\left({1+{4\pi^{2}m^{2}x^{2}\over
y^{2}}}\right)^{3/4}e^{-2y\left(\sqrt{1+{4\pi^{2}m^{2}x^{2}\over
y^{2}}}-1\right)},$ (III.26)

which allows me to write the Casimir energy per unit area by factoring out its
three dimensional form

${E\over L^{2}}=-{1\over 8\pi^{3/2}}\left({M\over
a}\right)^{3/2}e^{-2Ma}H(y,x).$ (III.27)

Figure 3, below, displays $H(y,x)$ as a function of $x$ for three different
values of $y$ showing a deviation from the standard three dimensional Casimir
effect when $y=1$ for $x\sim 0.7$, when $y=5$ for $x\sim 1$, and when $y=10$
for $x\sim 1.3$.

Figure 3:  Plot of $H(y,x)$ as a function of $x={a\over b}$, for $y=1$, $y=5$,
and $y=10$.

When comparing Figs. 1 and 3 it is clear that, in the large scalar mass
scenario, significant deviations from the standard three dimensional Casimir
effect happen for larger values of the plate distance $a$ than they do in the
low mass case and, in the large mass case, deviations grow large faster than
in the case of light mass.

##  IV Mixed boundary conditions

When investigating mixed boundary conditions, I use identity (III.1) in Eq.
(II.9), do the $k$-integration and obtain

$\zeta(s)={\mu^{2s}\over\Gamma(s)}{L^{2}\over(2\sqrt{\pi})^{3}}\sum_{n=0}^{\infty}\sum_{m=-\infty}^{\infty}\int_{0}^{\infty}dt\,\,t^{s-5/2}e^{-\left[(n+{1\over
2})^{2}({\pi\over a})^{2}+{4m^{2}\pi^{2}\over b^{2}}+M^{2}\right]t}.$ (IV.1)

Next, I do a Poisson resummation of the $n$-sum

$\sum_{n=0}^{\infty}e^{-(n+{1\over 2})^{2}({\pi\over a})^{2}t}={a\over
2\sqrt{\pi t}}\sum_{n=-\infty}^{\infty}(-1)^{n}e^{-{n^{2}a^{2}\over t}},$
(IV.2)

and obtain

$\zeta(s)={\mu^{2s}\over\Gamma(s)}{L^{2}a\over(4\pi)^{2}}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-M^{2}t}\left(1+2\sum_{n=1}^{\infty}(-1)^{n}e^{-{n^{2}a^{2}\over
t}}\right)\left(1+2\sum_{m=1}^{\infty}e^{-{4m^{2}\pi^{2}\over
b^{2}}t}\right).$ (IV.3)

Neglecting uniform energy density terms with only a linear dependence on $a$,
I write again $\zeta(s)=\zeta_{1}(s)+\zeta_{2}(s)$, where $\zeta_{1}(s)$
depends only on the plate distance $a$

$\zeta_{1}(s)={\mu^{2s}\over\Gamma(s)}{L^{2}a\over
8\pi^{2}}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-M^{2}t}\sum_{n=1}^{\infty}(-1)^{n}e^{-{n^{2}a^{2}\over
t}},$ (IV.4)

and $\zeta_{2}(s)$ depends also on the compactification length $b$ of the
large extra dimension

$\zeta_{2}(s)={\mu^{2s}\over\Gamma(s)}{L^{2}a\over
4\pi^{2}}\int_{0}^{\infty}dt\,\,t^{s-3}e^{-M^{2}t}\sum_{n=1}^{\infty}(-1)^{n}e^{-{n^{2}a^{2}\over
t}}\sum_{m=1}^{\infty}e^{-{4m^{2}\pi^{2}\over b^{2}}t}.$ (IV.5)

In the low mass scenario, I do a power series expansion of $e^{-M^{2}t}$ in
$\zeta_{1}(s)$ and retain the first three terms, change the integration
variable as in the previous section, and obtain

$\displaystyle\zeta_{1}(s)$ $\displaystyle=$ $\displaystyle{(\mu
a)^{2s}\over\Gamma(s)}{L^{2}\over
8\pi^{2}a^{3}}\biggl{[}(2^{2s-3}-1)\zeta_{R}(4-2s)\Gamma(2-s)-M^{2}a^{2}(2^{2s-1}-1)\zeta_{R}(2-2s)\Gamma(1-s)$
(IV.6) $\displaystyle+$ $\displaystyle\left.{M^{4}a^{4}\over
2}(2^{2s+1}-1)\zeta_{R}(-2s)\Gamma(-s)\right],$

which, using the following identities valid for $s\ll 1$

${(\mu
a)^{2s}\over\Gamma(s)}(2^{2s-3}-1)\zeta_{R}(4-2s)\Gamma(2-s)\simeq-{7\pi^{4}\over
720}s+{\cal O}(s^{2}),$ (IV.7) ${(\mu
a)^{2s}\over\Gamma(s)}(2^{2s-1}-1)\zeta_{R}(2-2s)\Gamma(1-s)\simeq-{\pi^{2}\over
12}s+{\cal O}(s^{2}),$ (IV.8) ${(\mu
a)^{2s}\over\Gamma(s)}(2^{2s+1}-1)\zeta_{R}(-2s)\Gamma(-s)\simeq\left({1\over
2s}+\gamma_{E}+\log{2\mu a\over\pi}\right)s+{\cal O}(s^{2}),$ (IV.9)

contributes an amount $E_{1}$, shown below, to the Casimir energy

$E_{1}={7\pi^{2}\over 5,760}{L^{2}\over a^{3}}\left[1-{60\over
7\pi^{2}}(Ma)^{2}-{360\over 7\pi^{4}}(Ma)^{4}\log(Ma)\right],$ (IV.10)

where a uniform energy density term has been neglected. The second part of the
zeta function, $\zeta_{2}(s)$, is obtained as it is done in the previous
section, with an extra factor of $(-1)^{n}$ being the only difference, and its
contribution to the Casimir energy is

$E_{2}=-{L^{2}\over\sqrt{2}(ab)^{3/2}}\sum_{m,n=1}^{\infty}(-1)^{n}{m^{3/2}\over
n^{5/2}}e^{-(M^{2}{nab\over 2\pi m})}e^{-{4\pi mna\over b}}.$ (IV.11)

The Casimir energy in the case of low mass, $E=E_{1}+E_{2}$, can be written as

${E\over L^{2}}={7\pi^{2}\over 5,760a^{3}}F_{m}(y,x),$ (IV.12)

where

$F_{m}(y,x)=1-{60\over 7\pi^{2}}y^{2}-{360\over 7\pi^{4}}y^{4}\log
y-{5,760\over
7\sqrt{2}\pi^{2}}x^{3/2}\sum_{m,n=1}^{\infty}(-1)^{n}{m^{3/2}\over
n^{5/2}}e^{-{ny^{2}\over 2\pi mx}}e^{-{4\pi mnx}}.$ (IV.13)

Figure 4 displays $F_{m}(y,x)$ as a function of $x$ for the same three values
of $y$ used previously and shows, one more time, that a deviation from the
standard mixed boundary conditions Casimir effect starts to happen when $x\sim
0.65$ for all three values of $y$. Notice that an increase of about $10\%$ of
the Casimir energy over its three-dimensional value is observed when $x\sim
0.4$, indicating again a rapid deviation of the Casimir effect with one
compactified extra dimension from the standard effect also for mixed boundary
conditions. Finally, it is important to point out that the repulsive nature of
the Casimir effect with mixed boundary conditions is maintained and, for small
$x$, the repulsion is enhanced. If this feature of the repulsive Casimir
effect is confirmed also for the spherical geometry, it could have interesting
cosmological implications. An investigation of the spherical Casimir effect in
the presence of a large compactified extra dimension is forthcoming.

As it happens in the case of Dirichlet boundary conditions, the three curves
of the ratio $F_{m}(y,x)\over F_{m}(y,\infty)$ as a function of $x$ are
identical to each other and to the curve $F_{m}(0,x)$ shown in Figure 4,
indicating that the Casimir energy in the presence of the compactified extra
dimension exceeds the three-dimensional Casimir energy by $10\%$ or more when
$a\sim 0.4b$ or less.

Figure 4:  Plot of $F_{m}(y,x)$ as a function of $x={a\over b}$, for $y=0$,
$y=0.1$, and $y=0.2$.

The large mass limit, $M>a^{-1}$, for mixed boundary conditions is obtained by
changing integration variable in $\zeta_{1}(s)$, $t\rightarrow{na\over M}t$,
and in $\zeta_{2}(s)$, $t\rightarrow{na\over\sqrt{M^{2}+{4\pi^{2}m^{2}\over
b^{2}}}}t$. Following the same steps taken in the case of Dirichlet boundary
conditions, I find

${E\over L^{2}}={1\over 8\pi^{3/2}}\left({M\over
a}\right)^{3/2}e^{-2Ma}H(y,x),$ (IV.14)

where the only difference with the Dirichlet boundary conditions case is the
lack of the overall negative sign, as expected for mixed boundary conditions.
The function $H(y,x)$ is displayed in Fig. 3 for three different values of
$y$. Notice that the overall positive sign for the case of mixed boundary
conditions produces a repulsive Casimir force, as it is well known. Notice
also that, for $x\rightarrow\infty$, the well known result for mixed boundary
conditions in the three-dimensional case is recovered.

##  V Casimir pressure

The Casimir pressure is obtained by taking a derivative of the Casimir energy

$P=-{1\over L^{2}}{\partial E\over\partial a}.$ (V.1)

For Dirichlet boundary conditions and light mass, I find

$P=-{\pi^{2}\over 240a^{4}}G_{D}(y,x),$ (V.2)

with

$G_{D}(y,x)=1-{5\over\pi^{2}}y^{2}-{15\over\pi^{4}}y^{4}(1+\log
y)+{240\over\sqrt{2}\pi^{2}}x^{3/2}\sum_{m,n=1}^{\infty}\left({m\over
n}\right)^{3/2}e^{-{ny^{2}\over 2\pi mx}}e^{-{4\pi mnx}}\left({3\over 2n}+4\pi
mx+{ny\over 2\pi mx}\right).$ (V.3)

Figure 5 shows $G_{D}(y,x)$ as a function of $x$ for the same three values of
$y$ used previously. As it is the case for the Casimir energy, we observe a
deviation from the standard three-dimensional Casimir pressure starting to
happen when $a\sim 0.65b$ for all three values of $y$ and reaching an increase
of more than $10\%$ over its three-dimensional value when the plate distance
$a\sim 0.40b$.

Figure 5:  Plot of $G_{D}(y,x)$ as a function of $x={a\over b}$, for $y=0$,
$y=0.1$, and $y=0.2$.

In the case of Dirichlet boundary conditions and heavy mass, $M>a^{-1}$, I
obtain the following

$P=-{M^{\frac{5}{2}}\over 4\pi^{\frac{3}{2}}}{e^{-2Ma}\over
a^{\frac{3}{2}}}J(y,x),$ (V.4)

with

$J(y,x)=1+2\sum_{m=1}^{\infty}\left({1+{4\pi^{2}m^{2}x^{2}\over
y^{2}}}\right)^{\frac{5}{4}}e^{-2y\left(\sqrt{1+{4\pi^{2}m^{2}x^{2}\over
y^{2}}}-1\right)}.$ (V.5)

Figure 6 shows $J(y,x)$ as a function of $x$ for the same three values of $y$
used in the previous section for the case of heavy mass. A deviation from the
standard three-dimensional Casimir pressure starts to happen for $x\sim 0.9$
when $y=1$, for $x\sim 1$ when $y=5$, and for $x\sim 1.3$ when $y=10$. Notice
that when $a\gg b$ and thus $x$ is large, the well known Casimir pressure for
the three-dimensional case of heavy mass is recovered.

Figure 6:  Plot of $J(y,x)$ as a function of $x={a\over b}$, for $y=1$, $y=5$,
and $y=10$.

When investigating mixed boundary conditions and light mass, I obtain

$P={7\pi^{2}\over 1,920a^{4}}G_{m}(y,x),$ (V.6)

where

$G_{m}(y,x)=1-{20\over 7\pi^{2}}y^{2}+{120\over 7\pi^{4}}y^{4}(1+\log
y)-{1,920\over
7\sqrt{2}\pi^{2}}x^{3/2}\sum_{m,n=1}^{\infty}(-1)^{n}\left({m\over
n}\right)^{3/2}e^{-{ny^{2}\over 2\pi mx}}e^{-{4\pi mnx}}\left({3\over 2n}+4\pi
mx+{ny\over 2\pi mx}\right).$ (V.7)

I display $G_{m}(y,x)$ as a function of $x$ in Figure 7, for the same three
values of $y$ used before. Again, we observe a deviation from the standard
three-dimensional Casimir pressure starting to happen when $a\sim 0.65b$ for
all three values of $y$. It happens also in this case that the well known
three dimensional Casimir pressure is obtained when $a\gg b$.

Figure 7:  Plot of $G_{m}(y,x)$ as a function of $x={a\over b}$, for $y=0$,
$y=0.1$, and $y=0.2$.

For mixed boundary conditions and heavy mass, $M>a^{-1}$, I obtain the
following

$P={M^{\frac{5}{2}}\over 4\pi^{\frac{3}{2}}}{e^{-2Ma}\over
a^{\frac{3}{2}}}J(y,x),$ (V.8)

where the pressure has the same magnitude as in the case of Dirichlet boundary
conditions and heavy mass, but opposite sign. $J(y,x)$ is the same function
shown in Eq. (V.5) for Dirichlet boundary conditions and is displayed in
Figure 6.

##  VI Discussion and conclusions

In this work I used the zeta function technique to study the Casimir effect
due to a massive complex scalar field that also permeates a compactified
fourth dimension. I investigated in detail Dirichlet and mixed boundary
conditions. Neumann boundary conditions have not been considered since they
produce the same results as Dirichlet.

I obtained simple and exact analytic expressions for the Casimir energy in the
case of Dirichlet boundary conditions and low mass, Eq. (III.19) with the
auxiliary Eq. (III.18) for $F_{D}(y,x)$, and heavy mass, Eq. (III.27) with the
auxiliary Eq. (III.26) for $H_{D}(y,x)$. Notice that, while the dimensionless
functions $F_{D}$ and $H_{D}$ contain each an infinite series, these series
are quickly convergent and the first two or three terms alone approximate the
full sum to within one percent or less. Figs. 1 and 2 display $F_{D}(y,x)$ as
a function of $x$ for three values of $y$ and show that in the case of low
mass, a deviation from the standard Casimir energy starts to happen when
$a\sim 0.65b$, independently of the value of $M$, and quickly increases. Fig.
3 displays $H_{D}(y,x)$ and shows that in the heavy mass scenario, deviations
from the standard Casimir energy start to happen when $a\sim 0.9b$ for $Ma=5$
and when $a\sim 1.2b$ for $Ma=10$ and increase rapidly. It is clear from the
results of this work that the Casimir effect due to a field with heavy mass is
affected more strongly by the presence of an accessible compactified fourth
dimension. My results for the Casimir pressure, the one quantity measurable in
the lab, are also exact simple analytic expressions. In the case of Dirichlet
boundary conditions these results are shown in Eqs. (V.2), (V.3), (V.4),
(V.5), and Figs. 5, 6, confirming fully the onset of a discrepancy with the
standard Casimir pressure at the values of $a$ that I report above, and
confirming that the Casimir pressure due to a heavily massive scalar field is
affected more strongly by the compactified fourth dimension. The infinite
series contained in the two dimensionless functions, $G_{D}(y,x)$ and
$J(y,x)$, that appear in my Casimir pressure results are also quickly
convergent and can accurately be approximated by their first two or three
terms.

I draw the same conclusions from my investigation of the Casimir energy and
pressure for the case of mixed boundary conditions. My results show that the
Casimir force is repulsive for a scalar field, massive or not, that permeates
a compactified extra dimension. Therefore the presence of a compactified extra
dimension does not alter the repulsive nature of the mixed boundary conditions
effect. The onset of Casimir energy and pressure deviations from their three-
dimensional values happens at the same distances $a$ reported in the previous
paragraph, both in the case of light mass and heavy mass, as can be seen from
Figs. 4 and 7 for $F_{m}$ and $G_{m}$ respectively.

## References

  * (1) H. B. G. Casimir, Indag. Math. 10, 261 (1948) [Kon. Ned. Akad. Wetensch. Proc. 51, 793 (1948)] [Front. Phys. 65, 342 (1987)] [Kon. Ned. Akad. Wetensch. Proc. 100N3-4, 61 (1997)].
  * (2) M. J. Sparnaay, Physica 24, 751 (1958).
  * (3) M. Bordag, U. Mohideen, and V. M. Mostepanenko, Phys. Rep. 353, 1 (2001).
  * (4) M. Bordag, G. L. Klimchitskaya, U. Mohideen, and V. M. Mostepanenko, International series of monographs on physics. 145.
  * (5) T. H. Boyer, Phys. Rev. A 9, 2078 (1974).
  * (6) T. H. Boyer, Phys. Rev. 174, 1764 (1968).
  * (7) T. Kaluza, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys. ) 1921, 966-972 (1921)
  * (8) O. Klein, Z. Phys. 37, 895-906 (1926)
  * (9) N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Lett. B 429, 263-272 (1998)
  * (10) I. Antoniadis, N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Lett. B 436, 257-263 (1998)
  * (11) N. Arkani-Hamed, S. Dimopoulos and G. R. Dvali, Phys. Rev. D 59, 086004 (1999)
  * (12) L. Randall and R. Sundrum, Phys. Rev. Lett. 83, 3370-3373 (1999)
  * (13) E. G. Adelberger, J. H. Gundlach, B. R. Heckel, S. Hoedl and S. Schlamminger, Prog. Part. Nucl. Phys. 62, 102-134 (2009)
  * (14) J. Murata and S. Tanaka, Class. Quant. Grav. 32, no.3, 033001 (2015)
  * (15) L. P. Teo, Nucl. Phys. B 819, 431-452 (2009)
  * (16) L. P. Teo, Phys. Lett. B 672, 190-195 (2009)
  * (17) S. W. Hawking, Commun. Math. Phys. 55, 133 (1977).
  * (18) E. Elizalde and A. Romeo, J. Math. Phys. 30, 1133 (1989) [erratum: J. Math. Phys. 31, 771 (1990)]
  * (19) E. Elizalde, J. Phys. A 41, 304040 (2008).