The success of deep learning over the past decade mainly relies on gradient-based optimisation and backpropagation. This paper focuses on analysing the performance of first-order gradient-based optimisation algorithms, gradient descent and proximal gradient, with time-varying non-convex cost function under (proximal) Polyak-Łojasiewicz condition. Specifically, we focus on using the forward mode of automatic differentiation to compute gradients in the fast-changing problems where calculating gradients using the backpropagation algorithm is either impossible or inefficient. Upper bounds for tracking and asymptotic errors are derived for various cases, showing the linear convergence to a solution or a neighbourhood of an optimal solution, where the convergence rate decreases with the increase in the dimension of the problem. We show that for a solver with constraints on computing resources, the number of forward gradient iterations at each step can be a design parameter that trades off between the tracking performance and computing constraints.

Online optimisation, forward gradient, over-parameterised systems, PL condition

§ INTRODUCTION

Gradient-based optimisation is the core of most machine learning algorithms <cit.>. Backpropagation, the reverse mode of Automatic Differentiation (AD) algorithms, has been the main algorithm for computing gradients. <cit.> has argued that in training neural networks, one can calculate gradients based on directional derivatives, which is faster than backpropagation. An unbiased estimate of gradients, named forward gradients, can be calculated in a single forward run of the function and speeds up training up to twice as fast in some examples. Forward differentiation is a mode of automatic differentiation algorithm. It is based on the mathematical definition of dual numbers (see <cit.> for a review).

Given a function $f: \mathbb{R}^n \to \mathbb{R}$, the forward gradient at point $\boldsymbol{x}$ is defined as

\begin{align}
\boldsymbol{g}(\boldsymbol{x}, \boldsymbol{u}) = \langle \nabla f(\boldsymbol{x}), \boldsymbol{u}\rangle \boldsymbol{u},
\end{align}

where $\boldsymbol{u}$ is a random vector with zero mean and unit variance. Note that $\langle \nabla f(\boldsymbol{x}), \boldsymbol{u} \rangle$ is the directional derivative of $f$ at point $\boldsymbol{x}$ in the direction $\boldsymbol{u}$. The forward gradient, $g(\boldsymbol{x}, \boldsymbol{u})$, is an unbiased estimator for $\nabla f(\boldsymbol{x})$ as long as the components of direction $\boldsymbol{u}$ are independent and sampled from a distribution with zero mean and unit variance; $\left[\mathbb{E}(g(\boldsymbol{x}, \boldsymbol{u}))\right]_i = \mathbb{E}(u_i \sum_{j=1}^n (\nabla f(\boldsymbol{x}))_j u_j) = (\nabla f(\boldsymbol{x}))_i$, where $a_i$ is the $i$-th element of vector $\boldsymbol{a}$.

Deep neural networks can be considered as a particular case of the over-parameterised system of nonlinear equations in which there are potentially a staggering number of trainable parameters <cit.>. In <cit.> authors argue that focusing on the convexity of the problem does not lead to a suitable framework for analysing such loss functions and instead show that the Ployak-Łojasiewicz (PL) condition <cit.> is satisfied on the most of parameter space. This in turn explains fast convergence of Gradient Descent (GD) algorithm to a global minimum.

Another area that has observed a flurry of activity is anchored at studying time-varying objective functions. These problems arise in many applications including but not limited to robotics <cit.>, smart grids <cit.>, communication systems <cit.> and signal processing <cit.>. An example is when a machine learning problem's training data (problem parameters) are observed during the learning process. These can be considered as a sequence of optimisation problems and a solver may solve each problem completely before approaching the next problem, which is not feasible in cases where each instance of the problem is of large scale (considering the rate of change and the speed of solver) or involves communication over a network.  Therefore, a more efficient way may be to solve each problem partially (e.g. performing only a few gradient descent steps) before receiving the next problem. These algorithms are called running method in <cit.>. Consider a neural network that has been trained over a data set $\mathcal{D} = \{\boldsymbol{a}_i, b_i \}_{i=0}^{n-1}$. Let $\mathcal{L}_0(\boldsymbol{x}) = \frac{1}{2}\sum_{i=0}^{n-1} \|f(\boldsymbol{x} ; \boldsymbol{a}_i) - b_i \|^2$. Assume that, while deploying the pre-trained neural network, new samples of data arrive at a high rate and one needs to update the training variable of the network, $\boldsymbol{x}$. This is an over-parameterised online optimisation problem. We assume that only the latest $n$ data samples to update the variables are used. The corresponding cost functions is

\begin{align}
\mathcal{L}_k(\boldsymbol{x}) = \frac{1}{2}\sum_{i=n_k}^{n_k+n-1} \| f(\boldsymbol{x} ; \boldsymbol{a}_i) - b_i \|^2,
\end{align}

where $n_k - n_{k-1}$ is the number of newly observed samples.
Note that we assume that the structure of the neural network does not change during time and we only try to update its variables after observing new data examples (See <cit.> and references therein for a review on expandable neural networks). We assume

* The algorithm uses the latest samples (because they represent the data distribution better);

* Only $n$ samples (e.g. because of memory limitations) can be used;

* Due to its fast computation, forward gradient is to be used in the optimisation steps.

Ideally, the solution generated by the method at each time step, $\boldsymbol{x}_k$, follows $\boldsymbol{x}_k^*$. In other words, the tracking error $\|\boldsymbol{x}_k - \boldsymbol{x}_k^* \|$ or the instantaneous sub-optimality gap $\|\mathcal{L}_k(\boldsymbol{x}_k) - \mathcal{L}_k(\boldsymbol{x}_k^*) \|$ either shrink or remains bounded as $k$ grows. We assume that $\ell$ optimisation step is applied to $\mathcal{L}_k$ before receiving the new batch of data, where $\ell$ is chosen based on the above constraints.

In this paper, we focus on the analysis of the performance of forward gradient-based method in over-parameterised time-varying settings characterised by the following cost function at time $k$:

\begin{align}\label{composite_sum_timevar}
\min_{\boldsymbol{x}\in \mathbb{R}^m} \ \{ \mathcal{L}_k(\boldsymbol{x}) := g_k(\boldsymbol{x}) + h_k(\boldsymbol{x}) \}.
\end{align}

Function $g_k$ is smooth and potentially nonconvex and $h_k(\boldsymbol{x})$ is possibly nonsmooth convex function that imposes some structure on the solution such as sparsity. We consider two cases with or without the presence of function $h$ in the cost function. Such problems arise in the online machine learning <cit.> and have been well studied under (strong) convexity assumption. The contributions of this paper are listed below:

* We prove the linear convergence in expectation of forward-gradient method to the minimiser for smooth nonconvex PL objective functions (Theorem <ref>).

* We prove that the expected tracking error of the forward-gradient method in a time-varying nonconvex setting, under smoothness and PL assumptions, converges to a neighbourhood of zero at a linear rate (Theorem <ref>). We show that the number of iterations at each step, $\ell$, is a design parameter that trades off between the tracking performance and computing resources constraints.

* In the case where the objective cost function is the sum of a PL and a convex (possibly non-smooth) function, we prove that the expected tracking error of the proximal forward-gradient-based algorithms converges to a neighbourhood of zero for time-varying cost functions under proximal-PL condition (Theorem <ref>). Error depend on the problem dimension, but the rate does not. The proof for quadratic growth of the proximal-PL functions is also novel.

The proofs are described in detail in the technical report associated with this paper <cit.>.

The rest of the paper is organised as follows: in Section <ref>, we discuss related works, and in Section <ref>, we review notations, definitions, basic assumptions and some background knowledge. In Section <ref>, we prove the linear convergence of forward gradient algorithm for the PL objective functions. In Section <ref>, we analyse the performance of gradient descent and proximal gradient algorithms under PL and Proximal-PL conditions. Section <ref> concludes the paper.

§ RELATED WORKS

PL condition in over-parameterised nonlinear systems: The nonlinear least square problem has been extensively studied in under-parameterised systems <cit.>. In <cit.> the sufficient condition for satisfying the PL condition in a wide neural network with least square loss function and linear convergence of (S)GD have been studied and the relation with neural tangent kernel has been discussed <cit.>.

Static cost functions: In <cit.> convergence of Randomised Coordinate Descent (RCD), Stochastic Gradient Descent, and SVRG algorithms under PL condition and Proximal coordinate descent algorithm under Proximal-PL has been studied. Inexact gradient descent algorithms have been analysed under a wide range assumptions (See <cit.> and references therein). In <cit.> $\mathcal{O}(1/k)$ convergence of SGD under weaker condition (Weak Growth Condition) has been proved.

Time-varying cost function: A closely related line of research to this paper is online convex optimisation (OCO) which was originally introduced in seminal work by <cit.>. Since the focus of this work is on nonconvex optimisation first order algorithm with non-exact gradients, we omit to discuss the vast literature on the online exact first order algorithms under strong convexity and smoothness assumptions (see <cit.> and references therein). <cit.> analysed the regret bounds for online first order optimisation algorithms under (strong) convexity and smoothness with noisy gradients. Regret analysis of online stochastic fixed and time-varying optimisation have been done in <cit.> and <cit.>. <cit.> and <cit.> have analysed the online (proximal) gradient methods with sub-Weibull gradient error under strong convexity and PL conditions, respectively. In <cit.> authors have proposed a method for estimating the gradient in a time-varying setting with zeroth-order distributed optimisation algorithm.

§ NOTATION AND PRELIMINARIES

Throughout this report, we show the set of real numbers as $\mathbb{R}$. For vectors $\boldsymbol{a}, \boldsymbol{b}\in \mathbb{R}^n$ the Euclidean inner product and its corresponding norms are denoted by $\langle \boldsymbol{a} , \boldsymbol{b} \rangle$ and $\| \boldsymbol{a}\|$ respectively. We denote the open ball centered at $\boldsymbol{x}$ with radius $r$ by $B(\boldsymbol{x}, r)$. With $\pi_{\mathcal{X}}(\boldsymbol{x})$ we indicate the projection of vector $\boldsymbol{x}$ into set $\mathcal{X}$. The Euclidean distance to  a set $\mathcal{X} \in \mathbb{R}^n$ is defined by $dist(\boldsymbol{x}, \mathcal{X}):= \inf_{\boldsymbol{u} \in \mathcal{X}} \| \boldsymbol{u} - \boldsymbol{x} \|$ for every $\boldsymbol{x}\in \mathbb{R}^n$. We denote the mapping of projection over set $\mathcal{X}$ with $\pi_{\mathcal{X}}$. With $\mathcal{X}^*$ we denote the set of minimisers of a function. For a matrix $A\in \mathbb{R}^{m\times n}$, $\| A\|$ and $\| A\|_F$ are the spectral and Frobenius norms of the matrix, respectively. We use $\mathcal{D} F(\cdot)$ and $\mathcal{D}^2 F(\cdot)$ to show the derivative and Hessian of a function $F:\mathbb{R}^n \to \mathbb{R}$, respectively. We use $\widetilde{\mathcal{O}}(\cdot)$ notation by ignoring the logarithmic factors in big-O notation, i.e. $\widetilde{\mathcal{O}}(f(n)) = \mathcal{O}(f(n)\log f(n))$.

In what follows, we introduce some definitions, assumptions and lemmas that will be utilised in the paper.

The function $\mathcal{L}: \mathbb{R}^n \to \mathbb{R}$ is $\beta-smooth$ if it is differentiable and

\begin{align*}
\|\nabla \mathcal{L}(\boldsymbol{y}) - \nabla \mathcal{L}(\boldsymbol{x}) \| \leq \beta\| \boldsymbol{y} - \boldsymbol{x} \|, \quad \forall \boldsymbol{x}, \boldsymbol{y} \in \mathbb{R}^n.
\end{align*}

The function $\mathcal{L}: \mathbb{R}^n \to \mathbb{R}$ satisfies the Polyak-Łojasiewicz (PL) condition on a set $\mathcal{X}$ with constant $\mu$ if

\begin{align}
\frac{1}{2} \|\nabla \mathcal{L}(\boldsymbol{x}) \|^2 \geq \mu (\mathcal{L}(\boldsymbol{x}) - \mathcal{L}^*), \quad \forall \boldsymbol{x} \in \mathcal{X}.
\end{align}

(Descent Lemma <cit.>) Let the function $\mathcal{L}: \mathbb{R}^n \to \mathbb{R}$ be a $\beta_{\mathcal{L}}$-smooth function. Then for every $\boldsymbol{x}$, $\boldsymbol{y} \in \mathbb{R}^n$ and for every $\boldsymbol{z} \in [\boldsymbol{x}, \boldsymbol{y}] := \{(1-\alpha)\boldsymbol{x} + \alpha \boldsymbol{y} : \alpha \in [0,1] \}$ the following holds

\begin{align*}
\mathcal{L}(\boldsymbol{y}) \leq \mathcal{L}(\boldsymbol{x}) + \langle \nabla \mathcal{L}(\boldsymbol{z}), \boldsymbol{y} - \boldsymbol{x}\rangle + \frac{\beta_{\mathcal{L}}}{2} \|\boldsymbol{y} - \boldsymbol{x} \|^2.
\end{align*}

(<cit.>) Let $\mathcal{L}: \mathbb{R}^n \to \mathbb{R}$ be a $\beta$-smooth and $\mu-PL$ function. Then

\begin{align*}
\frac{\mu}{2}\|\boldsymbol{x} - \pi_{\mathcal{X}^*}(\boldsymbol{x}) \|^2\leq \mathcal{L}(\boldsymbol{x}) - \mathcal{L}^* \leq \frac{\beta}{2}\|\boldsymbol{x} - \pi_{\mathcal{X}^*}(\boldsymbol{x})  \|^2, \quad \forall \boldsymbol{x}.
\end{align*}

(Proximal-PL condition <cit.>) Let $\mathcal{L}(\boldsymbol{x}) = g(\boldsymbol{x})+ h(\boldsymbol{x})$ where $g$ is $\beta-$smooth and $h$ is a convex function. The function $f$ satisfies $\mu-$Proximal PL-condition if the following holds

\begin{align}\label{prox PL ineq}
\frac{1}{2}\mathcal{D}_h(\boldsymbol{x}, \beta) \geq \mu(\mathcal{L}(\boldsymbol{x}) - \mathcal{L}^*), \quad \forall \boldsymbol{x},
\end{align}

\begin{align}
\mathcal{D}_h(\boldsymbol{x}, \alpha):= -2\alpha \min_{\boldsymbol{y}} \left\{ \langle \nabla g(\boldsymbol{x}), \boldsymbol{y}-\boldsymbol{x} \rangle + \frac{\alpha}{2}\|\boldsymbol{y}-\boldsymbol{x} \|^2 + h(\boldsymbol{y}) - h(\boldsymbol{x})\right\}
\end{align}

§.§ Over-parameterised systems and PL condition

For a system of $n$ nonlinear equations

\begin{align*}
f(\boldsymbol{x}; \boldsymbol{a}_i) = b_i, \quad i=1,2, \dots, n,
\end{align*}

where $\{\boldsymbol{a}_i, b_i \}_{i=1}^n$ is the set of model parameters and one aims at finding $\boldsymbol{x}\in \mathbb{R}^m$ that solves the system of equations. Aggregating all equations in a single map amounts to

\begin{align}\label{Single-map}
\mathcal{F}(\boldsymbol{x}) = \boldsymbol{y}, \ \text{where} \ \boldsymbol{x}\in \mathbb{R}^m, \mathcal{F}(\cdot):\mathbb{R}^m \to \mathbb{R}^n.
\end{align}

The system in (<ref>) is solved through minimising a certain loss function $\mathcal{L}(\boldsymbol{x})$:

\begin{align}
\mathcal{L}(\boldsymbol{x}):=\frac{1}{2} \|\mathcal{F}(\boldsymbol{x}) - \boldsymbol{b} \|^2=\frac{1}{2} \sum_{i=1}^n (f(\boldsymbol{x}; \boldsymbol{a}_i) - b_i)^2.
\end{align}

This problem has been studied extensively in under-parameterised setting (where $m<n$). We refer to exact solution of (<ref>) as interpolation.

Let $D\mathcal{F}(\boldsymbol{x})$ be the differential of the map $\mathcal{F}$ at $\boldsymbol{x}$ which can be represented as a $n\times m$ matrix. The tangent Kernel of $\mathcal{F}$ is defined as a $n\times n$ positive semidefinite matrix $K(\boldsymbol{x}):=D\mathcal{F}(\boldsymbol{x})D\mathcal{F}^T(\boldsymbol{x})$.

We say that a non-negative function $\mathcal{L}$ satisfies $\mu-PL^*$ condition on a set $\mathcal{X}\in \mathbb{R}^m$ for $\mu >0$, if

\begin{align}\label{PL}
\frac{1}{2} \|\nabla \mathcal{L}(\boldsymbol{x}) \|^2 \geq \mu \mathcal{L}(\boldsymbol{x}), \quad \forall \boldsymbol{x} \in \mathcal{X}.
\end{align}

* if a non-negative function satisfies $\mu-PL^*$, then it will satisfy $\mu-PL$ condition too, that is

\begin{align}
\frac{1}{2}\|\nabla \mathcal{L} (\boldsymbol{x} \|^2 \geq \mu(\mathcal{L}(\boldsymbol{x}) - \mathcal{L}^*)
\end{align}

* Every stationary point of a PL function, is a global minimum.

The following results are of great importance in the later discussions.
If $\mathcal{L}(\boldsymbol{x})$ is lower bounded with $\beta$-Lipschitz continuous gradients and satisfies $\mu-PL$ condition in the region $\mathcal{X} = \overline{B(\boldsymbol{x_0}, \rho)}$ where $\rho > \frac{\sqrt{2\beta(\mathcal{L}(\boldsymbol{x}_0) - \mathcal{L}^*)}}{\mu}$, then there exists a global minimum point $x^* \in \mathcal{X}$ and the gradient descent algorithm with a small enough step-size ($\frac{1}{\beta}$) converges to $x^*$ in a linear rate.

If $\mathcal{F}(\boldsymbol{x})$ is such that $\lambda_{min}(K(\boldsymbol{x}))\geq \mu > 0$ for all $\boldsymbol{x} \in \mathcal{X}$, then the square loss function $\mathcal{L}(\boldsymbol{x}) = \frac{1}{2}\|\mathcal{F}(\boldsymbol{x})- \boldsymbol{y} \|^2$ satisfies $\mu-PL^*$ condition on $\mathcal{X}$.

Note that $K(\boldsymbol{x})=D\mathcal{F}(\boldsymbol{x})D\mathcal{F}^T(\boldsymbol{x})$ implies that $rank(K(\boldsymbol{x})) = rank(D\mathcal{F}(\boldsymbol{x}))$. In an over-parameterised system, $m>n$, starting from a random initial point there is a high probability that $D\mathcal{F}(\boldsymbol{x})$ is full rank. However, from Theorem <ref>, one needs $\mu-PL^*$ condition in a large enough ball with radius $\mathcal{O}(\frac{1}{\mu})$ for linear convergence of the GD algorithm. For establishing such conditions, one intuitively can expect that in the cases of $\mathcal{C}^2$ function that if Hessian (curvature) is small enough in the neighbourhood of a point, then the Tangent Kernel should be almost constant and therefore, the conditions for linear convergence of the GD algorithm will be satisfied. It turns out that in the case of highly over-parameterised Neural Networks (wide NNs) with linear output layer, the Hessian matrix will have arbitrary small spectral norm (a transition to linearity). This is formulated as follows:
For an $L$ layer neural network with a linear output layer and minimum width $m$ of the hidden layers, for any $R>0$ and any $\boldsymbol{x}\in B_{R}(\boldsymbol{x}_0)$. the Hessian spectral norm satisfies the following with a high probability

\begin{align}
\|\mathcal{D}^2 \mathcal{F}(\boldsymbol{x}) \| = \widetilde{\mathcal{O}}\left (\frac{R^{3L}}{\sqrt{m}} \right).
\end{align}

With a nonlinear output layer, $\boldsymbol{x} \mapsto \varphi(\boldsymbol{x})$, the square loss function will satisfy $\mu\kappa^2$-PL condition where $\kappa:=\inf_{\boldsymbol{x} \in B(\boldsymbol{x}_0, \rho)} \|\mathcal{D}\varphi(\mathcal{F}(\boldsymbol{x})) \|$.

All in all, wide neural networks, as a particular case of over-parameterised systems, satisfy the $\mu-PL^*$ condition, which explains the fast convergence of (S)GD to a global minimum in square-loss problems. In the following sections, we theoretically analyse the performance of using forward gradient in the same setting (over-parameterised systems).

§ OPTIMISATION USING FORWARD GRADIENT

In this section, we analyse the performance of various gradient-based algorithms using forward gradient.

§.§ Forward gradient under PL condition for fixed cost functions

With the focus on the basic unconstrained optimisation problem

\begin{align}
\min_{\boldsymbol{x}\in \mathbb{R}^m} \mathcal{L}(x).
\end{align}

Consider the following iterations where the solver is not fast enough to do one gradient calculation, but it is able to do $\ell$ forward gradient updates at time step $k$

\begin{align}\label{GD with directional derivative}
\boldsymbol{x}_{k+1}^{(i+1)} = \boldsymbol{x}_{k+1}^{(i)} - \alpha_k^{(i)} \boldsymbol{v}_{k}(\boldsymbol{x}_{k+1}^{(i)}, U_{k+1}^{(i)}), \quad \textit{for}\  i=0,1, \dots, \ell -1,
\end{align}

where  $\boldsymbol{x}_{k+1}^{(\ell)} := \boldsymbol{x}_{k+1}$, $\boldsymbol{x}_{k+1}^{(0)} := \boldsymbol{x}_k$, $\alpha_k$ is step-size, $\boldsymbol{v}(\boldsymbol{x}, \boldsymbol{u}) = \langle \nabla \mathcal{L}(\boldsymbol{x}), \boldsymbol{u}\rangle \boldsymbol{u}$, and $U_k \sim \mathcal{N}(\boldsymbol{0}, I_m)$ are i.i.d random directions for all $k$ and $i$. It is extensively studied in <cit.> where the following properties have been proved:

\begin{align}
\mathbb{E} \boldsymbol{v}(\boldsymbol{x}, U_k) &= \nabla \mathcal{L}(\boldsymbol{x})\\
\mathbb{E}(\|\boldsymbol{v}(\boldsymbol{x}, U_k) \|^2) &\leq (m+4)\| \nabla \mathcal{L}(\boldsymbol{x})\|^2 \label{Var-Dir_Deriv}
\end{align}

In the following we analyse the convergence of this algorithm under PL condition. We refer the reader to the extended version of this paper <cit.> for proofs.

Assume that function $\mathcal{L}$ is $\beta-$smooth, has a non-empty solution set $\mathcal{X}^*$, and satisfies $\mu-$PL condition. Consider the algorithm (<ref>) with a step-size of $\frac{1}{\beta(m+4)}$. If random vector $U_k$ is chosen from a standard normal distribution, i.e. and $U_k \sim \mathcal{N}(\boldsymbol{0}, I_m)$, then the algorithm has an expected linear convergence rate

\begin{align}
\mathbb{E}\{\mathcal{L}(x_k) - \mathcal{L}^* \} \leq \left(1-\frac{\mu}{(m+4)\beta}\right)^k (\mathcal{L}(x_0) - \mathcal{L}^*)
\end{align}

§.§ Time varying Optimisation

In what follows, we prove a convergence property of online line-search methods for the objective functions that satisfy the following assumptions.
(Polyak-Łojasiewicz Condition) There exist a scalar $\mu>0$ such that $\mathcal{L}_{k}(\boldsymbol{x})$ satisfies (<ref>) for all $k$  and for all $\boldsymbol{x} \in \mathbb{R}^m$.

We also need to quantify the speed with which the objective function varies in each step.
(Drift in Time): There exist non-negative scalars $\eta_0$ and $\eta^*$ such that $\mathcal{L}_{k+1}(\boldsymbol{x}) - \mathcal{L}_k(\boldsymbol{x}) \leq \eta_0$ for all $\boldsymbol{x} \in \mathbb{R}^m$ and $\mathcal{L}_{k+1}^* - \mathcal{L}_k^* \leq \eta^*$ for all $k$.

In an over-parameterised setting with non-linear least square loss function, one can argue that $\mathcal{L}_k^* = 0$, for all $k$. Therefore, $\eta^* = 0$ in this setting.

§.§.§ Forward gradient convergence in Time varying setting

In this part, we analyse the performance of the forward-gradient algorithm with a fixed step size for the time-varying nonconvex cost functions under PL and smoothness assumptions. We prove that the tracking error linearly converges to a neighbourhood of zero where the rate of convergence depends on the dimension of the problem and the asymptotic error is independent of the fact that we are not using exact gradients.

Assume that $\mathcal{L}_k$ is $\beta-$smooth for each $k$, and Assumptions <ref> and <ref> hold. If we use directional derivative algorithm with constant step-size $\alpha_k^{(i)} = \alpha < \frac{2}{\beta(m+4)}$ and random direction $U_k^{(i)}$ is chosen from a standard normal distribution, then

\begin{align}\label{err_bdd_timevar}
\mathbb{E}\{\| \boldsymbol{x}_{k+1} - \pi_{\mathcal{X}_{k+1}^*}(\boldsymbol{x}_{k+1}) \|^2\} &\leq \frac{\eta_0 + \eta^*}{\mu^2 \gamma \ell} + \frac{2}{\mu}(1-2\mu\gamma\ell)^k (\mathcal{L}_0(\boldsymbol{x}_0) - \mathcal{L}_0^*),
\end{align}

where $\eta_0$ and $\eta^*$ are as in Assumption <ref>, $\gamma \in (0, \tilde{\gamma})$ where $\tilde{\gamma} = min \{\alpha(1-\frac{\beta}{2}(m+4)\alpha), \frac{1}{2\mu\ell} \}$, and $\mathcal{X}_k^*$ is the set of minimisers of $\mathcal{L}_k(\boldsymbol{x})$.

Using descent lemma and (<ref>) we have

\begin{align*}
\mathcal{L}_{k}(\boldsymbol{x}_{k+1}^{(i+1)}) &\leq \mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i)}) + \langle \nabla \mathcal{L}_{k}(\boldsymbol{x}_{k+1}^{(i)}), \boldsymbol{x}_{k+1}^{(i+1)} - \boldsymbol{x}_{k+1}^{(i)} \rangle + \frac{\beta}{2}\|\boldsymbol{x}_{k+1}^{(i+1)} - \boldsymbol{x}_{k+1}^{(i)} \|^2 \\
& =  \mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i)}) -\alpha_k^{(i)} \langle \nabla \mathcal{L}_{k}(\boldsymbol{x}_{k+1}^{(i)}), \boldsymbol{v}_k(\boldsymbol{x}_{k+1}^{(i)}, U_{k+1}^{(i)}) \rangle + \frac{\beta}{2}(\alpha_k^{(i)})^2\|\boldsymbol{v}_k(\boldsymbol{x}_{k+1}^{(i)}, U_{k+1}^{(i)}) \|^2.
\end{align*}

By taking conditional expectation given $\boldsymbol{x}_{k+1}^{(i)}$ with respect to $U_{k+1}^{(i)}$, and using (<ref>) we obtain

\begin{align}
\mathbb{E}\{\mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i+1)}) \mid \boldsymbol{x}_{k+1}^{(i)} \} &\leq \mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i)}) -\alpha_{k}^{(i)}  \|\nabla \mathcal{L}_{k}(\boldsymbol{x}_{k+1}^{(i)}) \|^2 +\frac{1}{2}\beta (\alpha_k^{(i)})^2(m+4)\|\nabla \mathcal{L}_{k}(\boldsymbol{x}_{k+1}^{(i)}) \|^2 \nonumber\\
& = \mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i)}) -\underbrace{\alpha_{k}^{(i)}(1-\frac{\beta}{2}(m+4)\alpha_{k}^{(i)})}_{:=\gamma_{k}^{(i)}}\|\nabla\mathcal{L}_{k}(\boldsymbol{x}_{k+1}^{(i)}) \|^2 \label{descent_lemma_implication}.
\end{align}

We now fix $\epsilon_1  <\alpha_{k}^{(i)}= \alpha < \frac{2}{\beta(m+4)} - \epsilon_2$ and $\gamma_{k}^{(i)} = \gamma >0$ for some $\epsilon_1 >0$ and $\epsilon_2 >0$. Assuming that $\mathcal{L}_k$ satisfies $\mu-PL$ condition at $\boldsymbol{x}_{k+1}^{(i)}$, for all $i$ and $k$, we have that

\begin{align}\label{bound_one_iteration}
\mathbb{E}\{\mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i+1)}) \mid \boldsymbol{x}_{k+1}^{(i)} \} &\leq \mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i)}) -2\mu \gamma (\mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i)}) - \mathcal{L}_k^*), \quad \text{for } i=0, \dots, \ell-1.
\end{align}

Applying the above inequality recursively for $i=0,1, \dots, \ell-1$ amounts to

\begin{align*}
\mathbb{E}\{\mathcal{L}_k(\boldsymbol{x}_{k+1}) \mid \boldsymbol{x}_{k} \} &\leq \mathcal{L}_k(\boldsymbol{x}_{k}) -2\mu \gamma \sum_{i=0}^\ell \left(\mathcal{L}_k(\boldsymbol{x}_{k+1}^{(i)}) - \mathcal{L}_k^*\right)\\
&\stackrel{(a)}{\leq} \mathcal{L}_k(\boldsymbol{x}_{k}) -2\mu \gamma \ell \left(\mathcal{L}_k(\boldsymbol{x}_{k}) - \mathcal{L}_k^*\right) ,
\end{align*}

where in (a) we have used the (<ref>). By adding and subtracting terms we have

\begin{align*}
\mathbb{E}\left\{\mathcal{L}_{k+1}(\boldsymbol{x}_{k+1}) - \mathcal{L}_{k+1}^* | \boldsymbol{x}_k\right\} = & \mathbb{E}\left\{\mathcal{L}_{k+1}(\boldsymbol{x}_{k+1}) - \mathcal{L}_{k}(\boldsymbol{x}_{k+1}) | \boldsymbol{x}_k\right\} + \mathbb{E}\left\{\mathcal{L}_{k}(\boldsymbol{x}_{k+1}) - \mathcal{L}_{k}(\boldsymbol{x}_{k}) | \boldsymbol{x}_k\right\} \\
&+ \left(\mathcal{L}_k(\boldsymbol{x}_k) - \mathcal{L}_k^*\right) + (\mathcal{L}_k^* - \mathcal{L}_{k+1}^*) \\
&\stackrel{(b)}{\leq} \eta_0 + \eta^* + (1-2\gamma \mu \ell) (\mathcal{L}_k(\boldsymbol{x}_k) - \mathcal{L}_k^*),
\end{align*}

where in (b) we have used the bounds on drift and (<ref>).
By using the tower rule of expectations and induction we have

\begin{align}
\mathbb{E}\{\mathcal{L}_k(\boldsymbol{x}_{k+1}) - \mathcal{L}_{k+1}^*\} &\leq (\eta_0 + \eta^*) \sum_{i=0}^{k} (1-2\gamma \mu \ell)^{k-i} + (1-2\mu\gamma\ell)^k (\mathcal{L}_0(\boldsymbol{x}_0) - \mathcal{L}_0^*) \\
&\leq \frac{\eta_0 + \eta^*}{2\mu \gamma \ell} + (1-2\mu\gamma\ell)^k (\mathcal{L}_0(\boldsymbol{x}_0) - \mathcal{L}_0^*).
\end{align}

Invoking (<ref>) and choosing $\gamma \in (0, \tilde{\gamma})$ completes the proof.

In the case of a pre-trained neural network where the term $\mathcal{L}_0(\boldsymbol{x}_{0}) - \mathcal{L}_{0}^*$ is small (or zero), the slower convergence rate in higher dimensions is not a problem while the asymptotic error bound, the first term in <ref>, is of great interest. By choosing the parameter $\ell$ properly, one can get better bounds based on the speed  of its solver or equivalently the number of iterations at each step.

Under the hypotheses of Theorem <ref> the following holds:

\begin{align}
\limsup_{k\to \infty}{\mathbb{E}\{\mathcal{L}_k(\boldsymbol{x}_{k+1}) - \mathcal{L}_{k+1}^*\} }&\leq \frac{\eta_0 + \eta^*}{2\mu\gamma\ell}.
\end{align}

This in turn results in $\limsup_{k \to \infty} \| \boldsymbol{x}_k - \pi_{\mathcal{X}_{k}^*}(\boldsymbol{x}_{k}) \|^2 \leq \frac{\eta_0 + \eta^*}{\mu^2\gamma\ell}$.

One can choose $\alpha_k = \alpha = \frac{1}{\beta(m+4)}$ to maximise the convergence rate, but this results in maximising the asymptotic optimality gap as well.

§.§.§ Convergence of Proximal Gradient with forward gradients

The proximal gradient algorithm, also known as forward-backward algorithm, applies to the optimisation problems where the cost function is the sum of two functions

\begin{align}\label{sum_problem}
\min_{\boldsymbol{x}\in \mathbb{R}^m}\left\{\mathcal{L}(\boldsymbol{x}) := g(\boldsymbol{x}) + h(\boldsymbol{x})\right\}
\end{align}

which satisfies the following assumption.

Function $g$ is $\beta$-smooth and $h$ is a possibly nonsmooth convex function for all $k$.

The proximal PL-condition in Definition <ref> has been used to analyse the convergence of the generalised proximal gradient algorithm <cit.> where authors have shown that the proximal PL condition is equivalent to the KL condition <cit.>. An important example of cost satisfying the proximal-PL inequality is the $\ell1$-regularized least squares problem.

Let $\mathcal{L}:\mathbb{R}^n \to \mathbb{R}$ be a function that satisfies Proximal-PL condition and Assumption <ref>. Then there exist a constant $\xi > 0$ such that for every $\boldsymbol{x}\in \mathbb{R}^n$ the following holds

\begin{align}
\frac{\xi}{2}\|\boldsymbol{x} - \pi_{\mathcal{X}^*}(\boldsymbol{x}) \|^2 \leq \mathcal{L}(\boldsymbol{x}) - \mathcal{L}^*
\end{align}

Consider the problem (<ref>) where the Assumption <ref> holds and $\mathcal{L}$ satisfies the proximal-PL condition with parameter $\mu$. Then  for proximal gradient algorithm with the step size $\frac{1}{\beta}$ the following holds for some $C>0$

\begin{align*}
\|\boldsymbol{x}_k - \pi_{\mathcal{X}^*}(\boldsymbol{x}_k) \| \leq C (1-\frac{\mu}{\beta})^{\frac{k}{2}}
\end{align*}

The proof is a result of Lemma <ref> and Theorem 5 in <cit.>.

The cost function does not need satisfy proximal-PL condition on the whole space. The condition should only be satisfied in a large enough ball around initial point, $B(\boldsymbol{x}_0, R)$ where $R \geq \frac{\sqrt{\frac{2}{\beta} (\mathcal{L}(\boldsymbol{x}_{0}) - \mathcal{L}^*)}}{1-\sqrt{1-\frac{\mu}{\beta}}}$.

We are interested in analysing the performance of proximal gradient algorithm in an online setting using forward gradient

\begin{align}\label{Prox_grad_iter}
\boldsymbol{x}_{k+1} = \textit{prox}_{\gamma_k h_k} (\boldsymbol{x}_k - \gamma_k \boldsymbol{v}_k(\boldsymbol{x}_k, U_k)) := \argmin_{\boldsymbol{y}} h_k(\boldsymbol{y}) + \frac{1}{2\gamma_k}\|\boldsymbol{y} - (\boldsymbol{x}_k - \gamma_k \boldsymbol{v}_k(\boldsymbol{x}_k, U_k)) \|^2
\end{align}

where $\boldsymbol{v}_k(\boldsymbol{x}, U)$ is the forward gradient of function $g_k$ at point $\boldsymbol{x}$. We aim at analysing the algorithm above in terms of expected optimality gap.

The functions $g_k$ has bounded gradients and functions $h_k$ sub-gradients for each $k$, i.e. $\| \nabla g_k(\boldsymbol{x})\| \leq c_1$ and $\| \partial h_k(\boldsymbol{x}) \| \leq c_2$.

We analyse the performance of the algorithm above, (<ref>), under the proximal PL assumption for time-varying objective functions. We show that the tracking error of the algorithm linearly converges to a neighbourhood of zero, where the asymptotic error depend on the dimension.

Assume that $\mathcal{L}_k$ is sum of a $\beta-$smooth function $g_k$ and a nonsmooth convex function $h_k$  which satisfies Proximal-PL condition with parameter $\mu$ and Assumptions <ref> and <ref> hold. If we use proximal algorithm with forward gradient and constant step-size $\gamma_k = \frac{1}{\beta}$ and random direction $U_k$ is chosen from a standard normal distribution, then

\begin{align*}
\mathbb{E}\{\| \boldsymbol{x}_{k+1} - \pi_{\mathcal{X}_{k+1}^*}(\boldsymbol{x}_{k+1}) \|^2\} &\leq   \frac{2}{ \xi (1-(1-\frac{\mu}{\beta})^{\ell})} \bigg(\eta_0 + \eta^* + 2G_1\sqrt{(m+3)}  \bigg) + \frac{2}{\xi}(1-\frac{\mu}{\beta})^{k\ell} (\mathcal{L}_{0}(\boldsymbol{x}_{0})-\mathcal{L}_{0}^*)
\end{align*}

where $\xi$ is as in Lemma <ref> and $G_1 = \frac{2c_1 (c_1 + c_2)}{\beta}$ with $c_1$ and $c_2$ as in Assumption <ref>.

Under the hypotheses of Theorem <ref> the following holds:

\begin{align}
\limsup_{k\to \infty} \mathbb{E} \bigg\{ \mathcal{L}_{k+1}(\boldsymbol{x}_{k+1})-\mathcal{L}_{k+1}^* \bigg\} \leq \frac{1}{1-(1-\frac{\mu}{\beta})^{\ell}} (\eta_0 + \eta^* + G_1\sqrt{(m+3)})
\end{align}

§ ILLUSTRATIVE NUMERICAL RESULTS

The theoretical error bounds of previous sections are illustrated here employing a numerical example. We consider a sequence of time-varying linear least square problems of the form

\begin{align}
\mathcal{L}_k(\boldsymbol{x}) = \frac{1}{2}\|A_k\boldsymbol{x} - \boldsymbol{b}_k \|^2
\end{align}

where $A_k\in \mathbb{R}^{n\times m}$ and $b_k \in \mathbb{R}^n$. The cost function $\mathcal{L}_k$ satisfies PL condition for every choice of $A_k$ and $\boldsymbol{b}_k$. We consider the over-parameterised case $m=60$ and $n = 10$. The vector $\boldsymbol{b}_k$ is generated as $\boldsymbol{b}_{k+1} = \boldsymbol{b}_k + \delta b_k$, where $\delta b_k$ follows a normal distribution $\mathcal{N}(\boldsymbol{0}, 10^{-2}I_n)$. The matrix $A_k$ is generated using its singular value decomposition, $A_k = U\Sigma_k V^T$ where $U \in \mathbb{R}^{n \times r}$ and $V \in \mathbb{R}^{m \times r}$ are orthogonal matrices and $\Sigma_{k+1} = \Sigma_k - 10^{-6}I_r$ that $\Sigma_0 = diag\ \{0.1, 0.2, \dots, 1\}$ and $r=10$. This setting assures that $\sup_{\boldsymbol{x}} \{\mathcal{L}_{k+1}(\boldsymbol{x}) - \mathcal{L}_k(\boldsymbol{x})\}$ is bounded. We have run 50 different experiments with the same starting point and have calculated the average loss. Results are plotted in Figure <ref>, showing the performance of the actual algorithm and the derived theoretical bounds, validating the convergence results.

Performance of online gradient descent algorithm using directional derivatives, and the theoretical bound

§ CONCLUSION

In this paper, we analysed the performance of gradient-based first-order algorithms focusing on faster algorithms for calculating gradients, namely forward gradients. We have exploited the fact that non-linear least square problems in over-parameterised settings will satisfy the PL condition. Based on this observation, we proved that the (proximal-) gradient-based algorithm based on the forward mode of automatic differentiation (forward gradient) with nonconvex objective functions convergences to optimal value function at a linear rate. We have also analysed the convergence of these algorithms in a time-varying setting and proved the linear convergence to the neighbourhood of a global minimiser. This paper gives new insights on using the forward mode of automatic differentiation in problems with limited resources and fast-changing cost functions where calculating full gradients using the backpropagation algorithm is either impossible or inefficient.

§ PROOF OF LEMMA <REF>

Using descent lemma we have

\begin{align}
\mathcal{L}(\boldsymbol{x}_{k+1}) &\leq \mathcal{L}(\boldsymbol{x}_k) + \langle \nabla \mathcal{L}(\boldsymbol{x}_k), \boldsymbol{x}_{k+1} - \boldsymbol{x}_k \rangle + \frac{\beta}{2}\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_k \|^2 \\
& =  \mathcal{L}(\boldsymbol{x}_k) -\alpha_k \langle \nabla \mathcal{L}(\boldsymbol{x}_k), \boldsymbol{v}_k(\boldsymbol{x}_k, U_k) \rangle + \frac{\beta}{2}\alpha_k^2\|\boldsymbol{v}_k(\boldsymbol{x}_k, U_k) \|^2
\end{align}

By taking conditional expectation given $\boldsymbol{x}_k$ with respect to $U_k$, and using (<ref>) we obtain

\begin{align*}
\mathbb{E}\{\mathcal{L}(\boldsymbol{x}_{k+1}) \mid \boldsymbol{x}_k \} &\leq \mathcal{L}(\boldsymbol{x}_k) -\alpha_k  \|\nabla \mathcal{L}(\boldsymbol{x}_k) \|^2 +\frac{1}{2}\beta \alpha_k^2(m+4)\|\nabla \mathcal{L}(\boldsymbol{x}_k) \|^2\\
& = \mathcal{L}(\boldsymbol{x}_k) -\underbrace{\alpha_k(1-\frac{\beta}{2}(m+4)\alpha_k)}_{:=\gamma_k}\|\nabla\mathcal{L}(\boldsymbol{x}_k) \|^2.
\end{align*}

We now fix $\alpha_k = \alpha = \frac{1}{\beta(m+4)}$ which implies $\gamma_k = \gamma = \frac{1}{2\beta(m+4)}$. Assuming that $\mathcal{L}$ satisfies $\mu-PL$ condition at $\boldsymbol{x}_k$, we have that

\begin{align}
\mathbb{E}\{\mathcal{L}(\boldsymbol{x}_{k+1}) \mid \boldsymbol{x}_k \} &\leq \mathcal{L}(\boldsymbol{x}_k) -2\mu \gamma (\mathcal{L}(\boldsymbol{x}_k) - \mathcal{L}^*).
\end{align}

By using the tower rule of expectations and induction we have

\begin{align}
\mathbb{E}\{\mathcal{L}(\boldsymbol{x}_{k+1}) - \mathcal{L}^*\} &\leq (\eta_0 + \eta^*)  (1-2\mu\gamma)^k (\mathcal{L}(\boldsymbol{x}_0) - \mathcal{L}_0^*) \\
&\leq  (1-2\mu\gamma)^k (\mathcal{L}(\boldsymbol{x}_0) - \mathcal{L}_0^*).
\end{align}

Invoking  $\gamma = \frac{1}{2\beta(m+4)}$ completes the proof.

§ PROOF OF LEMMA <REF>

Let $\boldsymbol{y} : = \argmin_{\boldsymbol{y}}\left\{ \langle \nabla g(\boldsymbol{x}), \boldsymbol{y}-\boldsymbol{x} \rangle + \frac{\beta}{2}\|\boldsymbol{y}-\boldsymbol{x} \|^2 + h(\boldsymbol{y}) - h(\boldsymbol{x})\right\} = \textit{prox}_{\frac{1}{\beta}h}(\boldsymbol{x} - \frac{1}{\beta} \nabla g(\boldsymbol{x}))$,

\begin{align*}
\mathcal{L}(\boldsymbol{y}) &= g(\boldsymbol{y}) + h(\boldsymbol{y}) +h(\boldsymbol{x}) - h(\boldsymbol{x}) \\
&\leq g(\boldsymbol{x}) + h(\boldsymbol{y}) + \langle \nabla g(\boldsymbol{x}), \boldsymbol{y} - \boldsymbol{x} \rangle + \frac{\beta}{2}\|\boldsymbol{y} - \boldsymbol{x}\|^2 +h(\boldsymbol{x}) - h(\boldsymbol{x}) \\
& \stackrel{(a)}{=} \mathcal{L}(\boldsymbol{x}) -\frac{1}{2\beta} \mathcal{D}_h(\boldsymbol{x}, \beta)
\end{align*}

where in (a) we have used the optimality condition and definition of Proximal-PL condition. Therefore,

\begin{align}
\mathcal{D}_{h}(\boldsymbol{x} , \beta) \leq 2\beta ( \mathcal{L}(\boldsymbol{x}) - \mathcal{L}^* )
\end{align}

Using the optimiality condition we have $\boldsymbol{y} - \boldsymbol{x} = -\frac{1}{\beta} \big( \nabla g(\boldsymbol{x}) + \boldsymbol{s} \big)$ where $\boldsymbol{s}\in \partial h(\boldsymbol{y})$. Plugging this back into (<ref>) gives

\begin{align*}
\mathcal{D}_h(\boldsymbol{x}, \beta) &= -2\beta \bigg(\langle \nabla g(\boldsymbol{x}), \boldsymbol{y} - \boldsymbol{x} \rangle + \frac{\beta}{2}\|\boldsymbol{y} - \boldsymbol{x}\|^2 +h(\boldsymbol{y}) - h(\boldsymbol{x}) \bigg)\\
&= -2\beta \bigg(\langle -\beta(\boldsymbol{y} - \boldsymbol{x}) - \boldsymbol{s}, \boldsymbol{y} - \boldsymbol{x} \rangle + \frac{\beta}{2}\|\boldsymbol{y} - \boldsymbol{x}\|^2 +h(\boldsymbol{y}) - h(\boldsymbol{x}) \bigg)\\
&=-2\beta \bigg( -\frac{\beta}{2}\|\boldsymbol{y} - \boldsymbol{x} \|^2 - \langle \boldsymbol{s}, \boldsymbol{y} - \boldsymbol{x} \rangle  +h(\boldsymbol{y}) - h(\boldsymbol{x}) \bigg)\\
&= \beta^2 \|\boldsymbol{y} - \boldsymbol{x} \|^2 + 2\beta \big(h(\boldsymbol{x}) - h(\boldsymbol{y}) -\langle \boldsymbol{s}, \boldsymbol{x} - \boldsymbol{y} \rangle \big)\\
&\stackrel{(a)}{\geq} \beta^2 \|\boldsymbol{y} - \boldsymbol{x} \|^2
\end{align*}

where in (a) we have used the convexity of $h$.

Using the fact that the conditions Proximal-EB (error bounds) and Proximal-PL are equivalent (see Appendix G in <cit.>) there exist $c>0$ such that

\begin{align}
\|\boldsymbol{x} - \pi_{\mathcal{X}^*}(\boldsymbol{x}) \| \leq c \left\|\boldsymbol{x} - \textit{prox}_{\frac{1}{\beta}h}(\boldsymbol{x} - \frac{1}{\beta} \nabla g(\boldsymbol{x})) \right\|
\end{align}

combining the last three inequalities yields in

\begin{align}
2\beta ( \mathcal{L}(\boldsymbol{x}) - \mathcal{L}^* ) \geq \mathcal{D}_{h}(\boldsymbol{x} , \beta) \geq \beta^2 \|\boldsymbol{y} - \boldsymbol{x} \|^2 \geq \frac{\beta^2}{c^2} \| \boldsymbol{x} - \pi_{\mathcal{X}^*}(\boldsymbol{x}) \|^2
\end{align}

setting $\xi = \frac{\beta}{c^2}$ completes the proof.

§ PROOF OF REMARK <REF>

By using Lipschitz continuity of gradient of $g$ we can see

\begin{align*}
\mathcal{L}(\boldsymbol{x}_{k+1}) &= g(\boldsymbol{x}_{k+1}) + h(\boldsymbol{x}_{k+1}) +h(\boldsymbol{x}_{k}) - h(\boldsymbol{x}_{k}) \\
&\leq g(\boldsymbol{x}_{k}) + \langle \nabla g(\boldsymbol{x}_{k}), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle + \frac{\beta}{2}\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k}\|^2 +h(\boldsymbol{x}_{k}) - h(\boldsymbol{x}_{k}) \\
& \stackrel{(a)}{=} \mathcal{L}(\boldsymbol{x}_{k}) -\frac{1}{2\beta} \mathcal{D}_h(\boldsymbol{x}_{k}, \beta)
\end{align*}

where in (a) we have used the Proximal PL definition (<ref>). Therefore

\begin{align} \label{upper_D}
\mathcal{D}_h(\boldsymbol{x}_{k}, \beta) \leq 2 \beta  \big( \mathcal{L}(\boldsymbol{x}_{k}) -  \mathcal{L}(\boldsymbol{x}_{k+1}) \big) \leq 2 \beta  \big( \mathcal{L}(\boldsymbol{x}_{k}) -  \mathcal{L}^* \big).
\end{align}

Using the optimiality condition for (<ref>) results in that $\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} = -\frac{1}{\beta} \big( \nabla g(\boldsymbol{x}_{k}) + \boldsymbol{s}^{k+1} \big)$ where $\boldsymbol{s}^{k+1}\in \partial h(\boldsymbol{x}_{k+1})$. Plugging this back into (<ref>) we have

\begin{align*}
\mathcal{D}_h(\boldsymbol{x}_{k}, \beta) &= -2\beta \bigg(\langle \nabla g(\boldsymbol{x}_{k}), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle + \frac{\beta}{2}\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k}\|^2 +h(\boldsymbol{x}_{k+1}) - h(\boldsymbol{x}_{k}) \bigg)\\
&= -2\beta \bigg(\langle -\beta(\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k}) - \boldsymbol{s}^{k+1}, \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle + \frac{\beta}{2}\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k}\|^2 +h(\boldsymbol{x}_{k+1}) - h(\boldsymbol{x}_{k}) \bigg)\\
&=-2\beta \bigg( -\frac{\beta}{2}\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \|^2 - \langle \boldsymbol{s}^{k+1}, \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle  +h(\boldsymbol{x}_{k+1}) - h(\boldsymbol{x}_{k}) \bigg)\\
&= \beta^2 \|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \|^2 + 2\beta \big(h(\boldsymbol{x}_{k}) - h(\boldsymbol{x}_{k+1}) -\langle \boldsymbol{s}^{k+1}, \boldsymbol{x}_{k} - \boldsymbol{x}_{k+1} \rangle \big)\\
&\stackrel{(a)}{\geq} \beta^2 \|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \|^2
\end{align*}

where in (a) we have used the convexity of $h$. By combining the last inequality with (<ref>) we have

\begin{align}
\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \| &\leq \sqrt{\frac{2}{\beta}\big( \mathcal{L}(\boldsymbol{x}_{k}) -  \mathcal{L}^* \big)} \\
& \leq \sqrt{\frac{2}{\beta} (\mathcal{L}(\boldsymbol{x}_{0}) - \mathcal{L}^*)}(1-\frac{\mu}{\beta})^{\frac{k}{2}}
\end{align}

which in turn gives a bound on the whole length of the sequence

\begin{align}
\|\boldsymbol{x}_{k+1}- \boldsymbol{x}_0 \| \leq \sum_{i=0}^{k}  \|\boldsymbol{x}_{i+1} - \boldsymbol{x}_{i} \| \leq \frac{\sqrt{\frac{2}{\beta} (\mathcal{L}(\boldsymbol{x}_{0}) - \mathcal{L}^*)}}{1-\sqrt{1-\frac{\mu}{\beta}}}:=R
\end{align}

§ PROOF OF THEOREM <REF>

Using Assumption <ref> we have

\begin{align*}
\mathcal{L}_{k+1}(\boldsymbol{x}_{k+1}) &\leq \eta_0 + g_k(\boldsymbol{x}_{k+1}) + h_k(\boldsymbol{x}_{k+1}) + h_k(\boldsymbol{x}_{k}) - h_k(\boldsymbol{x}_{k})\\
& \stackrel{(a)}{\leq} \eta_0 + g_k(\boldsymbol{x}_{k}) + h_k(\boldsymbol{x}_{k})+ \langle \nabla g_k(\boldsymbol{x}_{k}), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle + \frac{\beta}{2}\|\boldsymbol{x}_{k+1}- \boldsymbol{x}_{k} \|^2 + h_k(\boldsymbol{x}_{k+1})- h_k(\boldsymbol{x}_{k})\\
& = \eta_0 + \mathcal{L}_k(\boldsymbol{x}_{k}) + h_k(\boldsymbol{x}_{k+1})+ \langle \nabla g_k(\boldsymbol{x}_{k}), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle + \frac{\beta}{2}\|\boldsymbol{x}_{k+1}- \boldsymbol{x}_{k} \|^2 - h_k(\boldsymbol{x}_{k})\\
&= \eta_0 + \mathcal{L}_k(\boldsymbol{x}_{k}) + h_k(\boldsymbol{x}_{k+1})+ \langle \boldsymbol{v}_k(\boldsymbol{x}_k, U_k) +\nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle \\
&+ \frac{\beta}{2}\|\boldsymbol{x}_{k+1}- \boldsymbol{x}_{k} \|^2 - h_k(\boldsymbol{x}_{k})\\
&\stackrel{(b)}{=} \eta_0 + \mathcal{L}_k(\boldsymbol{x}_{k}) - \frac{1}{2\beta} \widetilde{\mathcal{D}}_{h_k}(\boldsymbol{x}_k, \beta, U_k) +
\langle \nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle
\end{align*}

where in (a) we have used the smoothness of $g_k$ and in (b) we have used the following definition

\begin{align}
\widetilde{\mathcal{D}}_h(\boldsymbol{x}, \alpha, U): = -2\alpha \min_{\boldsymbol{y}} \left\{ \langle \boldsymbol{v}({\boldsymbol{x}}, U), \boldsymbol{y}-\boldsymbol{x} \rangle + \frac{\alpha}{2}\|\boldsymbol{y}-\boldsymbol{x} \|^2 + h(\boldsymbol{y}) - h(\boldsymbol{x})\right\}.
\end{align}

By adding and subtracting term we have

\begin{align*}
\mathcal{L}_{k+1}(\boldsymbol{x}_{k+1})-\mathcal{L}_{k+1}^* &\leq \eta_0 + \eta^* + \mathcal{L}_k(\boldsymbol{x}_{k}) - \mathcal{L}_k^* - \frac{1}{2\beta} {\mathcal{D}}_{h_k}(\boldsymbol{x}_k, \beta)  \\
&+ \frac{1}{2\beta} {\mathcal{D}}_{h_k}(\boldsymbol{x}_k, \beta) - \frac{1}{2\beta} \widetilde{\mathcal{D}}_{h_k}(\boldsymbol{x}_k, \beta, U_k) + \langle \nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle \\
&\stackrel{(a)}{\leq} \eta_0 + \eta^* + (1-\frac{\mu}{\beta})(\mathcal{L}_k(\boldsymbol{x}_k) - \mathcal{L}_k^*) + \frac{1}{2\beta}\left(\widetilde{\mathcal{D}}_{h_k}(\boldsymbol{x}_k, \beta, U_k) - \mathcal{D}_{h_k}(\boldsymbol{x}_k, \beta) \right) \\
& + \langle \nabla g_k(\boldsymbol{x}_{k})  - v(\boldsymbol{x}_k, U_k), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle \\
&\stackrel{(b)}{\leq} \eta_0 + \eta^* + (1-\frac{\mu}{\beta})(\mathcal{L}_k(\boldsymbol{x}_k) - \mathcal{L}_k^*)\\
&+ \min_{\boldsymbol{y}} \left\{ \langle \nabla g_k(\boldsymbol{x}_k), \boldsymbol{y}-\boldsymbol{x}_k \rangle + \frac{\beta}{2}\|\boldsymbol{y}-\boldsymbol{x}_k \|^2 + h_k(\boldsymbol{y}) - h_k(\boldsymbol{x}_k)\right\}\\
& - \min_{\boldsymbol{y}} \left\{ \langle \boldsymbol{v}_k({\boldsymbol{x}_k}, U_k), \boldsymbol{y}-\boldsymbol{x}_k \rangle + \frac{\beta}{2}\|\boldsymbol{y}-\boldsymbol{x}_k \|^2 + h(\boldsymbol{y}) - h(\boldsymbol{x}_k)\right\}\\
& + \langle \nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle\\
&\stackrel{(c)}{\leq} \eta_0 + \eta^* + (1-\frac{\mu}{\beta})(\mathcal{L}_k(\boldsymbol{x}_k) - \mathcal{L}_k^*) +2 \langle \nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k), \boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \rangle
\end{align*}

where in (a) we have used (<ref>), in (b) we have used the definitions for $\widetilde{\mathcal{D}}_{h_k}$ and $\mathcal{D}_{h_k}$, and (c) is true because of the update step and sub-optimality condition.

Note that form optimality condition of (<ref>) we have

\begin{align*}
\boldsymbol{x}_{k+1}- \boldsymbol{x}_{k}\in -\frac{1}{\beta}\left({\boldsymbol{v}_k(\boldsymbol{x}_{k}, U_k) + \partial h_k(\boldsymbol{x}_{k+1})}\right).
\end{align*}

Therefore, $\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \|$ is bounded because of boundedness of (sub)gradients, $\|\boldsymbol{x}_{k+1} - \boldsymbol{x}_{k} \| \leq  \frac{c_1 + c_2}{\beta}$. Using this in the last inequality above results in

\begin{align}
\mathcal{L}_{k+1}(\boldsymbol{x}_{k+1})-\mathcal{L}_{k+1}^* \leq \eta_0 + \eta^* +  (1-\frac{\mu}{\beta})( \mathcal{L}_{k}(\boldsymbol{x}_{k})-\mathcal{L}_{k}^*) + 2\frac{c_1 + c_2}{\beta}\|\nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k) \|.
\end{align}

By taking conditional expectation given $\boldsymbol{x}_k$ with respect to $U_k$, we have

\begin{align*}
\mathbb{E} \bigg\{ \mathcal{L}_{k+1}(\boldsymbol{x}_{k+1})-\mathcal{L}_{k+1}^* \mid \boldsymbol{x}_k \bigg\} &\leq \eta_0 + \eta^* + (1-\frac{\mu}{\beta}) (\mathcal{L}_{k}(\boldsymbol{x}_{k})-\mathcal{L}_{k}^*) \\
&+ 2\frac{c_1 + c_2}{\beta}\ \mathbb{E} \bigg\{ \|\nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k) \| \mid \boldsymbol{x}_k \bigg\}.
\end{align*}

Invoking Lyapunov's inequality, $(\mathbb{E}\|Y\|^r)^{1/r}\leq (\mathbb{E}\|Y\|^s)^{1/s}$ for $r=1, s=2$, we have

\begin{align*}
\mathbb{E} \bigg\{ \mathcal{L}_{k+1}(\boldsymbol{x}_{k+1})-\mathcal{L}_{k+1}^* \mid \boldsymbol{x}_k \bigg\} &\leq \eta_0 + \eta^* + (1-\frac{\mu}{\beta}) (\mathcal{L}_{k}(\boldsymbol{x}_{k})-\mathcal{L}_{k}^*) \\
&+ 2\frac{c_1 + c_2}{\beta}\ \sqrt{ \mathbb{E} \bigg\{ \|\nabla g_k(\boldsymbol{x}_{k})  - \boldsymbol{v}_k(\boldsymbol{x}_k, U_k) \|^2 \mid \boldsymbol{x}_k \bigg\}}\\
& \stackrel{(a)}{=} \eta_0 + \eta_* + (1-\frac{\mu}{\beta}) (\mathcal{L}_{k}(\boldsymbol{x}_{k})-\mathcal{L}_{k}^*) + 2\frac{c_1 + c_2}{\beta}\sqrt{(m+3)} \| \nabla g_k(\boldsymbol{x}_k) \|\\
& \stackrel{(b)}{\leq} (1-\frac{\mu}{\beta}) (\mathcal{L}_{k}(\boldsymbol{x}_{k})-\mathcal{L}_{k}^*) + \eta_0 + \eta^* + 2c_1\frac{c_1 + c_2}{\beta}\sqrt{(m+3)},
\end{align*}

where in (a)  and (b) we have used (<ref>) and the boundedness of gradients, respectively. Using tower rule of expectations

\begin{align*}
\mathbb{E} \bigg\{ \mathcal{L}_{k+1}(\boldsymbol{x}_{k+1})-\mathcal{L}_{k+1}^*  \bigg\} &\leq (1-\frac{\mu}{\beta})^{k} (\mathcal{L}_{0}(\boldsymbol{x}_{0})-\mathcal{L}_{0}^*) \\
&+ \frac{1}{1-(1-\frac{\mu}{\beta})} \bigg(\eta_0 + \eta^* + 2c_1 \frac{c_1 + c_2}{\beta}\sqrt{(m+3)}  \bigg),
\end{align*}

Similarly, with $\ell$ iterations at each step we have

\begin{align*}
\mathbb{E} \bigg\{ \mathcal{L}_{k+1}(\boldsymbol{x}_{k+1})-\mathcal{L}_{k+1}^*  \bigg\} &\leq (1-\frac{\mu}{\beta})^{k\ell} (\mathcal{L}_{0}(\boldsymbol{x}_{0})-\mathcal{L}_{0}^*) \\
&+ \frac{1}{1-(1-\frac{\mu}{\beta})^{\ell}} \bigg(\eta_0 + \eta^* + 2c_1 \frac{c_1 + c_2}{\beta}\sqrt{(m+3)}  \bigg),
\end{align*}

Invoking Lemma <ref> completes the proof.


# Unextendible and uncompletable product bases in every bipartition

Fei Shi<EMAIL_ADDRESS>School of Cyber Security, University of
Science and Technology of China, Hefei, 230026, China  Department of Computer
Science, The University of Hong Kong, Pokfulam Road, Hong Kong, China  Mao-
Sheng Li<EMAIL_ADDRESS>School of Mathematics, South China
University of Technology, Guangzhou 510641, China  Xiande Zhang
<EMAIL_ADDRESS>School of Mathematical Sciences, University of Science
and Technology of China, Hefei, 230026, China  Qi Zhao<EMAIL_ADDRESS>Department of Computer Science, The University of Hong Kong, Pokfulam Road,
Hong Kong, China

###### Abstract

Unextendible product basis is an important object in quantum information
theory and features a broad spectrum of applications, ranging bound entangled
states, quantum nonlocality without entanglement, and Bell inequalities with
no quantum violation. A generalized concept called uncompletable product basis
also attracts much attention. In this paper, we find some unextendible product
bases that are uncompletable product bases in every bipartition, which answers
a 19 year-old open question proposed by DiVincenzo _et al._ [Commun. Math.
Phys. 238, 379 (2003)]. As a consequence, we connect such unextendible product
bases to local hiding of information and give a sufficient condition for the
existence of an unextendible product basis, that is still an unextendible
product basis in every bipartition. Our results advance the understanding of
the geometry of unextendible product bases.

##  I Introduction

An unextendible product basis (UPB) in a multipartite quantum system is an
incomplete orthogonal product basis whose complementary subspace contains no
product state bennett1999unextendible . UPBs have a lot of applications in
quantum information. The mixed state that is proportional to the projector on
the complementary subspace of any UPB is a positive-partial-transpose (PPT)
entangled state. PPT entangled states represent the so-called bound entangled
states from which no pure entanglement can be distilled under local operations
and classical communication (LOCC) bennett1999unextendible . Quantum
nonlocaltiy is another important application. UPBs can not be perfectly
distinguished under local positive operator-valued measures (POVMs) and
classical communication bennett1999unextendible , which shows the phenomenon
of quantum nonlocality without entanglement bennett1999quantum . For perfect
discrimination of UPBs, one can use entanglement resources
cohen2008understanding ; zhang2020locally . Some UPBs are locally irreducible
in every biparition, and showed the phenomenon of strong quantum nonlocality
without entanglement Halder2019Strong ; shi2021strong ; shi2022strongly . UPBs
also can be used to show more nonlocality with less purity
bandyopadhyay2011more . Bell nonlocality is from Bell inequalities, and UPBs
were connected to Bell inequalities with no quantum violation augusiak2011bell
; Augusiak2012tight .

In 2003, DiVincenzo _et al._ generalized the concept of UPBs
divincenzo2003unextendible . An uncompletable product basis (UCPB) in a
multipartite quantum system is an incomplete orthogonal product basis, which
can not be extended to a complete orthogonal product basis
divincenzo2003unextendible . An incomplete orthogonal product basis is a
strongly uncompletable product basis (SUCPB), if it is a UCPB in any locally
extended Hilbert space divincenzo2003unextendible . Actually, the set of all
UPBs is a proper subset of the set of all SUCPBs, and the set of all SUCPBs is
a proper subset of the set of all UCPBs. See also Fig. 2 for the inclusion
relation of these three sets. It is known that UPBs and SUCPBs cannot be
perfectly distinguished under local POVMs and classical communication, and
UCPBs cannot be perfectly distinguished under local projective measurements
and classical communication bennett1999unextendible ;
divincenzo2003unextendible . In divincenzo2003unextendible , DiVincenzo _et
al._ proposed an open question: whether there exists a UPB, which is a UCPB in
every bipartition? This open question exists for 19 years because there are
few constructions of UPBs in multipartite systems, and it is difficult to show
UCPBs in bipartite systems. Such UPBs can be used to understand the geometry
of UPBs. There exists another famous open question for UPBs
demianowicz2018unextendible : can we find a UPB, which is still a UPB in every
bipartition? Such UPBs cannot be perfectly distinguished under local POVMs and
classical communication in every bipartiton bennett1999unextendible , and can
be used to construct genuinely entangled subspaces demianowicz2018unextendible
. Recently, Demianowice showed that such UPBs with the minimum size do not
exist demianowicz2022genuinely . However, the existence of such UPBs is still
unknown.

In this work, we address the 19 year-old open question in
divincenzo2003unextendible , by presenting a UPB with a stronger property,
which is an SUCPB in every bipartition. We also show that such UPBs can be
used for local hiding of information. Tile structures in bipartite systems
provide an efficient method for constructing bipartite UPBs
divincenzo2003unextendible ; shi2020unextendible . We generalize the tile
structures to multipartite systems, and give a sufficient condition for the
existence a UPB that is still a UPB in every bipartition. This sufficient
condition is intuitive, and one can search such UPBs through computer under
this condition.

The rest of this paper is organized as follows. In Sec. II, we introduce the
concepts of UPBs, UCPBs, and SUCPBs. Next, in Sec. III, we find a UPB that is
an SUCPBs in every bipartition for arbitrary three-, and four-partite system.
In Sec. IV, we give a sufficient condition for the existence of a UPB that is
still a UPB in every bipartition. Finally, we conclude in Sec. V.

##  II Preliminaries

In this paper, we do not normalize product states for simplicity. We denote
$\mathbb{Z}_{n}:=\\{0,1,\ldots,n-1\\}$ and $w_{n}:=e^{\frac{2\pi i}{n}}$. For
a matrix $M$, let $\text{sum}(M)$ be the sum of all elements. Assume
$\\{|i\rangle_{A}\\}_{i\in\mathbb{Z}_{m}}$ and
$\\{|j\rangle_{B}\\}_{j\in\mathbb{Z}_{n}}$ are the computational bases of
${\cal H}_{A}$ and ${\cal H}_{B}$, respectively. For any bipartite state
$|\psi\rangle\in{\cal H}_{A}\otimes{\cal H}_{B}$, it can be expressed by

$|\psi\rangle=\sum_{i\in\mathbb{Z}_{m},j\in\mathbb{Z}_{n}}a_{i,j}|i\rangle_{A}|j\rangle_{B}.$
(1)

Then $|\psi\rangle$ corresponds to an $m\times n$ matrix

$M=(a_{i,j})_{i\in\mathbb{Z}_{m},j\in\mathbb{Z}_{n}}.$ (2)

If $\mathop{\rm rank}(M)=1$, then $|\psi\rangle$ is a product state; if
$\mathop{\rm rank}(M)\geq 2$, then $|\psi\rangle$ is an entangled state.
Assume $|\psi_{i}\rangle\in{\cal H}_{A}\otimes{\cal H}_{B}$ corresponds to an
$m\times n$ matrix $M_{i}$ for $i=1,2$, then
$\langle\psi_{1}|\psi_{2}\rangle=\mathop{\rm Tr}({M_{1}^{\dagger}}{M_{2}})$.
Let ${\cal H}=\otimes_{i=1}^{n}{\cal H}_{i}$ be an $n$-partite Hilbert space.
An orthogonal product set (OPS) in ${\cal H}$ is a set of orthogonal product
states, and an orthogonal product basis (OPB) in ${\cal H}$ is an OPS which
spans ${\cal H}$. Given ${\cal H}=\otimes_{i=1}^{n}{\cal H}_{i}$, let ${\cal
H}_{ext}=\otimes_{i=1}^{n}({\cal H}_{i}\oplus{\cal H}_{i}^{\prime})$ be a
locally extended Hilbert space of ${\cal H}$, where $H_{i}^{\prime}$ is a
local extension. Now, we review some definitions.

###### Definition 1

Let ${\cal S}$ be an OPS in ${\cal H}=\otimes_{i=1}^{n}{\cal H}_{i}$. The set
${\cal S}$ spans a subspace ${\cal H}_{{\cal S}}$ in ${\cal H}$, and
$\mathop{\rm Dim}({\cal H}_{{\cal S}})<\mathop{\rm Dim}({\cal H})$. If the
complementary subspace ${\cal H}_{{\cal S}}^{\bot}$ contains no product state,
then ${\cal S}$ is called an _unextendible product basis (UPB)_. If ${\cal S}$
cannot be extended to an OPB in ${\cal H}$, then ${\cal S}$ is called an
_uncompletable product basis (UCPB)_. Moreover, if ${\cal S}$ is a UCPB in any
locally extended Hilbert space ${\cal H}_{ext}=\otimes_{i=1}^{n}({\cal
H}_{i}\oplus{\cal H}_{i}^{\prime})$, then ${\cal S}$ is called a _strongly
uncompletable product basis (SUCPB)_.

From Definition 1, a UPB or an SUCPB must be a UCPB. We can always obtain a
UPB from a UCPB ${\cal S}$, by adding some orthogonal product states to ${\cal
S}$ from ${\cal H}_{{\cal S}}^{\bot}$ till the new OPS is a UPB. Moreover,
UPBs and SUCPBs cannot be perfectly distinguished under local POVMs and
classical communication, and UCPBs can not be perfectly distinguished under
local projective measurements and classical communication
bennett1999unextendible .

For an OPS ${\cal S}=\\{|\psi_{i}\rangle\\}_{i=1}^{s}$ in ${\cal
H}=\otimes_{i=1}^{n}{\cal H}_{i}$ with $\mathop{\rm Dim}({\cal H})=D$ (where
$s<D$), we can define a mixed state that is proportional to the projector on
${\cal H}_{{\cal S}}^{\bot}$,

$\overline{\rho}_{{\cal
S}}=\frac{1}{D-s}\left(\mathbb{I}-\sum_{i=1}^{s}|\psi_{i}\rangle\\!\langle\psi_{i}|\right).$
(3)

Applying partial transposition map to $\overline{\rho}_{{\cal S}}$ in any
bipartition, then we can find that $(\mathbb{I}\otimes
T)\overline{\rho}_{{\cal S}}\geq 0$. It means that $\overline{\rho}_{{\cal
S}}$ has the positive partial transpose (PPT) property in any bipartition. If
${\cal S}$ is a UPB, then $\overline{\rho}_{{\cal S}}$ must be entangled from
the definition. Thus $\overline{\rho}_{{\cal S}}$ is a PPT entangled state,
which is also a bound entangled state (no pure entanglement can be distilled)
bennett1999unextendible ; divincenzo2003unextendible . However, if ${\cal S}$
is a UCPB or an SUCPB, $\overline{\rho}_{{\cal S}}$ is either separable or
entangled bennett1999unextendible ; divincenzo2003unextendible .

It is difficult to show that an OPS is an SUCPB from the definition. This
exists a sufficient condition.

###### Lemma 1

Let ${\cal S}$ be an OPS in ${\cal H}=\otimes_{i=1}^{n}{\cal H}_{i}$. If all
the product states in ${\cal H}_{{\cal S}}^{\bot}$ cannot span ${\cal
H}_{{\cal S}}^{\bot}$, then ${\cal S}$ is an SUCPB.

Proof.  If all the product states in ${\cal H}_{{\cal S}}^{\bot}$ cannot span
${\cal H}_{{\cal S}}^{\bot}$, then $\overline{\rho}_{{\cal S}}$ must be
entangled by Theorem 2(ii) in horodecki1997separability . Further, according
to Proposition 1 in divincenzo2003unextendible , ${\cal S}$ is an SUCPB.

Figure 1: Tile structure with $6$ tiles in
$\mathbb{C}^{3}\otimes\mathbb{C}^{4}$.

By Lemma 1, a UPB must be an SUCPB. However, the converse is not true. We will
give an example of SUCPB, which is not a UPB.

Tile structures can be used to construct UPBs divincenzo2003unextendible ;
halder2019family ; shi2020unextendible . Next, we show that tile structures
can also be used to construct SUCPBs. A tile structure ${\cal T}$ in
$\mathbb{C}^{m}\otimes\mathbb{C}^{n}$ is an $m\times n$ rectangle, which can
be partitioned into $s$ disjoint tiles $\\{t_{i}\\}_{i=1}^{s}$. Each tile
$t_{i}$ is a rectangle. We denote ${\cal T}:=\cup_{i=1}^{s}t_{i}$. For
example, Fig. 1 gives a tile structure ${\cal T}=\cup_{i=1}^{6}t_{i}$ with $6$
tiles in $\mathbb{C}^{3}\otimes\mathbb{C}^{4}$. Any tile $t_{i}$ of ${\cal
T}=\cup_{i=1}^{s}t_{i}$ has row coordinates
$\\{p_{0},p_{1},\ldots,p_{k-1}\\}_{A}$ and column coordinates
$\\{q_{0},q_{1},\ldots,q_{\ell-1}\\}_{B}$, and we denote it as
$t_{i}=\\{p_{0},p_{1},\ldots,p_{k-1}\\}_{A}\times\\{q_{0},q_{1},\ldots,q_{\ell-1}\\}_{B}$,
where $\\{p_{0},p_{1},\ldots,p_{k-1}\\}$ and
$\\{q_{0},q_{1},\ldots,q_{\ell-1}\\}$ are subsets of $\mathbb{Z}_{m}$ and
$\mathbb{Z}_{n}$, respectively. For tile $t_{i}$, we can construct an OPS of
size $k\ell$ in $\mathbb{C}^{m}\otimes\mathbb{C}^{n}$,

$\displaystyle{\cal A}_{i}=\\{$
$\displaystyle|\psi_{i}{(a,b)}\rangle:=\left(\sum_{e\in\mathbb{Z}_{k}}m_{a,e}|p_{e}\rangle\right)_{A}\left(\sum_{e\in\mathbb{Z}_{\ell}}n_{b,e}|q_{e}\rangle\right)_{B}$
(4) $\displaystyle\mid(a,b)\in\mathbb{Z}_{k}\times\mathbb{Z}_{\ell}\\}.$

Here the coefficient matrix $M=(m_{a,e})_{a,e\in\mathbb{Z}_{k}}$ is a $k\times
k$ row orthogonal matrix (row vectors are mutually orthogonal), and
$m_{0,e}=1$ for $e\in\mathbb{Z}_{k}$, and the coefficient matrix
$N=(n_{b,e})_{b,e\in\mathbb{Z}_{\ell}}$ is an $\ell\times\ell$ row orthogonal
matrix, and $n_{0,e}=1$ for $e\in\mathbb{Z}_{\ell}$. For example, we can
choose $M=(w_{k}^{ae})_{a,e\in\mathbb{Z}_{k}}$, and
$N=(w_{\ell}^{be})_{b,e\in\mathbb{Z}_{\ell}}$. Since those tiles in ${\cal
T}=\cup_{i=1}^{s}t_{i}$ are disjoint, we can obtain an OPB

${\cal B}:=\cup_{i=1}^{s}{\cal A}_{i}$ (5)

in $\mathbb{C}^{m}\otimes\mathbb{C}^{n}$. Further, we define the “stopper”
state as

$|S\rangle=\left(\sum_{i\in\mathbb{Z}_{m}}|i\rangle\right)_{A}\left(\sum_{j\in\mathbb{Z}_{n}}|j\rangle\right)_{B}.$
(6)

We mainly consider the following OPS,

${\cal S}:=\cup_{i=1}^{s}({\cal
A}_{i}\setminus\\{|\psi_{i}{(0,0)}\rangle\\})\cup\\{|S\rangle\\}.$ (7)

For example, using the tile structure ${\cal T}=\cup_{i=1}^{6}t_{i}$ in Fig.
1, we obtain an OPB ${\cal B}=\cup_{i=1}^{6}{\cal A}_{i}$ in
$\mathbb{C}^{3}\otimes\mathbb{C}^{4}$, where

$\displaystyle{\cal A}_{1}$
$\displaystyle=\\{|\psi_{1}{(0,b)}\rangle=|0\rangle_{A}(|0\rangle+(-1)^{b}|1\rangle)_{B}\mid
b\in\mathbb{Z}_{2}\\},$ (8) $\displaystyle{\cal A}_{2}$
$\displaystyle=\\{|\psi_{2}{(0,0)}\rangle=|0\rangle_{A}|2\rangle_{B}\\},$
$\displaystyle{\cal A}_{3}$
$\displaystyle=\\{|\psi_{3}{(a,0)}\rangle=(|0\rangle+(-1)^{a}|1\rangle)_{A}|3\rangle_{B}\mid
a\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{4}$
$\displaystyle=\\{|\psi_{4}{(0,b)}\rangle=|2\rangle_{A}(|1\rangle+w_{3}^{b}|2\rangle+w_{3}^{2b}|3\rangle)_{B}\mid
b\in\mathbb{Z}_{3}\\},$ $\displaystyle{\cal A}_{5}$
$\displaystyle=\\{|\psi_{5}{(a,0)}\rangle=(|1\rangle+(-1)^{a}|2\rangle)_{A}|0\rangle_{B}\mid
a\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{6}$
$\displaystyle=\\{|\psi_{6}{(0,b)}\rangle=|1\rangle_{A}(|1\rangle+(-1)^{b}|2\rangle)_{B}\mid
b\in\mathbb{Z}_{2}\\}.$

The “stopper” state is

$|S\rangle=(|0\rangle+|1\rangle+|2\rangle)_{A}(|0\rangle+|1\rangle+|2\rangle+|3\rangle)_{B}.$
(9)

Next, we show that

${\cal S}:=\cup_{i=1}^{6}({\cal
A}_{i}\setminus\\{|\psi_{i}{(0,0)}\rangle\\})\cup\\{|S\rangle\\}$ (10)

is an SUCPB in $\mathbb{C}^{3}\otimes\mathbb{C}^{4}$.

###### Example 1

In $\mathbb{C}^{3}\otimes\mathbb{C}^{4}$, the OPS ${\cal S}$ given by Eq. (10)
is an SUCPB.

Proof.  Let ${\cal S}_{1}:=\cup_{i=1}^{6}({\cal
A}_{i}\setminus\\{|\psi_{i}{(0,0)}\rangle\\})$ and ${\cal
S}_{2}:=\cup_{i=1}^{6}|\psi_{i}{(0,0)}\rangle$. We know that ${\cal
S}_{1}\cup{\cal S}_{2}$ is an OPB in $\mathbb{C}^{3}\otimes\mathbb{C}^{4}$.
Since ${\cal H}_{{\cal S}_{1}}\subset{\cal H}_{{\cal S}}$, it implies ${\cal
H}_{{\cal S}}^{\bot}\subset{\cal H}_{{\cal S}_{1}}^{\bot}={\cal H}_{{\cal
S}_{2}}$. Then for any product state $|\psi\rangle\in{\cal H}_{{\cal
S}}^{\bot}$, there exists $a_{i}\in\mathbb{C}$ for $1\leq i\leq 6$, such that

$|\psi\rangle=\sum_{i=1}^{6}a_{i}|\psi_{i}{(0,0)}\rangle.$

Next, $|\psi\rangle$ corresponds to a $3\times 4$ matrix,

$M=\begin{pmatrix}a_{5}&a_{4}&a_{4}&a_{4}\\\ a_{5}&a_{6}&a_{6}&a_{3}\\\
a_{1}&a_{1}&a_{2}&a_{3}\\\ \end{pmatrix}.$

Note that $M$ has a similar structure to the tile structure in Fig. 1. The
“stopper” state $|S\rangle$ corresponds to a all-ones matrix $J$, where every
element is equal to one. Since $|\psi\rangle$ is a product state and $\langle
S|\psi\rangle=0$, we have $\mathop{\rm rank}(M)=1$ and $\text{sum}(M)=0$. This
is only possible for

$M=\begin{pmatrix}0&0&0&0\\\ 0&0&0&0\\\ a_{1}&a_{1}&a_{2}&0\\\
\end{pmatrix},\quad 2a_{1}+a_{2}=0.$

It means that ${\cal H}_{{\cal S}}^{\bot}$ contains only one product state
$|0\rangle(|0\rangle+|1\rangle-2|2\rangle)$. Since $\mathop{\rm Dim}({\cal
H}_{{\cal S}}^{\bot})=5$, the OPS ${\cal S}$ is an SUCPB by Lemma 1.

Since there exists a product state
$|\psi\rangle=|0\rangle(|0\rangle+|1\rangle-2|2\rangle)\in{\cal H}_{{\cal
S}}^{\bot}$, ${\cal S}$ is not a UPB. However, if we add $|\psi\rangle$ to
${\cal S}$, then ${\cal S}^{\prime}={\cal S}\cup\\{|\psi\rangle\\}$ must be a
UPB in $\mathbb{C}^{3}\otimes\mathbb{C}^{4}$. In fact, for any product state
$|\phi\rangle\in{\cal H}_{{\cal S}^{\prime}}^{\bot}$, $|\phi\rangle$
corresponds to a $3\times 4$ matrix,

$M=\begin{pmatrix}a_{5}&a_{4}&a_{4}&a_{4}\\\ a_{5}&a_{6}&a_{6}&a_{3}\\\
a_{1}&a_{1}&a_{2}&a_{3}\\\ \end{pmatrix},$

where $2a_{1}=a_{2}$, $\mathop{\rm rank}(M)=1$ and $\text{sum}(M)=0$. Such a
matrix $M$ does not exist. From the above discussion, we can obtain a
sufficient condition for the construction of UPBs by tile structures.

###### Lemma 2

For a tile structure ${\cal T}=\cup_{i=1}^{s}t_{i}$ ($s\geq 5$) in
$\mathbb{C}^{m}\otimes\mathbb{C}^{n}$, if any $r$ ($2\leq r\leq s-1$) tiles
cannot form a rectangle, then the OPS ${\cal S}$ given by Eq. (7) is a UPB in
$\mathbb{C}^{m}\otimes\mathbb{C}^{n}$.

The tile structure in Lemma 2 is the U-tile structure proposed in Ref.
shi2020unextendible . In Ref. bennett1999unextendible , the authors gave a
UCPB, which is not an SUCPB. Let ${\cal D}(\text{UPB})$ be the set of all
UPBs; ${\cal D}(\text{SUCPB})$ be the set of all SUCPBs; and ${\cal
D}(\text{UCPB})$ be the set of all UCPBs. Then following set inclusion
relation is obtained,

${\cal D}(\text{UPB})\subsetneq{\cal D}(\text{SUCPB})\subsetneq{\cal
D}(\text{UCPB}),$

See also Fig. 2 for the inclusion relation of these three sets.

Figure 2: A set inclusion relation among the set of all UPBs ${\cal
D}(\text{UPB})$, the set of all SUCPBs ${\cal D}(\text{SUCPB})$, and the set
of all UCPBs ${\cal D}(\text{UCPB})$. The set ${\cal D}(\text{UPB})$ is a
proper subset of ${\cal D}(\text{SUCPB})$, and ${\cal D}(\text{SUCPB})$ is a
proper subset of ${\cal D}(\text{UCPB})$.

In Ref. divincenzo2003unextendible , the authors proposed an open question:
can we find a UPB which is a UCPB in every bipartition? We will give a
positive answer, by showing a stronger UPB, which is an SUCPB in every
bipartition.

##  III The existence of a UPB that is an SUCPB in every bipartition

In this section, we show that there exists a UPB which is an SUCPB in every
bipartition in any three, and four-partite system. Since any OPS in
$\mathbb{C}^{2}\otimes\mathbb{C}^{n}$ can be extended to an OPB
bennett1999unextendible ; divincenzo2003unextendible , the minimum system for
the existence of such UPBs is
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$. The following UPB
in $\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$ is from
agrawal2019genuinely , which is constructed from the tile structure in
tripartite system (we will introduce tile structures in multipartite systems
in Sec. IV). Consider an OPB $\cup_{i=1}^{9}{\cal A}_{i}$ in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$,

$\displaystyle{\cal A}_{1}:=$
$\displaystyle\\{|\psi_{1}(i,j)\rangle=|\xi_{i}\rangle_{A}|0\rangle_{B}|\eta_{k}\rangle_{C}\mid(i,j)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
(11) $\displaystyle{\cal A}_{2}:=$
$\displaystyle\\{|\psi_{2}(i,j)\rangle=|\xi_{i}\rangle_{A}|\eta_{j}\rangle_{B}|2\rangle_{C}\mid(i,j)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{3}:=$
$\displaystyle\\{|\psi_{3}(i,j)\rangle=|2\rangle_{A}|\xi_{i}\rangle_{B}|\eta_{j}\rangle_{C}\mid(i,j)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{4}:=$
$\displaystyle\\{|\psi_{4}(i,j)\rangle=|\eta_{i}\rangle_{A}|2\rangle_{B}|\xi_{j}\rangle_{C}\mid(i,j)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{5}:=$
$\displaystyle\\{|\psi_{5}(i,j)\rangle=|\eta_{i}\rangle_{A}|\xi_{j}\rangle_{B}|0\rangle_{C}\mid(i,j)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{6}:=$
$\displaystyle\\{|\psi_{6}(i,j)\rangle=|0\rangle_{A}|\eta_{i}\rangle_{B}|\xi_{j}\rangle_{C}\mid(i,j)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{7}:=$
$\displaystyle\\{|0\rangle_{A}|0\rangle_{B}|0\rangle_{C}\\},$
$\displaystyle{\cal A}_{8}:=$
$\displaystyle\\{|1\rangle_{A}|1\rangle_{B}|1\rangle_{C}\\},$
$\displaystyle{\cal A}_{9}:=$
$\displaystyle\\{|2\rangle_{A}|2\rangle_{B}|2\rangle_{C}\\},$

where $|\eta_{s}\rangle_{X}=|0\rangle_{X}+(-1)^{s}|1\rangle_{X}$,
$|\xi_{s}\rangle_{X}=|1\rangle_{X}+(-1)^{s}|2\rangle_{X}$ for
$s\in\mathbb{Z}_{2}$, and $X\in\\{A,B,C\\}$. The “stopper” state is,

$|S\rangle=(|0\rangle+|1\rangle+|2\rangle)_{A}(|0\rangle+|1\rangle+|2\rangle)_{B}(|0\rangle+|1\rangle+|2\rangle)_{C}.$
(12)

Then

${\cal U}:=\cup_{i=1}^{6}({\cal
A}_{i}\setminus\\{|\psi_{i}(0,0)\rangle\\})\cup|S\rangle$ (13)

is a UPB in $\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$
agrawal2019genuinely . Now, we have the following lemma.

Figure 3: The corresponding tile structure in
$\mathbb{C}^{3}\otimes\mathbb{C}^{9}$ of the OPB $\cup_{i=1}^{9}{\cal A}_{i}$
(Eq. (11)) in $A|BC$ bipartition.

###### Lemma 3

In $\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$, the UPB ${\cal
U}$ given by Eq. (13) is an SUCPB in every bipartition.

Proof.  First, we consider the bipartition $A|BC$. The OPB
$\cup_{i=1}^{9}{\cal A}_{i}$ given by Eq. (11) in $A|BC$ bipartition
corresponds the tile structure in Fig. 3. Next, we show that the OPS ${\cal
U}_{A|BC}$ is an SUCPB in $\mathbb{C}^{3}\otimes\mathbb{C}^{9}$. For the same
discussion as Example 1, we can assume that $|\psi\rangle\in{\cal H}_{{\cal
U}_{A|BC}}^{\bot}$ is a product state. By Fig. 3, $|\psi\rangle$ corresponds
to a $3\times 9$ matrix

$M=\begin{pmatrix}a_{1}&a_{1}&a_{2}&a_{2}&a_{3}&a_{3}&a_{3}&a_{3}&a_{9}\\\
a_{1}&a_{1}&a_{2}&a_{2}&a_{8}&a_{5}&a_{5}&a_{4}&a_{4}\\\
a_{7}&a_{6}&a_{6}&a_{6}&a_{6}&a_{5}&a_{5}&a_{4}&a_{4}\\\ \end{pmatrix},$

where $a_{i}\in\mathbb{C}$ for $1\leq i\leq 9$, $\mathop{\rm rank}(M)=1$, and
$\text{sum}(M)=0$. There are only four cases,

  1. (i)

$a_{1}+a_{2}=0$, $a_{i}=0$ for $1\leq i\leq 9$ and $i\neq 1,2$;

  2. (ii)

$4a_{3}+a_{9}=0$, $a_{i}=0$ for $1\leq i\leq 9$ and $i\neq 3,9$;

  3. (iii)

$a_{4}+a_{5}=0$, $a_{i}=0$ for $1\leq i\leq 9$ and $i\neq 4,5$;

  4. (iv)

$4a_{6}+a_{7}=0$, $a_{i}=0$ for $1\leq i\leq 9$ and $i\neq 6,7$.

It means that there are only four product states in ${\cal H}_{{\cal
U}_{A|BC}}^{\bot}$:
$(|1\rangle+|2\rangle)_{A}(|00\rangle+|01\rangle-|02\rangle-|12\rangle)_{BC}$,
$|2\rangle_{A}(|11\rangle+|10\rangle+|20\rangle+|21\rangle-4|22\rangle)_{BC}$,
$(|0\rangle+|1\rangle)_{A}(|10\rangle+|20\rangle-|21\rangle-|22\rangle)_{BC}$,
and
$|0\rangle_{A}(|01\rangle+|02\rangle+|12\rangle+|11\rangle-4|00\rangle)_{BC}$.
Since $\mathop{\rm Dim}({\cal H}_{{\cal U}_{A|BC}}^{\bot})=8$, the OPS ${\cal
U}_{A|BC}$ is an SUCPB by Lemma 1.

Further, since the OPB $\cup_{i=1}^{9}{\cal A}_{i}$ given by Eq. (11) in any
bipartition of $\\{A|BC,B|AC,C|AB\\}$ corresponds to a similar tile structure
in Fig. 3, we obtain that ${\cal U}_{A|BC}$, ${\cal U}_{B|AC}$, and ${\cal
U}_{C|AB}$ are all SUCPBs. Thus ${\cal U}$ is an SUCPB in every bipartition.

A similar construction of UPB in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\mathbb{C}^{d_{3}}$ for
$d_{1},d_{2},d_{3}\geq 3$ was given in shi2021strong . For the same discussion
as Lemma 3, we have the following theorem.

###### Theorem 1

In $\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\mathbb{C}^{d_{3}}$,
$d_{1},d_{2},d_{3}\geq 3$, there exists a UPB which is an SUCPB in every
bipartition.

Next, we consider the four-partite UPB. The following UPB in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$
was given in shi2021strong , which is constructed from the tile structure in
four-partite system. Consider an OPB in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$,

$\displaystyle{\cal A}_{1}:=\\{$
$\displaystyle|\psi_{1}(i,j,k)\rangle=|\xi_{i}\rangle_{A}|\eta_{j}\rangle_{B}|0\rangle_{C}|\xi_{k}\rangle_{D}$
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{2}:=\\{$
$\displaystyle|\psi_{2}(i,j,k)\rangle=|\xi_{i}\rangle_{A}|2\rangle_{B}|\eta_{j}\rangle_{C}|\eta_{k}\rangle_{D}$
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{3}:=\\{$
$\displaystyle|\psi_{3}(i,j,k)\rangle=|\xi_{i}\rangle_{A}|\xi_{j}\rangle_{B}|\xi_{k}\rangle_{C}|2\rangle_{D}$
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{4}:=\\{$
$\displaystyle|\psi_{4}(i,j,k)\rangle=|\xi_{i}\rangle_{A}|2\rangle_{B}|0\rangle_{C}|2\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{5}:=\\{$
$\displaystyle|\psi_{5}(i,j,k)\rangle=|2\rangle_{A}|\eta_{i}\rangle_{B}|\xi_{j}\rangle_{C}|\eta_{k}\rangle_{D}$
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{6}:=\\{$
$\displaystyle|\psi_{6}(i,j,k)\rangle=|2\rangle_{A}|\eta_{i}\rangle_{B}|0\rangle_{C}|0\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{7}:=\\{$
$\displaystyle|\psi_{7}(i,j,k)\rangle=|2\rangle_{A}|0\rangle_{B}|\xi_{i}\rangle_{C}|2\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{8}:=\\{$
$\displaystyle|\psi_{8}(i,j,k)\rangle=|2\rangle_{A}|2\rangle_{B}|2\rangle_{C}|\eta_{i}\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{9}:=\\{$
$\displaystyle|\psi_{9}(i,j,k)\rangle=|\eta_{i}\rangle_{A}|\xi_{j}\rangle_{B}|2\rangle_{C}|\eta_{k}\rangle_{D}$
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{10}:=\\{$
$\displaystyle|\psi_{10}(i,j,k)\rangle=|\eta_{i}\rangle_{A}|0\rangle_{B}|\xi_{j}\rangle_{C}|\xi_{k}\rangle_{D}$
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{11}:=\\{$
$\displaystyle|\psi_{11}(i,j,k)\rangle=|\eta_{i}\rangle_{A}|\eta_{j}\rangle_{B}|\eta_{k}\rangle_{C}|0\rangle_{D}$
(14)
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{12}:=\\{$
$\displaystyle|\psi_{12}(i,j,k)\rangle=|\eta_{i}\rangle_{A}|0\rangle_{B}|2\rangle_{C}|0\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{13}:=\\{$
$\displaystyle|\psi_{13}(i,j,k)\rangle=|0\rangle_{A}|\xi_{i}\rangle_{B}|\eta_{j}\rangle_{C}|\xi_{k}\rangle_{D}$
$\displaystyle\mid(i,j,k)\in\mathbb{Z}_{2}\times\mathbb{Z}_{2}\times\mathbb{Z}_{2}\\},$
$\displaystyle{\cal A}_{14}:=\\{$
$\displaystyle|\psi_{14}(i,j,k)\rangle=|0\rangle_{A}|\xi_{i}\rangle_{B}|2\rangle_{C}|2\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{15}:=\\{$
$\displaystyle|\psi_{15}(i,j,k)\rangle=|0\rangle_{A}|2\rangle_{B}|\eta_{i}\rangle_{C}|0\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{16}:=\\{$
$\displaystyle|\psi_{16}(i,j,k)\rangle=|0\rangle_{A}|0\rangle_{B}|0\rangle_{C}|\xi_{i}\rangle_{D}\mid
i\in\mathbb{Z}_{2}\\},$ $\displaystyle{\cal A}_{17}:=\\{$
$\displaystyle|1\rangle_{A}|1\rangle_{B}|1\rangle_{C}|1\rangle_{D}\\},$

where $|\eta_{s}\rangle_{X}=|0\rangle_{X}+(-1)^{s}|1\rangle_{X}$,
$|\xi_{s}\rangle_{X}=|1\rangle_{X}+(-1)^{s}|2\rangle_{X}$ for
$s\in\mathbb{Z}_{2}$, and $X\in\\{A,B,C,D\\}$. The “stopper” state is,

$\displaystyle|S\rangle=$
$\displaystyle(|0\rangle+|1\rangle+|2\rangle)_{A}(|0\rangle+|1\rangle+|2\rangle)_{B}(|0\rangle+|1\rangle+|2\rangle)_{C}$
(15) $\displaystyle(|0\rangle+|1\rangle+|2\rangle)_{D}.$

Then

${\cal V}:=\cup_{i=1}^{16}({\cal
A}_{i}\setminus\\{|\psi_{i}(0,0,0)\rangle\\})\cup\\{|S\rangle\\}$ (16)

is a UPB in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$
agrawal2019genuinely . We can show that this UPB is an SUCPB in every
bipartition.

###### Lemma 4

In
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$,
the UPB ${\cal V}$ given by Eq. (16) is an SUCPB in every bipartition.

Proof.  We need to consider the bipartition set
$\\{A|BCD,B|ACD,C|ABD,D|ABC,AB|CD,AC|BD,$ $AD|BC\\}$. Since the OPB
$\cup_{i=1}^{17}{\cal A}_{i}$ given by Eqs. (14) in any bipartition of
$\\{A|BCD,B|ACD,C|ABD,D|ABC\\}$ (or $\\{AB|CD,AC|BD,$ $AD|BC\\}$) has a
similar structure, we only need to consider ${\cal V}_{A|BCD}$ and ${\cal
V}_{AB|CD}$.

For ${\cal V}_{A|BCD}$, we assume that $|\psi\rangle\in{\cal H}_{{\cal
V}_{A|BCD}}^{\bot}$ is a product state, then $|\psi\rangle$ corresponds to a
$3\times 27$ matrix,

where $a_{i}\in\mathbb{C}$ for $1\leq i\leq 17$, $\mathop{\rm rank}(M)=1$, and
$\text{sum}(M)=0$. There are only four cases,

  1. (i)

$4a_{1}+4a_{2}+4a_{3}+a_{4}=0$, and $a_{i}=0$ for $1\leq i\leq 17$ and $i\neq
1,2,3,4$;

  2. (ii)

$4a_{5}+a_{6}+a_{7}+a_{8}=0$, and $a_{i}=0$ for $1\leq i\leq 17$ and $i\neq
5,6,7,8$;

  3. (iii)

$4a_{9}+4a_{10}+4a_{11}+a_{12}=0$, and $a_{i}=0$ for $1\leq i\leq 17$ and
$i\neq 9,10,11,12$;

  4. (iv)

$4a_{13}+a_{14}+a_{15}+a_{16}=0$, and $a_{i}=0$ for $1\leq i\leq 17$ and
$i\neq 13,14,15,16$.

Then $|\psi\rangle$ must belong to one of the four subspaces,

  1. (i)

$O_{1}=\\{(|1\rangle+|2\rangle)_{A}(a_{1}(|000\rangle+|002\rangle+|101\rangle+|102\rangle)+a_{2}(|200\rangle+|210\rangle+|201\rangle+|211\rangle)+a_{3}(|122\rangle+|222\rangle+|112\rangle+|212\rangle)+a_{4}|202\rangle)_{BCD}\mid
4a_{1}+4a_{2}+4a_{3}+a_{4}=0\\}$;

  2. (ii)

$O_{2}=\\{(|2\rangle)_{A}(a_{5}(|111\rangle+|020\rangle+|010\rangle+|110\rangle+|011\rangle+|021\rangle+|120\rangle+|121\rangle)+a_{6}(|000\rangle+|100\rangle)+a_{7}(|012\rangle+|022\rangle)+a_{8}(|220\rangle+|221\rangle))_{BCD}\mid
4a_{5}+a_{6}+a_{7}+a_{8}=0\\}$;

  3. (iii)

$O_{3}=\\{(|0\rangle+|1\rangle)_{A}(a_{9}(|221\rangle+|220\rangle+|121\rangle+|120\rangle)+a_{10}(|022\rangle+|012\rangle+|021\rangle+|011\rangle)+a_{11}(|100\rangle+|000\rangle+|110\rangle+|010\rangle)+a_{12}|020\rangle)_{BCD}\mid
4a_{9}+4a_{10}+4a_{11}+a_{12}=0\\}$;

  4. (iv)

$O_{4}=\\{(|0\rangle)_{A}(a_{13}(|111\rangle+|202\rangle+|212\rangle+|112\rangle+|211\rangle+|201\rangle+|102\rangle+|101\rangle)+a_{14}(|222\rangle+|122\rangle)+a_{15}(|210\rangle+|200\rangle)+a_{16}(|002\rangle+|000\rangle))_{BCD}\mid
4a_{13}+a_{14}+a_{15}+a_{16}=0\\}$,

where $\mathop{\rm Dim}(O_{i})=3$ for $1\leq i\leq 4$, and $O_{i}\bot O_{j}$
for $1\leq i\neq j\leq 4$. Then $\mathop{\rm
Dim}(O_{1}+O_{2}+O_{3}+O_{4})=12$. Since $\mathop{\rm Dim}({\cal H}_{{\cal
V}_{A|BCD}}^{\bot})=16$, the OPS ${\cal V}_{A|BCD}$ is an SUCPB by Lemma 1.

For ${\cal V}_{AB|CD}$, we assume that $|\phi\rangle\in{\cal H}_{{\cal
V}_{AB|CD}}^{\bot}$ is a product state, then $|\phi\rangle$ corresponds to a
$9\times 9$ matrix,

where $b_{i}\in\mathbb{C}$ for $1\leq i\leq 17$, $\mathop{\rm rank}(N)=1$, and
$\text{sum}(N)=0$. There are only eight cases,

  1. (i)

$4b_{1}+b_{16}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 1,16$;

  2. (ii)

$4b_{2}+b_{4}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 2,4$;

  3. (iii)

$4b_{3}+b_{7}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 3,7$;

  4. (iv)

$4b_{5}+b_{6}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 5,6$;

  5. (v)

$4b_{9}+b_{8}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 8,9$;

  6. (vi)

$4b_{10}+b_{12}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 10,12$;

  7. (vii)

$4b_{11}+b_{15}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 11,15$;

  8. (viii)

$4b_{13}+b_{14}=0$, and $b_{i}=0$ for $1\leq i\leq 17$ and $i\neq 13,14$.

It means that there are only eight product states in ${\cal H}_{{\cal
V}_{AB|CD}}^{\bot}$:
$(|10\rangle+|11\rangle+|20\rangle+|21\rangle-4|00\rangle)_{AB}(|01\rangle+|02\rangle)_{CD}$,
$(|12\rangle+|22\rangle)_{AB}(|00\rangle+|01\rangle+|10\rangle+|11\rangle-4|02\rangle)_{CD}$,
$(|11\rangle+|12\rangle+|21\rangle+|22\rangle-4|20\rangle)_{AB}(|12\rangle+|22\rangle)_{CD}$,
$(|20\rangle+|21\rangle)_{AB}(|10\rangle+|11\rangle+|20\rangle+|21\rangle-4|00\rangle)_{CD}$,
$(|01\rangle+|02\rangle+|11\rangle+|12\rangle-4|22\rangle)_{AB}(|20\rangle+|21\rangle)_{CD}$,
$(|00\rangle+|10\rangle)_{AB}(|11\rangle+|12\rangle+|21\rangle+|22\rangle-4|20\rangle)_{CD}$,
$(|00\rangle+|01\rangle+|10\rangle+|11\rangle-4|02\rangle)_{AB}(|00\rangle+|10\rangle)_{CD}$,
$(|01\rangle+|02\rangle)_{AB}(|01\rangle+|02\rangle+|11\rangle+|12\rangle-4|22\rangle)_{CD}$.
Since $\mathop{\rm Dim}({\cal H}_{{\cal V}_{AB|CD}}^{\bot})=16$, the OPS
${\cal V}_{AB|CD}$ is an SUCPB by Lemma 1.

Above all, ${\cal V}$ is an SUCPB in every bipartition.

The construction of UPB in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$
was generalized to any four-partite system
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\mathbb{C}^{d_{3}}\otimes\mathbb{C}^{d_{4}}$
for $d_{1},d_{2},d_{3},d_{4}\geq 3$ shi2021strong . Obviously, for the same
discussion as Lemma 3, we have the following theorem.

###### Theorem 2

In
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\mathbb{C}^{d_{3}}\otimes\mathbb{C}^{d_{4}}$,
$d_{1},d_{2},d_{3},d_{4}\geq 3$, there exists a UPB which is an SUCPB in every
bipartition.

However, not all UPBs have this property. For example, let

$\displaystyle|\psi_{1}\rangle$
$\displaystyle=|0\rangle_{A}(|0\rangle-|1\rangle)_{B},$
$\displaystyle|\psi_{2}\rangle$
$\displaystyle=(|0\rangle-|1\rangle)_{A}|2\rangle_{B},$
$\displaystyle|\psi_{3}\rangle$
$\displaystyle=|2\rangle_{A}(|1\rangle-|2\rangle)_{B},$
$\displaystyle|\psi_{4}\rangle$
$\displaystyle=(|1\rangle-|2\rangle)_{A}|0\rangle_{B},$
$\displaystyle|\psi_{5}\rangle$
$\displaystyle=(|0\rangle+|1\rangle+|2\rangle)_{A}(|0\rangle+|1\rangle+|2\rangle)_{B},$

then $\cup_{i=1}^{5}|\psi_{i}\rangle$ is a UPB in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}$ divincenzo2003unextendible . We can
construct a UPB in $\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$
from $\cup_{i=1}^{5}|\psi_{i}\rangle$ as follows,

$\displaystyle{\cal A}_{1}$
$\displaystyle:=\\{|\psi_{i}\rangle|0\rangle_{C}\mid 1\leq i\leq 5\\},$
$\displaystyle{\cal A}_{2}$
$\displaystyle:=\\{|i\rangle_{A}|j\rangle_{B}|1\rangle_{C}\mid
i,j\in\mathbb{Z}_{3}\\},$ $\displaystyle{\cal A}_{3}$
$\displaystyle:=\\{|i\rangle_{A}|j\rangle_{B}|2\rangle_{C}\mid
i,j\in\mathbb{Z}_{3}\\}.$

Let ${\cal W}:=\cup_{i=1}^{3}{\cal A}_{i}$. For any product state
$|\varphi\rangle=|\varphi_{1}\rangle_{A}|\varphi_{2}\rangle_{B}|\varphi_{3}\rangle_{C}\in{\cal
W}^{\bot}$, since $|\varphi\rangle$ is orthogonal to any state in ${\cal
A}_{2}\cup{\cal A}_{3}$, $|\varphi_{3}\rangle_{C}$ must be $|0\rangle_{C}$.
Further, since $\cup_{i=1}^{5}|\psi_{i}\rangle$ is a UPB, it means that
$|\varphi_{1}\rangle_{A}|\varphi_{2}\rangle_{B}$ cannot be a product state.
Thus ${\cal W}$ is a UPB in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$. Nevertheless,
${\cal W}$ is not a UCPB in every bipartition. There must exist four
orthogonal states
$\\{|\psi_{6}\rangle,|\psi_{7}\rangle,|\psi_{8}\rangle,|\psi_{9}\rangle\\}$
such that $\cup_{i=1}^{9}|\psi_{i}\rangle$ is an orthogonal basis in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}$. Then
$\\{\cup_{i=6}^{9}|\psi_{i}\rangle|0\rangle_{C}\\}\cup{\cal W}_{AB|C}$ is an
OPB in $AB|C$ bipartition. Thus, ${\cal W}$ is not a UCPB in $AB|C$
bipartition.

Next, we consider the application. Note that all UPBs in Theorems 1 and 2
cannot be perfectly distinguished under local POVMs and classical
communication in any bipartition. These UPBs can be used for local hiding of
information shi2022strong . For example, assume the information is encoded in
the UPB ${\cal U}$ given by Eq. (13), and the boss send it to his three
subordinates: A, B and C. These three subordinates are from different offices.
They can only perform local POVMs, and communicate classic information by
telephones. In this case, the three subordinates cannot obtain the full
information, even if any two of them are collusive. A and B are collusive
means that A and B are from the same office and can perform joint
measurements. See also Fig. 4. Further, Ref. shi2021strong  showed a stronger
property. Any UPB in Theorems 1 and 2 is locally irreducible111A set of
multipartite orthogonal states is locally irreducible if it is not possible to
eliminate one or more states from the set by orthogonality-preserving local
POVMs Halder2019Strong .  in every bipartition, which shows the phenomenon of
strong quantum nonlocality without entanglement Halder2019Strong .

Figure 4: The information is encoded to a tripartite UPB which is an SUCPB is
every bipartition, and the boss send it to his three subordinates: A, B, and
C. Even if any two of them are collusive, the three subordinates cannot obtain
the full information under local POVMs and classical communication.

##  IV UPBs in every bipartition

There exists another open question for UPBs demianowicz2018unextendible : can
we find a UPB, which is still a UPB in every bipartition? Such a UPB can be
used to construct genuinely entangled subspace, and it cannot be perfectly
distinguished under local POVMs and classical communication in any
bipartition. Unfortunately, any UPB in Sec. III is not a UPB in every
bipartition. We will give a sufficient condition for the existence of such a
UPB.

We can also generalize the tile structures to multipartite systems. A tile
structure ${\cal T}=\cup_{i=1}^{s}t_{i}$ in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\cdots\otimes\mathbb{C}^{d_{n}}$
is a $d_{1}\times d_{2}\times\cdots\times d_{n}$ hypercube, which can be
partitioned into $s$ disjoint tiles $\\{t_{i}\\}_{i=1}^{s}$. Each tile $t_{i}$
is a hypercube, and it can be expressed by

$\displaystyle t_{i}=$
$\displaystyle\\{x_{0}^{(1)},x_{1}^{(1)},\ldots,x_{k_{1}-1}^{(1)}\\}_{A_{1}}\times\\{x_{0}^{(2)},x_{1}^{(2)},\ldots,x_{k_{2}-1}^{(2)}\\}_{A_{2}}$
(17)
$\displaystyle\times\cdots\times\\{x_{0}^{(n)},x_{1}^{(n)},\ldots,x_{k_{n}-1}^{(n)}\\}_{A_{n}},$

where $\\{x_{0}^{(j)},x_{1}^{(j)},\ldots,x_{k_{j}-1}^{(j)}\\}$ is a subset of
$\mathbb{Z}_{d_{j}}$ for $1\leq j\leq n$. For tile $t_{i}$, we can construct
an OPS of size $k_{1}k_{2}\cdots k_{n}$ in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\cdots\otimes\mathbb{C}^{d_{n}}$,

$\displaystyle{\cal A}_{i}=\\{$
$\displaystyle\otimes_{j=1}^{n}\left(\sum_{e_{j}\in\mathbb{Z}_{k_{j}}}m_{a_{j},e_{j}}^{(j)}|x_{e_{j}}^{(j)}\rangle\right)_{A_{j}}$
(18) $\displaystyle\mid a_{j}\in\mathbb{Z}_{k_{j}},1\leq j\leq n\\},$

where each coefficient matrix
$M^{(j)}=(m_{a_{j},e_{j}}^{(j)})_{a_{j},e_{j}\in\mathbb{Z}_{k_{j}}}$ is a
$k_{j}\times k_{j}$ row orthogonal matrix, and $m_{0,e_{j}}^{(j)}=1$ for
$e_{j}\in\mathbb{Z}_{k_{j}}$. Then we obtain an OPB ${\cal
B}:=\cup_{i=1}^{s}{\cal A}_{i}$ in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\cdots\otimes\mathbb{C}^{d_{n}}$.
Note that

$|\psi_{i}\rangle=\otimes_{j=1}^{n}\left(\sum_{e_{j}\in\mathbb{Z}_{k_{j}}}|x_{e_{j}}^{(j)}\rangle\right)_{A_{j}}\in{\cal
A}_{i}.$ (19)

The “stopper” state is

$|S\rangle=\otimes_{j=1}^{n}\left(\sum_{r\in\mathbb{Z}_{d_{j}}}|r\rangle\right)_{A_{j}}.$
(20)

We may wonder whether the OPS

${\cal X}:=\cup_{i=1}({\cal
A}_{i}\setminus\\{|\psi_{i}\rangle\\})\cup\\{|S\rangle\\}$ (21)

is a UPB that is still a UPB in every bipartition.

For any bipartition $C\mid D$, where $C,D\subset\\{1,2,\ldots,n\\}$, and
$C\cup D=\\{1,2,\ldots,n\\}$, the OPB ${\cal B}$ must correspond to a tile
structure ${\cal T}_{C\mid D}$ in
$\mathbb{C}^{h_{1}}\otimes\mathbb{C}^{h_{2}}$, where $h_{1}=\prod_{g\in
C}d_{g}$ and $h_{2}=\prod_{g\in D}d_{g}$ (For example, see Fig. 3). By using
Lemma 2, if any $r$ ($2\leq r\leq s-1$) tiles in ${\cal T}_{C\mid D}$ cannot
form a rectangle, then ${\cal X}_{C\mid D}$ is a UPB in
$\mathbb{C}^{h_{1}}\otimes\mathbb{C}^{h_{2}}$. Note that in this case, ${\cal
X}$ is also a UPB in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\cdots\otimes\mathbb{C}^{d_{n}}$.
This is because if ${\cal X}$ is not a UPB, then there exists a product state
$|\psi\rangle$ in ${\cal X}^{\bot}$, and the product state
$|\psi\rangle_{C\mid D}$ in bipartition $C|D$ belongs to ${\cal X}_{C\mid
D}^{\bot}$, which contradicts ${\cal X}_{C\mid D}$ being a UPB in
$\mathbb{C}^{h_{1}}\otimes\mathbb{C}^{h_{2}}$. Now, we have the following
theorem.

###### Theorem 3

Consider a tile structure ${\cal T}=\cup_{i=1}^{s}t_{i}$ ($s\geq 5$) in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\cdots\otimes\mathbb{C}^{d_{n}}$.
For any bipartition $C\mid D$, if any $r$ ($2\leq r\leq s-1$) tiles in ${\cal
T}_{C\mid D}$ cannot form a rectangle, then the OPS ${\cal X}$ given by Eq.
(21) is a UPB in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\cdots\otimes\mathbb{C}^{d_{n}}$,
which is still a UPB in every bipartition.

One can use computer to search the tile structure in Theorem 3. By exhaustive
search, we show that such a tile structure does not exist in
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$. However, we
conjecture that such a tile structure may exist in a higher multipartite
system.

##  V conclusion and discussion

In this paper, we showed that there exist some unextendible product bases that
are uncompletable product bases in every bipartition in
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\mathbb{C}^{d_{3}}$ and
$\mathbb{C}^{d_{1}}\otimes\mathbb{C}^{d_{2}}\otimes\mathbb{C}^{d_{3}}\otimes\mathbb{C}^{d_{4}}$
for $d_{1},d_{2},d_{3},d_{4}\geq 3$, and
$\mathbb{C}^{3}\otimes\mathbb{C}^{3}\otimes\mathbb{C}^{3}$ achieved the
minimum system for the existence of such unextendible product bases. This
result answers an open question proposed in divincenzo2003unextendible . We
also showed that such unextendible product bases can be used for local hiding
of information. Finding an unextendible product basis that is still an
unextendible product basis in every bipartition can be challenging, and we
gave a sufficient condition for the existence of such an unextendible product
basis.

There are some interesting open questions left. How to find an unextendible
product basis that is still an unextendible product basis in every bipartition
by using Theorem 3? What is the minimum size of unextendible product basis
that is an uncompletable product basis in every bipartition?

## Acknowledgments

We thank Lin Chen for discussing this problem. The research of X.Z. and F.S.
were supported by the NSFC under Grants No. 12171452 and No. 11771419, the
Anhui Initiative in Quantum Information Technologies under Grant No.
AHY150200, the National Key Research and Development Program of China
(2020YFA0713100), and the Innovation Program for Quantum Science and
Technology (2021ZD0302904). M.-S.L. was supported by the NSFC under Grant No.
12005092, the China Postdoctoral Science Foundation (2020M681996).

## References and Notes

  * (1) C. H. Bennett, D. P. DiVincenzo, T. Mor, P. W. Shor, J. A. Smolin, and B. M. Terhal, “Unextendible product bases and bound entanglement,” _Physical Review Letters_ , vol. 82, no. 26, p. 5385, 1999.
  * (2) C. H. Bennett, D. P. DiVincenzo, C. A. Fuchs, T. Mor, E. Rains, P. W. Shor, J. A. Smolin, and W. K. Wootters, “Quantum nonlocality without entanglement,” _Physical Review A_ , vol. 59, no. 2, p. 1070, 1999.
  * (3) S. M. Cohen, “Understanding entanglement as resource: Locally distinguishing unextendible product bases,” _Physical Review A_ , vol. 77, no. 1, p. 012304, 2008.
  * (4) Z.-C. Zhang, X. Wu, and X. Zhang, “Locally distinguishing unextendible product bases by using entanglement efficiently,” _Physical Review A_ , vol. 101, no. 2, p. 022306, 2020.
  * (5) S. Halder, M. Banik, S. Agrawal, and S. Bandyopadhyay, “Strong quantum nonlocality without entanglement,” _Physical Review Letters_ , vol. 122, no. 4, p. 040403, 2019.
  * (6) F. Shi, M.-S. Li, L. Chen, and X. Zhang, “Strong quantum nonlocality for unextendible product bases in heterogeneous systems,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 55, no. 1, p. 015305, 2021.
  * (7) F. Shi, M.-S. Li, M. Hu, L. Chen, M.-H. Yung, Y.-L. Wang, and X. Zhang, “Strongly nonlocal unextendible product bases do exist,” _Quantum_ , vol. 6, p. 619, 2022.
  * (8) S. Bandyopadhyay, “More nonlocality with less purity,” _Physical Review Letters_ , vol. 106, no. 21, p. 210402, 2011.
  * (9) R. Augusiak, J. Stasińska, C. Hadley, J. Korbicz, M. Lewenstein, and A. Acin, “Bell inequalities with no quantum violation and unextendable product bases,” _Physical Review Letters_ , vol. 107, no. 7, p. 070401, 2011\.
  * (10) R. Augusiak, T. Fritz, M. Kotowski, M. Kotowski, M. Pawlowski, M. Lewenstein, and A. Acin, “Tight Bell inequalities with no quantum violation from qubit unextendible product bases,” _Physical Review A_ , vol. 85, no. 4, p. 042113, 2012.
  * (11) D. P. DiVincenzo, T. Mor, P. W. Shor, J. A. Smolin, and B. M. Terhal, “Unextendible product bases, uncompletable product bases and bound entanglement,” _Communications in Mathematical Physics_ , vol. 238, no. 3, pp. 379–410, 2003.
  * (12) M. Demianowicz and R. Augusiak, “From unextendible product bases to genuinely entangled subspaces,” _Physical Review A_ , vol. 98, no. 1, p. 012313, 2018\.
  * (13) M. Demianowicz, “Genuinely entangled subspaces of maximal dimensions (and many smaller ones) cannot be constructed from orthogonal unextendible product bases,” _arXiv preprint arXiv:2202.08356_ , 2022.
  * (14) F. Shi, X. Zhang, and L. Chen, “Unextendible product bases from tile structures and their local entanglement-assisted distinguishability,” _Physical Review A_ , vol. 101, no. 6, p. 062329, 2020.
  * (15) P. Horodecki, “Separability criterion and inseparable mixed states with positive partial transposition,” _Physics Letters A_ , vol. 232, no. 5, pp. 333–339, 1997.
  * (16) S. Halder, M. Banik, and S. Ghosh, “Family of bound entangled states on the boundary of the peres set,” _Physical Review A_ , vol. 99, no. 6, p. 062329, 2019.
  * (17) S. Agrawal, S. Halder, and M. Banik, “Genuinely entangled subspace with all-encompassing distillable entanglement across every bipartition,” _Physical Review A_ , vol. 99, no. 3, p. 032335, 2019.
  * (18) F. Shi, Z. Ye, L. Chen, and X. Zhang, “Strong quantum nonlocality in $N$-partite systems,” _Physical Review A_ , vol. 105, no. 2, p. 022209, 2022.