#  Generation of New Exciting Regressors for Consistent On-line Estimation of
Unknown Constant Parameters

Alexey Bobtsov  Senior Member, IEEE  Bowen Yi  Romeo Ortega  Life Fellow, IEEE
Alessandro Astolfi  Fellow, IEEE This work has been partially supported by the
European Union’s Horizon 2020 Research and Innovation Programme under grant
agreement No. 739551 (KIOS CoE), the Italian Ministry for Research in the
framework of the 2017 Program for Research Projects of National Interest
(PRIN), Grant no. 2017YKXYXJ, the Open Research Project of the State Key
Laboratory of Industrial Control Technology, Zhejiang University, China (No.
ICT2022B68), and the Russian Science Foundation, project no. 18-19-00627,
https://rscf.ru/project/18-19-00627. _(Corresponding author: Bowen Yi)_ A.
Bobtsov is with Hangzhou Dianzi University (HDU), 310018, Hangzhou, China and
the Faculty of Control Systems and Robotics, ITMO University, Kronverkskiy av.
49, St. Petersburg, 197101, Russia (bobtsov@mail.ru)B. Yi is with Australian
Centre for Field Robotics, The University of Sydney, Sydney, NSW 2006,
Australia. (bowen.yi@sydney.edu.au)R. Ortega is with Departamento Académico de
Sistemas Digitales, ITAM, Ciudad de México, México (romeo.ortega@itam.mx)A.
Astolfi is with the Department of Electrical and Electronic Engineering,
Imperial College London, London SW7 2AZ, UK and the DICII, Universita di Roma
“Tor Vergata”, Via del Politecnico 1, 00133 Roma, Italy<EMAIL_ADDRESS>
###### Abstract

The problem of estimating constant parameters from a standard vector linear
regression equation in the absence of sufficient excitation in the regressor
is addressed. The first step to solve the problem consists in transforming
this equation into a set of scalar ones using the well-known dynamic regressor
extension and mixing technique. Then a novel procedure to generate new scalar
exciting regressors is proposed. The superior performance of a classical
gradient estimator using this new regressor, instead of the original one, is
illustrated with comprehensive simulations.

parameter estimation, system identification, persistent excitation

##  1 Introduction

It is well known that consistent estimation of constant parameters from a
linear regression equation (LRE) with a gradient (or least squares) estimator
is possible only if the regressor satisfies certain excitation conditions. A
classical result shows that exponential convergence is possible if and only if
the regressor verifies the persistence of excitation (PE) condition [18]—which
is a uniform observability property for the associated linear time-varying
(LTV) system. It has recently been shown in [3, 16] that asymptotic (but not
exponential) convergence is guaranteed under the strictly weaker condition of
generalized PE, the definition of which may be found in [3, Proposition 6].

Unfortunately, the PE (or the generalized PE) property is rarely satisfied in
applications, hence the interest to propose new adaptation algorithms that
ensure parameter convergence without PE. This research line has been
intensively pursued in the last few years and some recent adaptive schemes, in
which the PE assumption is obviated via the incorporation of some off-line
data manipulation, have been reported in [5, 15, 4, 17]—see also [14] for a
recent survey.

In this paper we are interested in on-line estimation using recursive
algorithms. It is well known that, in contrast with off-line estimation
schemes, on-line estimation provides, via the accumulation of past
measurements and noise averaging, stronger robustness properties. Moreover, if
the adaptation gain of the estimator remains bounded away from zero—a property
usually referred as alertness111It is well known [18, Section 2.3.2] that, due
to the so-called covariance wind-up problem, the alertness property is lost in
standard least-squares estimators. Therefore, we concentrate on gradient-
descent schemes.—it has the ability of tracking slowly time-varying
parameters. We concentrate our attention to the case of a single uncertain
parameter. Our main motivation to study the scalar case stems from the recent
development of the dynamic regressor extension and mixing (DREM) estimator
[2], which is a procedure that generates, from a $q$-dimensional LRE, $q$
scalar LREs for each of the unknown parameters. It has been observed in
several applications that the absence of excitation stymies the successful use
of DREM. For instance, in [1] it is shown that consistent estimation of the
parameters of a linear time-invariant (LTI) system with DREM is possible if
and only if the original regressor is PE. Actually, since the key scalar
function that defines the convergence properties of the gradient estimator in
DREM is the determinant of the extended regressor, this converges in many
cases to zero, hence we only have excitation on a finite interval.

Our main contribution is to propose a procedure to generate, from a scalar
LRE, a new scalar LRE in which the new regressor satisfies the excitation
property of non-square integrability, even in the case in which the original
regressor is not sufficiently exciting—for instance an exponentially decaying
signal. It is shown in [2] that non-square integrability of the (scalar)
regressor is necessary and sufficient to ensure asymptotic convergence, which
becomes exponential imposing the PE condition. Instrumental for the
construction of the new LREs, that include some free signals, is to borrow the
key idea of the parameter estimation based observer (PEBO) proposed in [11],
later generalized in [12], which is a constructive procedure to design state
observers for state-affine nonlinear systems. Then, applying the energy
pumping-and-damping injection principle of [20], which was proposed as a
passivity-based orbital stabilization technique for port-Hamiltonian systems,
we select these signals to guarantee the desired excitation properties of the
new regressor.

The remainder of the paper is organized as follows. In Section 2 we briefly
recall the DREM procedure and review the problem of parameter estimation for a
scalar LRE. Section 3 presents the procedure to generate the new LRE with some
free signals, which are selected in Section 4 to comply with the excitation
injection requirement. Simulation results, which illustrate the superior
performance of the classical gradient estimator using the new regressor,
instead of the original one, are presented in Section 5. The paper is wrapped-
up with concluding remarks and a discussion on future research in Section 6.

Notation. $I_{n}$ is the $n\times n$ identity matrix. For a vector
$x\in\mathbb{R}^{n}$ we denote the Euclidean norm as $|x|:=\sqrt{x^{\top}x}$.
${\cal L}_{1}$ and ${\cal L}_{2}$ denote the absolute integrable and square
integrable function spaces, respectively, and ${\cal L}_{\infty}$ represents
the vector space of essentially bounded functions.

##  2 Gradient Estimation of a Single Parameter

In this paper we deal with the problem of on-line estimation of the unknown
constant parameters $\theta\in\mathbb{R}^{q}$ appearing in an LRE of the form

$w=\psi^{\top}\theta,$ (1)

where $w(t)\in\mathbb{R}$ and $\psi(t)\in\mathbb{R}^{q}$ are measurable
signals. This problem appears in several applications including system
identification [9] and adaptive control [7, 10, 18] in which, as discussed in
Section 1, a key requirement to achieve their objectives is that the regressor
$\psi$ is PE, a condition that is rarely satisfied in practice. Our task is
then to generate new LREs that satisfy the PE requirement.

###  2.1 Generation of scalar LRE

The first step in our design is to apply the DREM procedure [2] to obtain $q$
scalar LREs, one for each of the unknown parameters. Towards this end, we
introduce a linear, single-input $q$-output, bounded-input bounded-output
(BIBO)–stable operator ${\cal H}$ and define the vector $W\in\mathbb{R}^{m}$
and the matrix $\Psi\in\mathbb{R}^{q\times q}$ as

$W:={\cal H}[w],\;\Psi:={\cal H}[\psi^{\top}].$

Clearly, because of linearity and BIBO stability, these signals satisfy

$W=\Psi\theta.$

At this point the key step of regressor “mixing” of the DREM procedure is used
to obtain a set of $q$ scalar equations multiplying from the left the vector
equation above by the adjunct matrix, denoted $\mbox{adj}\\{\cdot\\}$, of
$\Psi$ to get

$y_{i}=\Delta\theta_{i},\;i\in\bar{q}:=\\{1,2,\dots,q\\},$ (2)

in which we have defined

This fundamental modification has numerous advantages and DREM has been
instrumental to solve many, previously open, problems—see [14, 13, 12] for a
recent account of some of these results.

###  2.2 Properties of the gradient estimator

Motivated by the developments above, in the remaining part of the paper we
consider scalar LREs of the form (2). The following property of the gradient
estimator is easy to establish [2].

###### Proposition 1

Consider the scalar LREs (2) with222To simplify the notation we omit the
subindex $(\cdot)_{i}$. $y(t)\in\mathbb{R}$ and $\Delta(t)\in\mathbb{R}$
bounded, measurable, signals and $\theta\in\mathbb{R}$ an unknown parameter
estimated on-line via the gradient descent adaptation algorithm

$\dot{\hat{\theta}}=-\gamma\Delta(\Delta\hat{\theta}-y),$ (3)

with $\gamma>0$ the adaptation gain.

  1. P1

The following equivalence holds true:

$\lim_{t\to\infty}|\hat{\theta}(t)-\theta|=0~{}~{}\iff~{}~{}\Delta\notin{\cal
L}_{2}.$

  2. P2

The convergence of the estimate $\hat{\theta}$ to $\theta$ is exponential if
and only if $\Delta$ is PE, that is, if there exist $T>0$ and $\delta>0$ such
that

$\int_{t}^{t+T}\Delta^{2}(\tau)d\tau\geq\delta,\;\forall t\geq 0.$

##  3 Generation of the New LRE

In this section, in order to generate a new LRE, we apply some filters to both
sides of the original LRE (2) with some free terms to comply with the
excitation requirement. Then, in Section 4 we study how to design these terms
to fulfill this task.

In the sequel we apply the construction used in generalized parameter
estimation based observer (GPEBO) [12] to create a new LRE from the scalar LRE
(2).

###### Proposition 2

Consider the scalar LRE (2). There exists a measurable signal ${\cal
Y}_{2}(t)\in\mathbb{R}$ such that the new LRE

${\cal Y}_{2}=\Phi_{21}\theta$ (4)

holds, in which the new regressor $\Phi_{21}(t)\in\mathbb{R}$ is obtained from
the solution of the ordinary differential equation

$\begin{bmatrix}\dot{\Phi}_{11}\\\
\dot{\Phi}_{21}\end{bmatrix}=\begin{bmatrix}0&u_{1}\\\
u_{2}\Delta&u_{3}\end{bmatrix}\begin{bmatrix}\Phi_{11}\\\
\Phi_{21}\end{bmatrix},$ (5)

with initial conditions

$\begin{bmatrix}\Phi_{11}(0)\\\ \Phi_{21}(0)\end{bmatrix}=\begin{bmatrix}1\\\
0\end{bmatrix},$ (6)

and arbitrary bounded signals $u_{i}(t)\in\mathbb{R},$ $i=1,2,3$.

######  Proof 3.1.

Define the scalar dynamics

$\dot{z}=u_{2}y+u_{3}z,\;z(0)=0,$ (7)

and note that, since $\theta$ is constant, we can write

$\dot{\theta}=u_{1}(z-z).$ (8)

Combining (7) and (8), and using (2), we can write the “virtual” LTV system

$\dot{x}=A(t)x+b(t)$ (9)

with

$\displaystyle x$ $\displaystyle:=\begin{bmatrix}\theta\\\ z\end{bmatrix}$
$\displaystyle A(t)$ $\displaystyle:=\begin{bmatrix}0&u_{1}(t)\\\
u_{2}(t)\Delta(t)&u_{3}(t)\end{bmatrix},$ $\displaystyle b(t)$
$\displaystyle:=\begin{bmatrix}-u_{1}(t)x_{2}(t)\\\ 0\end{bmatrix},$

and initial conditions

$x(0)=\begin{bmatrix}\theta\\\ 0\end{bmatrix}.$ (10)

Following the GPEBO approach define the dynamics

$\quad\quad\quad~{}\dot{\xi}=A(t)\xi+b(t),\;\xi(0)={\bf 0}_{2\times 1}$ (11a)
$\dot{\Phi}=A(t)\Phi,\;\Phi(0)=I_{2}.$ (11b)

Clearly, the state transition matrix $\Omega(t,s)$ of $A(t)$ is given by
$\Omega(t,s)=\Phi(t)\Phi^{-1}(s)$.

The error signal

$e:=\xi-x$ (12)

satisfies $\dot{e}=A(t)e$. Consequently, from (12) and the properties of the
system matrix $\Phi$ defined in (11b), we get

$\displaystyle x$ $\displaystyle=\xi-\Phi e(0)$
$\displaystyle=\xi+\Phi\begin{bmatrix}\theta\\\ 0\end{bmatrix}$
$\displaystyle=\xi+\begin{bmatrix}\Phi_{11}\\\ \Phi_{21}\end{bmatrix}\theta,$

where, to get the second identity, we have taken into account (10) and the
initial conditions in (11a), and introduced the notation

$\Phi=\begin{bmatrix}\Phi_{11}&\Phi_{12}\\\ \Phi_{21}&\Phi_{22}\end{bmatrix}.$

Note now that

$\begin{bmatrix}y\\\ z\end{bmatrix}=\begin{bmatrix}\Delta&0\\\
0&1\end{bmatrix}x=\begin{bmatrix}\Delta&0\\\
0&1\end{bmatrix}\Big{(}\xi+\begin{bmatrix}\Phi_{11}\\\
\Phi_{21}\end{bmatrix}\theta\Big{)}.$

The proof is then completed defining

${\cal Y}=\begin{bmatrix}{\cal Y}_{1}\\\ {\cal
Y}_{2}\end{bmatrix}:=\begin{bmatrix}y\\\
z\end{bmatrix}-\begin{bmatrix}\Delta&0\\\ 0&1\end{bmatrix}\xi,$ (13)

noting that

${\cal Y}_{2}=z-\xi_{2}$

is a measurable signal, and computing the dynamics of the first column of
$\Phi$.

The procedure to generate the new LRE described in Proposition 2 is summarized
in the diagram in Fig. 1.

Figure 1: Block diagram of the proposed LRE generator

##  4 Improving the Excitation of the Regressor

###  4.1 Main Result

In this section we follow the basic idea of the energy pumping-and-damping
injection construction used in [20] for orbital stabilization to select the
signals $u_{i}$ in (5) to provide excitation to the new regressor $\Phi_{21}$.

To articulate the main result we need the following assumption.

###### Assumption 1

The bounded scalar signal $\Delta$ in (2) is interval exciting (IE), [8, 19],
that is there exist a time instant $t_{c}>0$ and a constant $\delta>0$ such
that

$\int_{0}^{{t_{c}}}\Delta^{2}(\tau)d\tau\geq\delta.$

######  Proposition 4.2.

Assume $\Delta$ verifies Assumption 1. Define the signals

$\displaystyle u_{1}$ $\displaystyle=-\alpha\Delta$ (14) $\displaystyle u_{2}$
$\displaystyle=\alpha$ $\displaystyle u_{3}$
$\displaystyle=-\tilde{V}(\Phi_{11},\Phi_{21}),$

with $\alpha(t)\in\mathbb{R}$ a bounded signal such that
$\alpha(t)\Delta(t)\not\equiv 0$ for some $t$,

$\alpha\Delta\in{\cal L}_{1},$ (15)

and

$\tilde{V}(\Phi_{11},\Phi_{21}):={1\over
2}(\Phi_{11}^{2}+\Phi_{21}^{2})-\beta,$

where $\Phi_{11},\Phi_{21}$ are the solutions of (5) and $0<\beta<{1\over 2}$.
The resulting dynamics verifies the following properties.

  1. F1

Either $\Phi_{21}\not\in{\cal L}_{2}$ or

$\Phi^{2}_{11}(t)+\Phi^{2}_{21}(t)\geq 2\beta+\varepsilon,\quad\forall t\geq
0,$ (16)

for some (sufficiently small) $\varepsilon>0$.

  2. F2

The full state of the system—$\Phi,\xi$ and $z$—is bounded.

######  Proof 4.3.

Replacing the signals $u_{i}$ in (14) to (5) we get

$\displaystyle\dot{\Phi}$ $\displaystyle={A(t)}\Phi$
$\displaystyle=\begin{bmatrix}0&-{\alpha}\Delta\\\
{\alpha}\Delta&-\tilde{V}\end{bmatrix}\Phi.$

The dynamics of the remaining states of the regressor generator, that is,
$\xi$ and $z$, is given by

$\begin{bmatrix}\dot{\xi}_{1}\\\ \dot{\xi}_{2}\\\
\dot{z}\end{bmatrix}=\begin{bmatrix}0&-\alpha\Delta&\alpha\Delta\\\
\alpha\Delta&-\tilde{V}&0\\\
0&0&-\tilde{V}\end{bmatrix}\begin{bmatrix}\xi_{1}\\\ \xi_{2}\\\
z\end{bmatrix}+\begin{bmatrix}0\\\ 0\\\ \alpha y\end{bmatrix}.$

First, note that the IE Assumption 1 rules out the extreme case
$y(t)=\Delta(t)\equiv 0$. Indeed, in this case all signals of the regressor
generator remain at their initial conditions, which is an equilibrium. To
avoid this situation we also require the technical assumption that
$\alpha(t)\Delta(t)\not\equiv 0$ for some $t\geq 0$.

From the first column of the dynamic equation of $\Phi$ we immediately get

$\dot{\tilde{V}}=-\Phi^{2}_{21}\tilde{V}\leq 0,$ (17)

from which we conclude boundedness of $(\Phi_{11},\Phi_{21})$ as well as the
invariance of the set

$\Omega:=\\{q\in\mathbb{R}^{2}|\tilde{V}(q)=0\\}.$

Now, invoking the initial condition constraint (6) yields

$\displaystyle\quad{1\over 2}(\Phi^{2}_{11}(0)+\Phi^{2}_{21}(0))={1\over 2}$
$\displaystyle\iff$
$\displaystyle\quad\tilde{V}(\Phi_{11}(0),\Phi_{21}(0))+\beta={1\over 2}$
$\displaystyle\iff$
$\displaystyle\quad\tilde{V}(\Phi_{11}(0),\Phi_{21}(0))={1\over 2}-\beta$
$\displaystyle\implies$
$\displaystyle\quad\tilde{V}(\Phi_{11}(0),\Phi_{21}(0))>0,$

where we have used the fact that $\beta<{1\over 2}$ to get the last
implication. The latter inequality implies that the trajectory starts outside
the disk described by the set $\Omega$. This, together with the invariance of
the set implies that the whole trajectory $(\Phi_{11}(t),\Phi_{21}(t))$ is
outside this disk, that is,

$\tilde{V}(t):=\tilde{V}(\Phi_{11}(t),\Phi_{21}(t))\geq 0,\;\forall t\geq 0.$
(18)

Let us, now, consider the two possible, mutually exclusive, cases:

  1. (i)

$\lim_{t\to\infty}\tilde{V}(t)=0$

  2. (ii)

$\lim_{t\to\infty}\tilde{V}(t)\neq 0$.

Consider, first, the second case and note that, due to the invariance of the
set $\Omega$, this case also rules out the possibility that the convergence to
zero happens in finite time. In this case the inequality (16) holds true.

In the first case, integrating (17) we see that the following equivalence
holds true:

$\lim_{t\to\infty}\tilde{V}(t)=0\quad\iff\quad\Phi_{21}\not\in{\cal L}_{2},$

completing the proof of the claim F1.

Now we show boundedness of the remaining states. For the second column of the
matrix $\Phi$, note that it satisfies the dynamics

$\displaystyle\begin{bmatrix}\dot{\Phi}_{12}\\\ \dot{\Phi}_{22}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}0&-{\alpha}\Delta\\\
{\alpha}\Delta&-\tilde{V}\end{bmatrix}\begin{bmatrix}\Phi_{12}\\\
\Phi_{22}\end{bmatrix},\;\begin{bmatrix}\Phi_{12}(0)\\\
\Phi_{22}(0)\end{bmatrix}=\begin{bmatrix}0\\\ 1\end{bmatrix}.$

Consider the function

$W(\Phi_{12},\Phi_{22}):={1\over 2}(\Phi_{12}^{2}+\Phi_{22}^{2}),$

the derivative of which along the trajectories of the system (5) satisfies

$\dot{W}=-\Phi^{2}_{22}\tilde{V}\leq 0,$

where the bound—from which we conclude the boundedness of
$(\Phi_{12},\Phi_{22})$—follows from (18).

Finally, from (12) we have

$\displaystyle\xi$ $\displaystyle=x+e$ $\displaystyle=\begin{bmatrix}\theta\\\
z\end{bmatrix}+\Phi\begin{bmatrix}\theta\\\ 0\end{bmatrix}.$

Hence, if $z$ is bounded we conclude that $\xi$ is also bounded completing the
proof.

To prove boundedness of $z$ we, again, consider the two scenarios (i) and (ii)
described above. For the second case we consider the function

$V_{z}(z):={1\over 2}z^{2}.$

Its time derivative along the trajectories of (5) satisfies

$\displaystyle\dot{V}_{z}$ $\displaystyle=-\tilde{V}z^{2}+\alpha yz$
$\displaystyle\leq-(\tilde{V}-{\epsilon\over 2})z^{2}+{1\over
2\epsilon}(\alpha y)^{2}$
$\displaystyle\leq-(\varepsilon-\epsilon)V_{z}+{1\over 2\epsilon}(\alpha
y)^{2}$ $\displaystyle=-{\varepsilon\over 2}V_{z}+{1\over 4\varepsilon}(\alpha
y)^{2},$

where we have used the fact that, for all $\epsilon>0$, we have that

$\alpha yz\leq{1\over 2\epsilon}(\alpha y)^{2}+{\epsilon\over 2}z^{2},$

to get the first inequality, the fact that $\tilde{V}\geq{\varepsilon\over
2}$, which follows from (16), to get the second bound and selected
$\epsilon={\varepsilon\over 2}$ in the final identity. The proof of
boundedness of $z$ follows, then, from the inequality above and the fact that
$\alpha y$ is bounded.

Now, consider the scenario (i). It is clear that the solution of the equation

$\dot{z}=-\tilde{V}z+\alpha\Delta\theta,$ (19)

is given by

$z(t)=\left(\int_{0}^{t}e^{-\int_{\sigma}^{t}\tilde{V}(s)ds}\alpha(\sigma)\Delta(\sigma)d\sigma\right)\theta.$
(20)

From (20) we immediately obtain

$\displaystyle|z(t)|$
$\displaystyle\leq\bigg{|}\int_{0}^{t}e^{-\int_{\sigma}^{t}\tilde{V}(s)ds}\alpha(\sigma)\Delta(\sigma)d\sigma\bigg{|}|\theta|$
$\displaystyle\leq\int_{0}^{t}\Big{|}e^{-\int_{\sigma}^{t}\tilde{V}(s)ds}\Big{|}|\alpha(\sigma)\Delta(\sigma)|d\sigma|\theta|$
$\displaystyle\leq\int_{0}^{t}|\alpha(\sigma)\Delta(\sigma)|d\sigma|\theta|,$

where, to obtain the last bound, we have used the inequality (18), which
implies that

$\left|\exp\left(-\int_{\sigma}^{t}\tilde{V}(s)ds\right)\right|\leq 1.$ (21)

Boundedness of $z$ is concluded, from the inequality above, if
$\alpha\Delta\in{\cal L}_{1}$.

Clearly, invoking the equivalence in the claim P1 of Proposition 1, we have
that the condition $\Phi_{21}\not\in{\cal L}_{2}$ of the claim F1 of
Proposition 2 ensures global convergence of a gradient estimator with the new
LRE. On the other hand, the inequality (16) in the claim F1 guarantees the
following excitation properties for the new regressor $\Phi_{21}$. If we can
ensure that $\Phi_{11}\not\to\sqrt{2\beta}$ then $\Phi_{21}\not\to 0$ and,
consequently, $\Phi_{21}$ is PE.333Recall that if a scalar signal is PE then
it is not square integrable, but the converse is not true [2]. In this sense,
the new regressor $\Phi_{21}$ is “more exciting” than the original IE
regressor $\Delta$. Although we have not been able to prove this property, it
has systematically been true in all our simulations, some of which are given
in Section 5. Besides, in the following corollary we identify a scenario which
guarantees that the new regressor $\Phi_{21}\notin{\cal L}_{2}$.

######  Corollary 4.4.

Consider the filter design in Proposition 4.2. If
$\lim_{t\to\infty}\alpha(t)\Delta(t)=0$, then $\Phi_{21}\notin{\cal
L}_{2}$.444It equivalently means that the second case in (16) does not occur.

######  Proof 4.5.

We prove this fact by contradiction assuming that $\Phi_{21}\in{\cal L}_{2}$.
According to the analysis in Proposition 4.2, $\tilde{V}$ has a limit when
$t\to+\infty$, which can be zero or a positive constant, invoking non-
negativeness of $\tilde{V}$ from (18). Clearly, for the case
$\Phi_{21}\in{\cal L}_{2}$, we have

$\lim_{t\to\infty}\tilde{V}(t)=\varepsilon$

for some $\varepsilon>0$. From $\Phi_{21}\in{\cal L}_{2}$ and
$\dot{\Phi}_{21}\in{\cal L}_{\infty}$, we have $\Phi_{21}\to 0$, and
consequently

$\lim_{t\to\infty}\Phi_{11}(t)=\sqrt{2(\beta+\varepsilon)}.$

Rewrite now $(\Phi_{11},\Phi_{21})$ in polar coordinates, that is

$\Phi_{11}=r\cos\rho,\quad\Phi_{21}=r\sin\rho,$

with

$\displaystyle\rho$
$\displaystyle:=\arctan\left({\Phi_{21}\over\Phi_{11}}\right),$ (22)
$\displaystyle r$ $\displaystyle:=\sqrt{\Phi_{11}^{2}+\Phi_{21}^{2}}.$

The transformation is well-posed and bijective almost globally except of the
origin, which the above analysis proves that is not possible. In terms of the
convergence of $\Phi_{21}$ and $\Phi_{22}$, we have

$\lim_{t\to\infty}\rho(t)=\rho_{\star},\quad\lim_{t\to\infty}r(t)=\sqrt{2(\beta+\varepsilon)},$
(23)

for some constant $\rho_{\star}\in[0,2\pi)$. On the other hand, the dynamics
of $\rho$ is given by

$\displaystyle\dot{\rho}$ $\displaystyle={1\over
1+({\Phi_{21}\over\Phi_{11}})^{2}}\cdot{\dot{\Phi}_{21}\Phi_{11}-\dot{\Phi}_{11}\Phi_{22}\over\Phi_{11}^{2}}$
$\displaystyle=\alpha\Delta-{1\over 2}+{\beta\over r^{2}}.$

Invoking the convergence condition $\rho\to\rho_{\star}$ as $t\to\infty$, the
boundedness of $\dot{\rho}$ and the $C^{1}$-smoothness of the solution
$\rho(t)$, we have

$\lim_{t\to\infty}\dot{\rho}(t)=\lim_{t\to\infty}\left[\alpha(t)\Delta(t)-{1\over
2}+{\beta\over r^{2}(t)}\right]=0.$

From the convergence condition $\alpha\Delta\to 0$ and
$r\to\sqrt{2(\beta+\varepsilon)}$, we have $\varepsilon=0$, which contradicts
the assumption $\varepsilon>0$. Therefore, $\Phi_{21}\notin{\cal L}_{2}$,
completing the proof.

Two interesting observations are, first, that the condition
$\alpha\Delta\in{\cal L}_{1}$ cannot be replaced by $\alpha\Delta\in{\cal
L}_{2}$. Indeed, a counterexample to this claim is given by the selection

$\tilde{V}(t)=e^{-t},\;\alpha(t)\Delta(t)={1\over t+1},$

for which $\alpha\Delta\in{\cal L}_{2}$, but $\alpha\Delta\not\in{\cal
L}_{1}$. In this case the system (19) has unbounded solutions. Second, the
condition $\alpha\Delta\in{\cal L}_{1}$ is sufficient, but not necessary with
a counterexample to this claim given by the selection

$\tilde{V}(t)=e^{-t},\;\alpha(t)\Delta(t)={\sin t\over t+1},$

that does not satisfy the ${\cal L}_{1}$ condition but with associated
solutions (19) bounded. Note also that ${\sin t\over t+1}\in{\cal L}_{2}$!

###  4.2 Robustness

We now consider robustness of the proposed scheme.555Robustness of the filter
in Proposition 2 is studied in the BIBO sense from the perturbed term to the
filtered output, which is a property of the filter regardless of collected
data. For simplicity we gather all perturbations in a term $\delta(t)$, i.e.,
the original LRE becomes

$y_{\tt N}(t)=\Delta(t)\theta+\delta(t),$

in which $(\cdot)_{\tt N}$ is used to represent the perturbed signal of the
signal $(\cdot)$. Assuming that the selection of $\alpha$ is independent of
the regression output $y$, the filter dynamics (7) becomes

$\dot{z}_{\tt N}=u_{2}y_{\tt N}+u_{3}z_{\tt N},~{}z_{\tt N}(0)=0.$ (24)

From linearity we have

$z_{\tt N}=z+\delta_{z},$

with $z$ the solution of $\dot{z}=u_{2}y+u_{3}z$ from $z(0)=0$ and

$\delta_{z}(t)=\int_{0}^{t}\exp\left(\int_{s}^{t}-\tilde{V}(\tau)d\tau\right)\alpha(s)\delta(s)ds.$
(25)

Due to the above parameterization of the perturbation, $\delta$ does not
affect $A(t)$, equivalently, the solution of $\Phi(t)$. On the other hand, the
$\xi$-dynamics is cascaded to the filter (7), i.e., for the perturbed case

$\dot{\xi}_{\tt N}=A(t)\xi_{\tt N}+\begin{bmatrix}\alpha\Delta z_{\tt N}\\\
0\end{bmatrix}=A(t)\xi_{\tt N}+\begin{bmatrix}\alpha\Delta(z+\delta_{z})\\\
0\end{bmatrix},$

with $\xi_{\tt N}(0)=\mathbf{0}_{2\times 1}$. From superposition we have
$\xi_{\tt N}=\xi+\delta_{\xi}$ with

$\delta_{\xi}(t):=\int_{0}^{t}\Phi(t)\Phi^{-1}(s)\alpha(s)\Delta(s)\begin{bmatrix}\delta_{z}(s)\\\
0\end{bmatrix}ds.$

In terms of ${\cal Y}_{2}=z-\xi_{2}$, we have the perturbed new LRE as

${\cal Y}_{2}=\Phi_{21}\theta+\delta_{\tt N},$ (26)

with $\delta_{\tt N}:=\delta_{z}-\delta_{\xi,2}$. Some observations are in
order.

  * -

The perturbation $\delta$ from the original LRE (25) affects the new LRE (26)
in a different way.

  * -

Consider the perturbation $\delta$ as a zero-average high-frequency
measurement noise. Noting that $\tilde{V}(t)\geq 0$ for all $t\geq 0$, roughly
speaking, the filter (7) plays a role similarly to a low-pass filter (which
may be analyzed via standard averaging analysis), making the new LRE robust to
the high-frequency component in $\delta$.

  * -

For biased but bounded perturbations, e.g., environmental disturbances and
unmodelled dynamics, we have

$\displaystyle|\delta_{z}|$
$\displaystyle\leq\int_{0}^{t}\left|\exp\left(\int_{s}^{t}-\tilde{V}(\tau)d\tau\right)\right||\alpha(s)||\delta(s)|ds$
$\displaystyle\leq\int_{0}^{t}|\alpha(s)|\|\delta\|_{\infty}ds,$

in which we have used (21). Clearly, selecting $\alpha\in{\cal L}_{1}$ can
guarantee $\delta_{z}\in{\cal L}_{\infty}$, and a vanishing signal $\alpha$
may also yield $\delta_{\xi}\in{\cal L}_{\infty}$, thus providing a BIBO
stability property.

##  5 Simulations

In this section we provide comprehensive simulation results to verify our main
claims. In all simulations the parameter for pumping-and-damping injection is
selected as $\beta=0.4$ and—following the construction of the filter proposed
above—the initial conditions are fixed at

$z(0)=0,\;\Phi(0)=I_{2},\;\xi(0)=0.$

We first consider a constant unknown parameter $\theta=-5$, and test the
proposed regressor generator under different excitation conditions for the
original regressor $\Delta$. Namely, they are chosen as follows:

$\displaystyle\Delta_{1}(t)$ $\displaystyle={1\over 2}e^{-t}$
$\displaystyle\Delta_{2}(t)$ $\displaystyle=\left\\{\begin{aligned} {1\over
2}&&t\in[0,5]\mbox{s}\\\ 0&&t>5\mbox{s}.\end{aligned}\right.$

Note that these signals are clearly not PE and they belong to ${\cal L}_{2}$,
hence the gradient estimator (3) does not ensure parameter convergence. On the
other hand, they satisfy the extremely weak condition that they are IE. For
these regressor we use $\alpha(t)\equiv 1$, which verifies
$\alpha\Delta\in{\cal L}_{1}$, as well as the additional condition in
Corollary 4.4.

To estimate the parameter we use the standard gradient descent adaptation with
the old and the new regressors, that is,

$\dot{\hat{\theta}}_{\tt old}=\gamma\Delta(y-\Delta\hat{\theta}_{\tt old}),$

and

$\dot{\hat{\theta}}_{\tt
new}=\gamma\Phi_{21}(\mathcal{Y}_{2}-\Phi_{21}\hat{\theta}_{\tt new}),$

selecting, in both cases, $\gamma=2$. The initial conditions are selected as
$\hat{\theta}_{\tt new}(0)=0$ and $\hat{\theta}_{\tt old}(0)=0.$ As we can see
from Figs. 2-3, the proposed regressor generator transforms the interval
exciting regressors into PE regressors ensuring that the estimate
$\hat{\theta}_{\tt new}$ exponentially converges to the true value. On the
other hand, the estimates generated with the original regressors
$\hat{\theta}_{\tt old}$ exhibit a steady state error. In these figures, we
also plot the evolution of the Frobenius norm of $\Phi(t)$, that is
$\|\Phi\|_{F}=\sqrt{{\mbox{trace}}(\Phi^{\top}\Phi)}$. It is also observed
that all the states are bounded.

(a) Original and new regressors

(b) Estimates using the new and original LRE

(c) Internal states

(d) The value of $\int_{0}^{t}|\alpha(s)\Delta(s)|$ds

Figure 2: Simulation results for the signal $\Delta_{1}(t)={1\over 2}\exp(-t)$

(a) Original and new regressors

(b) Estimates using the new and original LREs

(c) Internal states

(d) The value of $\int_{0}^{t}|\alpha(s)\Delta(s)|$ds

Figure 3: Simulation results for the step signal $\Delta_{2}(t)$

Now, let us consider another vanishing signal

$\Delta_{3}(t)={1\over 2t+1}.$

It is clear that a non-zero constant $\alpha(t)$ cannot guarantee the
condition $\alpha\Delta\in{\cal L}_{1}$. Therefore, we parameterize $\alpha$
as

$\alpha(t)=\alpha_{0}(t)-ky(t)z(t),$

with $k>0$. Then, the key equation (7) takes the form

$\dot{z}=-(\tilde{V}+ky^{2})z+\alpha_{0}y,$

with the second term in $\alpha(t)$ providing additional damping in the
dynamics of $z$. By carefully selecting $\alpha_{0}$, we may guarantee that
$\alpha\Delta\in{\cal L}_{1}$. For $\Delta_{3}$ we select $k=0.1$ and
$\alpha_{0}(t)=e^{-{t\over 10}}$ in simulation, see Fig. 5 for the
corresponding results. It also illustrates the theoretical analysis.

We then add high-frequency noise in the LRE output $y$ via the “Uniform Random
Number” block in Simulink/Matlab${}^{\text{TM}}$. The signal-to-noise ratio
(SNR) is selected around 20 dB with sampling time 0.01s. In this example we
consider the signal $\Delta_{2}$. The simulation results are given in Fig. 6,
from which we observe that the new regressor output ${\cal Y}_{2}$ is hardly
affected by the noise, as analyzed in Subsection 4.2.

(a) Regression outputs $y(t)$ and ${\cal Y}_{2}$

(b) Estimate using the new LRE

(c) Internal states (d) The value of $\int_{0}^{t}|\alpha(s)\Delta(s)|$ds
Figure 4: Simulation results for the signal $\Delta_{3}(t)={1\over 2t+1}$

(a) Outputs in the original and new regressors

(b) Estimates using the new LRE

Figure 5: Simulation results in the presence of high-frequency measurement
noise for the signal $\Delta_{2}(t)$

As is well known, one of the main motivations of on-line identification is to
track (slowly) time-varying parameters. To verify the alertness to the
variation of the unknown parameter of the proposed procedure a simulation in
which the parameter $\theta$ jumps from $-5$ to $-4$ at $t=10$s has been
carried out with the IE regressor

$\Delta_{4}(t)=\left\\{\begin{aligned} \sin\left({\pi\over
10}t\right),&&t\in[0,12]\mbox{s}\\\ 0,&&t>12\mbox{s}\end{aligned}\right.$

and $\gamma=2$ and $\alpha(t)\equiv 1$.666We need to consider the signal being
IE after the parameter variation for a while. Otherwise, the change of
parameter cannot be captured by the output $y$. The simulation results are
shown in Fig. 6 from which we empirically conclude that the proposed new
regressor is capable of dealing with parameter variations. The estimation
$\hat{\theta}_{\tt old}$ from the original regressor is also presented, which,
as we can see, fails to get satisfactory performance due to the lack of PE.

We underscore that there is an ultimate error from $\hat{\theta}_{\tt new}$,
because Proposition 2 only holds when $\theta$ is constant. For a time-varying
$\theta$, a perturbation term appears in the new LRE (4), inducing the
observed estimation error. As future work, it would be interesting to study
the ultimate error quantitatively. To overcome this issue, re-initialization
of the filter in Proposition 2 after parameter variation may be implemented,
and it would be interesting to study how to re-initialize it automatically.

Figure 6: Alertness to time-varying parameter $\theta$ for the IE signal
$\Delta_{4}(t)$

##  6 Conclusions and Future Research

We have presented a procedure to generate a new LRE in which the new regressor
has some excitation properties, even if the original regressor does not. The
key idea to carry out this task is the use of the virtual dynamics
$\dot{\theta}=0$ given in (8). As seen in (5), the new LRE has some free
signals $u_{i}$ that we selected in a particular way in Proposition 4.2 to
enforce the excitation-like condition of F1. However, further research is
needed to select other signals $u_{i}$ that would ensure bona fide excitation
for a well-defined class of regressors $\Delta$—for instance, IE. It would be
particularly interesting to relax the standing assumption (15)—that imposes
constraints on the tuning function $\alpha$.

A very simple extension of the result is obtained combining the elements of
the signal ${\cal Y}$ in (13) as

$Y:={\cal H}_{1}({\cal Y}_{1})+{\cal H}_{2}({\cal Y}_{2}),$

with ${\cal H}_{1}$ and ${\cal H}_{2}$ BIBO stable, scalar, LTV operators, to
generate a new LRE

$Y=[{\cal H}_{1}(\Delta\Phi_{11})+{\cal H}_{2}(\Phi_{21})]\theta.$

A topic of future research is the selection of suitable operators ${\cal
H}_{i}$ to ensure excitation of the new regressor and boundedness of all
signals. Note also that the extension of the results—including the application
of the DREM procedure—to the case of nonlinear, separable parameterizations is
straightforward.

Our LRE generator is restricted to the case of scalar LREs. However, mimicking
the derivations of Section 3 it is possible to generate a new LRE also for the
vector case, that is, when the original LRE is of the form (1). Indeed, using
GPEBO it is possible to construct a measurable signal $e(t)\in\mathbb{R}$ and
a new regressor $\mu(t)\in\mathbb{R}^{q}$, such that the new vector LRE

$e=\mu^{\top}\theta$

holds. In this case $\mu$ consists of the first $q$ elements of the last row
of the $(q+1)\times(q+1)$ matrix solution of the differential equation

$\dot{\Phi}=\begin{bmatrix}0_{q\times q}&u_{1}\\\
u_{2}\Delta^{\top}&u_{3}\end{bmatrix}\Phi,\;\Phi(0)=I_{q+1},$

with arbitrary functions $u_{1}(t)\in\mathbb{R}^{q}$, $u_{2}(t)\in\mathbb{R}$
and $u_{3}(t)\in\mathbb{R}$. The problem is that selecting the signals $u_{i}$
to stabilize this new LTV system seems a daunting task. In any case, given the
availability of DREM, the interest of such a study is highly questionable.

Another interesting line of research that we are pursuing now is to apply, to
the general problem of state observation, the idea of virtual signal
injection—that combined with GPEBO—is exploited in this paper. Indeed, the key
step in the construction of our LRE generator is captured in (8), that may be
interpreted as a virtual signal injection. Some encouraging preliminary
results have been developed and we expect to report them in the near future.

## References

  * [1] S. Aranovskiy, A. Belov, R. Ortega, N. Barabanov and A, Bobtsov, Parameter identification of linear time-invariant systems using dynamic regressor extension and mixing. International Journal of Adaptive Control and Signal Processing, vol. 33, no. 6, pp. 1016–1030, 2019.
  * [2] S. Aranovskiy, A. Bobtsov, R. Ortega and A. Pyrkin, Performance enhancement of parameter estimators via dynamic regressor extension and mixing, IEEE Trans. Automatic Control, vol. 62, pp. 3546–3550, 2017. (See also arXiv:1509.02763 for an extended version.)
  * [3] N. E. Barabanov and R. Ortega, On global asymptotic stability of $\dot{x}=\phi(t)\phi^{\top}(t)x$ with $\phi(t)$ bounded and not persistently exciting, Systems & Control Letters, vol. 109, pp. 24–27, 2017.
  * [4] N. Cho, H. Shin, Y. Kim and A. Tsourdos, Composite model reference adaptive control with parameter convergence under finite excitation, IEEE Trans. Automatic Control, vol. 63, pp. 811–818, 2018.
  * [5] G. Chowdhary, T. Yucelen, M. Muhlegg and E. N. Johnson, Concurrent learning adaptive control of linear systems with exponentially convergent bounds, International Journal of Adaptive Control and Signal Processing, vol. 27, no. 4, pp. 280–301, 2013.
  * [6] Y. Cui, J. Gaudio and A. Annaswamy, A new algorithm for discrete-time parameter estimation, ArXiv Preprint, 2021. (arXiv:2103.16653)
  * [7] P. Ioannou and J. Sun, Robust Adaptive Control, Prentice-Hall, NJ, 1996.
  * [8] G. Kreisselmeier and G. Rietze-Augst, Richness and excitation on an interval—with application to continuous-time adaptive control, IEEE Trans. Automatic Control, vol. 35, no. 2, pp. 165–171, 1990.
  * [9] L. Ljung, System Identification: Theory for the User, Prentice Hall, New Jersey, 1987.
  * [10] K. Narendra and A. Annaswamy, Stable Adaptive Systems, Prentice Hall, New Jersey, 1989.
  * [11] R. Ortega, A. Bobtsov, A. Pyrkin and A. Aranovskyi: A parameter estimation approach to state observation of nonlinear systems, Systems & Control Letters, vol. 85, pp 84–94, 2015.
  * [12] R. Ortega, A. Bobtsov, N. Nikolaev, J. Schiffer and D. Dochain, Generalized parameter estimation-based observers: Application to power systems and chemical-biological reactors, Automatica, vol. 129, 109635, 2021.
  * [13] R. Ortega, S. Aranovskiy, A. Pyrkin, A Astolfi and A. Bobtsov, New results on parameter estimation via dynamic regressor extension and mixing: Continuous and discrete-time cases, IEEE Trans. Automatic Control, vol. 66, no. 5, pp. 2265–2272, 2021.
  * [14] R. Ortega, V. Nikiforov and D. Gerasimov, On modified parameter estimators for identification and adaptive control: A unified framework and some new schemes, IFAC Annual Reviews in Control, vol. 50, pp. 278–293, 2020.
  * [15] Y. Pan and H. Yu, Composite learning robot control with guaranteed parameter convergence, Automatica, vol. 89, pp. 398–406, 2018.
  * [16] L. Praly, Convergence of the gradient algorithm for linear regression models in the continuous and discrete-time cases, Int. Rep. MINES ParisTech, Centre Automatique et Systèmes, December 26, 2017.
  * [17] S.B. Roy, S. Bhasin and I.N. Kar, Combined MRAC for unknown MIMO LTI systems with parameter convergence, IEEE Trans. Automatic Control, vol. 63, pp. 283–290, 2018.
  * [18] S. Sastry and M. Bodson, Adaptive Control: Stability, Convergence and Robustness, Prentice-Hall, New Jersey, 1989.
  * [19] G. Tao, Adaptive control design and analysis. vol. 37. John Wiley & Sons, New Jersey, 2003.
  * [20] B. Yi, R. Ortega, D. Wu and W. Zhang, Orbital stabilization of nonlinear systems via Mexican sombrero energy pumping-and-damping injection, Automatica, vol. 112, 108861, 2020.


# Combining and Steganography of 3D Face Textures

Mohsen Moradi and Mohammad-Reza Rafsanjani-Sadeghi M. Moradi and M.-R. Sadeghi
are with the Department of Mathematics, Amirkabir University of Technology -
Tehran Polytechnic , Tehran, Iran, ( e-mail<EMAIL_ADDRESS>).

###### Abstract

One of the serious issues in communication between people is hiding
information from others, and the best way for this is deceiving them. Since
nowadays face images are mostly used in a three-dimensional format, in this
paper we are going to steganography 3-D face images, detecting which by
curious people will be impossible. As in detecting the face, only its texture
is important, we separate texture from shape matrices, for eliminating half of
the extra information, steganography is done only for face texture, and for
reconstructing the 3-D face, we can use any other shape. Moreover, we will
indicate that, by using two textures, how two 3-D faces can be combined. For a
complete description of the process, first, 2-D faces are used as an input for
building 3-D faces, and then 3-D face and texture matrices are extracted
separately from the constructed 3-D face. Finally, 3-D textures are hidden
within other images.

###### Index Terms:

steganography, shape, texture, face image, combining images.

##  I Introduction

In cryptography, encrypted message is at the center of people’s attention, and
its security is based on the difficulties on access to the message key. In
steganography, other’s unawareness is used for hiding the message to send it
in the safest way. For this reason, first, the essential information of one
image which we want to send to others will be embedded in another host image
in a way that others cannot understand outward discrepancy of the initial host
image and the embedded host image, therefore legal receiver can extract and
reconstruct the initial embedded image. The most important researches
conducted in this area are [2, 15, 1, 7, 10].

Embedding capacity and visual quality are the two essential parameters in the
stego or cover images [4]. Embedding capacity refers maximum amount of
secretive message which can be embedded in the host image, and visual quality
is embedding a message in the host image in a way that human eye cannot notice
any difference between the new form and the original one. One criterion which
is usually used for evaluating visual quality is the peak of the signal to
noise ratio (PSNR) between stego image and original host image and expressed
in dB unit. The bigger the PSNR, the higher visual quality of stego image. In
other words, it is more difficult for the eye to detect stego image than to do
so for host image [5].
.

In this paper, to construct the 3-D image of the intended face for
steganography within other images first, a colored 2-D face image will be
processed using Basel database as [9] to the 3-D model is fitted. 3-D face and
texture matrices are extracted separately from the constructed 3-D face. Now,
if this 3-D texture is used for any other 3-D shapes, it is possible to
recognize the person as well. For this reason, in steganography, the shape
matrix is not too important and we only use texture matrix. A wavelet-based
watermarking algorithm is used to enhance the secrecy. By using singular value
decomposition (SVD) and discrete wavelet transform (DWT), the information of a
3-D texture of an image will be hidden in a label image. This labeled image
can be any 2-D image e.g., a fingerprint, a shape, or a texture of other 3-D
images. First, label image is converted to the frequency domain, and SVD is
used on both of the original cover image and 3-D texture. In the following,
the two obtained singular values will be replaced with each other.

##  II  Morphable model

At first, 3-D morphable models had been introduced by Blanz and Vetter [2, 3].
They were applied successfully in computer’s images and graphics. A 3-DMM
includes separate shape and texture models which by themselves can build each
person’s shape and texture’s changes, respectively. An 3-DMM is used for one
group of 3-D scanned images.

Constructing 3-D models is very hard and time-consuming which Paysan and his
colleagues provided their 3-DMM for public using; Basel model uses 3-D models
of 100 men and 100 women who ranged from 8 to 64 [9].
An Iterative Multiresolution Dense 3-D Registration which has been introduced
by Rodriguez [14] will be studied in this section. Suppose that the ith vertex
of registered image be in the $(x_{i},y_{i},z_{i})$ point and its RGB color is
$(R_{i},G_{i},B_{i})$. It will be assumed that for any face N 3-D points had
been registered; therefore, one registered face in the face and texture
language can be shown as:

$S^{{}^{\prime}}=\begin{pmatrix}x_{1},\cdots,x_{N}\\\ y_{1},\cdots,y_{N}\\\
z_{1},\cdots,z_{N}\end{pmatrix}$ (1)

$T^{{}^{\prime}}=\begin{pmatrix}R_{1},\cdots,R_{N}\\\ G_{1},\cdots,G_{N}\\\
B_{1},\cdots,B_{N}\par\end{pmatrix}.$ (2)

Also, these points can be shown by one row or column vector of 3N length. For
example, the vertical vector of shape and texture are:

$S^{{}^{\prime}}=\begin{pmatrix}x_{1},\cdots,x_{N},y_{1},\cdots,y_{N},z_{1},\cdots,z_{N}\par\end{pmatrix}^{T},$
(3)
$T^{{}^{\prime}}=\begin{pmatrix}R_{1},\cdots,R_{N},G_{1},\cdots,G_{N},B_{1},\cdots,B_{N}\par\end{pmatrix}^{T}$
(4)

which N is the number of registered faces. We suppose that these two kinds of
point’s display are equivalent, and we will use any of these two formats when
needed.

Now, for a database which has 200 shapes and texture matrices, a linear
combination for shapes and textures will be used. These linear combinations
can be as follows:

$S=\sum_{i=1}^{200}\alpha_{i}.S_{i}\ \ \ \ \ \ \ \ \ \ \
T=\sum_{i=1}^{200}\beta_{i}.T_{i}\ .$ (5)

It is almost impossible that these combinations, despite the fact that they
include all the possible faces, be similar to a real face. If a convex
combination ($\sum\alpha_{i}=1,\sum\beta_{i}=1,\alpha_{i},\beta_{i}\in[0,1]$)
of faces be supposed, a face can also be obtained, but again, it is not
possible that the points which are far from convex area form actual face
points. Therefore, for any vector, a coefficient is needed to be allocated to
a probability distribution for the description of a face. This probability is
modeled by a Gaussian distribution, in which shapes and textures are
decorrelated, with a diagonal matrix. Suppose a Gaussian distribution allows
that subspace of face to be estimated by a smaller set of orthogonal basis
vectors which is calculated using principle component analysis (PCA) of the
test samples.

Principle component analysis is a statistical tool which transform shape or
texture so that covariance matrix will be diagonal (it means data are
decorrelated). In this section, using PCA for the shape will be studied. Using
it for the texture is done in a similar way. PCA is a transform in the vector
space which is used mostly for decreasing the dimension of the data sets. At
first, principle component analysis had been used by Karl Pearson in 1901
[11]. This analysis includes the decomposition of the eigenvalues of the
covariance matrix.
The average of the shapes is calculated as:

$\overline{S}=\dfrac{1}{200}\sum_{i=1}^{200}S_{i}\ .$ (6)

By subtracting each sample shape from average shape matrix, the vertical
vector $a_{i}$ can be calculated as:

$a_{i}=vec(S_{i}-\overline{S})\ .$ (7)

These vertical vectors are used as the columns of matrix A, and the eigenvalue
vectors of a covariance C is calculated by a singular value decomposition
[12]. So we have:

$A=(a_{1},a_{2},\cdots,a_{200})=UWU^{T}$ (8)

and

$C=\dfrac{1}{200}AA^{T}=\dfrac{1}{200}UW^{2}U^{T},$ (9)

$Vec(S)$ will change matrix $S$ to a column vector by concatenating its
columns in a vertical order. 200 columns of orthogonal matrix $U$ are
eigenvectors of the covariance matrix $C$ and
$\sigma_{i}^{2}=\dfrac{\lambda_{i}^{2}}{200}$ is its eigenvalues which
$\lambda_{i}$’s are diagonal elements of matrix $W$ that are arranged in
decreasing order. The $i$th column of $U$ is shown by $U_{0,i}$ , and the
principle component of $i$ will be changed to a $3\times n$ matrix by
$S^{i}=U_{0,i}^{(3)}$. The notation $a_{m\times 1}^{(n)}$, changes the
$m\times 1$ vector a to a $n\times(m/n)$ matrix [8].

Now, instead of describing a new shape or texture as the linear combination of
the samples as equation 5, they can be expressed by the linear combination of
$n_{s}$ shapes and $n_{t}$ textures as principle components:

$S=\overline{S}+\sum_{i=1}^{n_{s}}\alpha_{i}S^{i},$ (10)

and

$T=\overline{T}+\sum_{i=1}^{n_{t}}\beta_{i}T^{i}\ .$ (11)

Therefore, the combination of an arbitrary number of shapes and textures for
constructing 3D morphable faces is used.
The 3DMMs should be so limited that improbable faces could be rarely sampled.

Experiments on the real data shows that supposing Gaussian distribution on the
face when sufficient information in the extraction is available, the result
will be satisfying. Anyhow, supposing Gaussian distribution is unexpected, so
Patel and Smith [10] presented another assumption. They observed that the
length of parameters’ vectors has Chebyshev distribution. In other words, the
real faces are on the spherical manifolds. In the model, the average model,
which its vector length is zero, is improbable, and the faces should have a
certain level of distinction.

##  III Extracting 3-D Shape and Texture From 2-D Face Images

In this section extracting three-dimensional shapes and textures of 2-D
colored images will be explained and we would show that how we can very well
combine two of them. The flowchart of this procedure is shown in Fig. 1.

Figure 1: Flowchart of extracting 3-D shapes and textures and combining them.

Bas [1] used also Basel database for fitting 2-D images on the 3-D model which
fitting $N$ 2-D $X_{i}=\begin{pmatrix}x_{i}\\\ y_{i}\end{pmatrix}$ points for
$i=1,\cdots,N$ to the 3D $V=\begin{pmatrix}u\\\ v\\\ w\end{pmatrix}$ vectors
on the model for getting the best shape and pose which can minimize the
Euclidean distance of the $N$ $X_{i}$ points from scaled orthogonal projection
of four components $(V,R,t,s)$ which $R$ is a 3 by 3 rotation matrix of real
numbers, $t$ is an translation 2-tuple of real numbers, and also $s$ scale is
a real number which is a positive real number, are three pose parameters. By
fitting 2-D image on a 3-D image, the 3-D shape $S$ and texture $T$ of the 3-D
face are obtained as equation 1 and 2 which by using shapes and textures, the
3-D faces can be constructed. Now suppose that two 2-D face images, which can
be chosen completely by accident, and their obtained 3-D images are as Fig. 2
.

Figure 2: 2-D face images and their obtained 3-D models.

The shape and texture of Obama’s face are $S_{1},T_{1}$, respectively, and the
shape and texture of Trump’s image are $S_{2}.T_{2}$ , respectively.

Now, if we use Obama’s texture which is a 3-D matrix $T_{1}$, and Trump’s 3-D
shape matrix which is shown by $S_{2}$, the Fig. 3 will be obtained.

Figure 3: Obama’s texture on Trump’s shape matrix

As it can be seen in the Fig. 3, the 3-D shape for detecting a person’s face
is not so important, and only the face’s texture for reconstructing human face
is sufficient. For this reason, if we want to steganography a 3-D face image,
it would be better to ignore half of the information, which is the shape
matrix and only use the texture matrix which will be described in the section
IV.

Moreover, if we want to combine the obtained 3-D images from the 2-D images
with each other, an appropriate way is using coefficients of principle
components. For $53490\times 3$ dimensional shape matrix $S_{1}$, by using
$PCA$ with Alternating Least Squares $(ALS)$ algorithm [6] we have:

$[coeffS_{1},scoreS_{1},\mu]=pca(S_{1})$ (12)

that $3\times 3$ dimensional coefficient $coeffS_{1}$ is obtained which its
columns are arranged in an increasing order and include any of the 3 principle
component coefficients. The $53490\times 3$ dimensional score matrix is also
shown by $scoreS_{1}$. The $1\times 3$ dimensional average of any variables’
matrix is shown by $\mu_{S_{1}}$. In the same way, by principle component
analysis and using $ALS$ algorithm, it can be done on the $S_{2},T_{1},T_{2}$.
Now, if the new shapes and textures are obtained as equations 13, the 3D face
models of Fig. 4 can be constructed.

$\displaystyle T_{1}=scoreT_{1}\ast
coeffT_{2}^{{}^{\prime}}+repmat(\mu_{2},53490,1),$ (13) $\displaystyle
T_{2}=scoreT_{2}\ast coeffT_{2}^{{}^{\prime}}+repmat(\mu_{2},53490,1),$
$\displaystyle S_{1}=scoreS_{2}\ast coefS_{2}^{{}^{\prime}}+repmat(\mu
S_{2},53490,1),$ $\displaystyle S_{2}=scoreS_{2}\ast
coefS_{2}^{{}^{\prime}}+repmat(\mu S_{2},53490,1)\ .$

Figure 4: Obama’s and Trump’s transformed shapes and textures as equations 13

To construct the left face in Fig. 4, the new shape and texture matrices
$S_{1},T_{1}$ have been used, and to construct the right face in Fig. 5, the
new shape and texture matrices $S_{2},T_{2}$ have been used.
Now by the transform which has been applied in 13, only by taking a simple
average in 14, we can see a satisfactory behavior of two politicians as the
face in the Fig. 5.

$\displaystyle S_{m}=(S_{1}+S_{2})/2$ (14) $\displaystyle
T_{m}=(T_{1}+T_{2})/2$

Figure 5: Combined Obama and Trump 3-D models

In Fig. 6 some other input 2-D colored face images and their combined 3-D
models are shown.

Figure 6: 3-D combining some arbitrary face images

##  IV  3-D Face Image Steganography

As it is shown in the previous section, what is essential in the constructing
of a 3-D face images and needs to be masked, is half of the face information
videlicet, the texture of the 3-D face image, and by using any other shape the
face can be reconstructed significantly same as the original one.

As it has been said, a texture matrix is a $53490\times 3$ dimensional matrix.
To hide this 3-D image, we add some zeros to the end of this matrix to obtain
a $53824\times 3$ dimension matrix which we have $232\times 232=53824$. By
using the second dimension of the resulted matrix, we change it to three
vertical $53824\times 1$ vectors, and after that we reshape each of these
vectors to a $232\times 232$ square matrix as the symbols of red, green, and
blue colors, respectively. The values of these matrices which might not be
number will be replaced by zero (also for better results, other techniques
could be used by adjacent points to recover the missing texture points ). Now
by using singular value decomposition, any of the colors is decomposed to a
3-tuples like $(U_{t},S_{t},V_{t})$, by subindices $i=r,g,b$ which are related
to their colors.

Also, a colored image is used as the cover image which its size will be
changed into to a $464\times 464\times 3$ dimensional matrix that we show it
by $C$ ( the cover matrix can also be any other 3-D texture matrix that in
this way we do as mentioned above). The matrix $C$, by using a single-level
discrete 2-D wavelet transform, is decomposed to a 4-tuple $(CA,CH,CV,CD)$ by
same $232\times 232\times 3$ dimensions which $CA$ is the approximation
coefficient matrix which has low frequency, and 3 details coefficient matrices
by the names horizontal matrix $CH$, vertical matrix $CV$ which has medium
frequency, and diagonal matrix $CD$ which has high frequency. Now, because of
the third dimension of the $CD$, we decompose it to 3 $CD_{r},CD_{g},CD_{b}$
matrices by $232\times 232$ dimension for the symbol of red, green, and blue
colors, respectively. In the following, by using singular decomposition, each
color matrices, like the matrices which were obtained by texture, is
decomposed to 3-tuples $(Uc,Sc,Vc)$ by subindices belongs to the three colors.

Now, for red, green, and blue colors, the obtained singular value matrix from
the cover image is replaced by the $S_{c_{i}}$ of that color $i$, added to the
tenth of $S_{t_{i}}$, for example for the red color we can obtain new value of

$CD_{new_{r}}=Uc_{r}\times(Sc_{r}+0.1\ast St_{r})\times Vc_{r}^{T}.$ (15)

Then by concatenation of these 3 red, green, and blue matrices, the new value
of $D_{new}$ will be obtained and by inverse wavelet decomposition on the
4-tuple $(CA,CH,CV,CD_{new})$ the stego image will be constructed. The
procedure of steganography which is described in this section can be seen in
Fig. 7.

Figure 7: Steganography of 3-D texture of face images which $i$ is equal to
three red, green, and blue colors.

To extract the hidden image we should do as the inverse of abovementioned
process, and at first 2-D wavelet decomposition for rebuilding 4-tuple
$(wA,wH,wV,wD)$ is performed on the stego image. Because of the third
dimension of $wD$, we decompose $wD$ to the three red, green, and blue color
images, and on any color by singular value decomposition, the 3-tuple
$(U_{e_{i}},S_{e_{i}},V_{e_{i}})$ by three subindices $i=r,g,b$ which belong
to any color, is obtained which for calculating $S_{t_{i}}$ we should subtract
that value of $S_{c_{i}}$ from obtained singular value $S_{e_{i}}$ and
multiply it to 10 ( the inverse of what mentioned before). Next, the new
texture column components from the result of $U_{e_{i}}\times S_{t_{i}}\times
V_{e_{i}}^{T}$ for each color will be calculated. This procedure is shown in
the Fig. 8.

Figure 8: Extracting hidden texture for constructing 3-D face image which $i$
is equal to three red, green, and blue colors.

We can implement the resulted texture on any shape to reconstruct the desired
3-D face image. As an example, the images for the mentioned process are as
Fig. 9 which from left to right are the original cover image, 3-D used
original image, stego image, and desired extracted the image, respectively. A
known used criteria for measuring image quality is the peak of the signal to
noise ratio [13]. The peak of the signal to noise ratio of two original cover
image and stego image of Fig. 9 is $81.51$ (dB) by equation 16 which MSE is
the Mean Square Error and S is the maximum pixel value. A known used criteria
for measuring image quality is peak of signal to noise ratio [13]. The peak of
signal to noise ratio of two original cover image and stego image of Fig. 9 is
$81.51$ (dB) by equation 16 which MSE is the Mean Square Error and $S$ is the
maximum pixel value.

$PSNR=-10log_{10}\dfrac{MSE}{S^{2}}$ (16)

Figure 9: Steganography of 3-D texture which from left to right are the
original cover image, 3-D used original image, stego image, and desired
extracted the image, respectively.

Also, some other covers and hidden images with their PSNR values are shown in
the Fig. 10. In all results, the destruction of cover image by hiding texture
in it from our method is very insignificant which all the PSNR values are near
to 80. A similar result is achieved by Chao [16].

Figure 10: Steganography and extracting of three other textures with PSNR
equal to 79.4391, 81.3087, and 80.9286 dB, respectively.

##  V Conclusion

Hiding the secretive information from the enemies’ view in a way that they do
not be at all skeptical to the existence of hidden messages helps to the safe
communication. In this paper, by using the idea of extracting 3D face texture,
this matrix is embedded in other images, and it has been indicated that how by
using only half of the face information, which is texture, we can
appropriately reconstruct face images. In addition, it has been shown that how
by combining two textures and without any stricture on shape matrix, we can
combine two faces pleasantly.

## References

  * [1] Bas, Anil, et al. "Fitting a 3D morphable model to edges: A comparison between hard and soft correspondences." arXiv preprint arXiv:1602.01125 (2016).
  * [2] Blanz, Volker, and Thomas Vetter. "A morphable model for the synthesis of 3D faces." Proceedings of the 26th annual conference on Computer graphics and interactive techniques. ACM Press/Addison-Wesley Publishing Co., 1999.
  * [3] Blanz, Volker, and Thomas Vetter. "Face recognition based on fitting a 3D morphable model." IEEE Transactions on pattern analysis and machine intelligence 25.9 (2003): 1063-1074.
  * [4] Chang, C. C., T. S. Chen, and K. F. Hwang. "Electronic Image Techniques." Taipei: Unalis (2000).
  * [5] Johnson, N.F., Jajodia, S., 1998. Exploring steganography: Seeing the unseen. IEEE Comput. (February), 26–64
  * [6] Kiers, Henk AL. "Setting up alternating least squares and iterative majorization algorithms for solving various matrix optimization problems." Computational statistics and data analysis 41.1 (2002): 157-170.
  * [7] Kittler, Josef, et al. "3D Morphable Face Models and Their Applications." International Conference on Articulated Motion and Deformable Objects. Springer International Publishing, 2016.
  * [8] Minka, Thomas P. "Automatic choice of dimensionality for PCA." Nips. Vol. 13. 2000.
  * [9] P. Paysan, R. Knothe, B. Amberg, S. Romdhani, and T. Vetter, “A 3D face model for pose and illumination invariant face recognition,” in Proc. IEEE Intl. Conf. on Advanced Video and Signal based Surveillance, 2009.
  * [10] Patel, Ankur, and William AP Smith. "3d morphable face models revisited." Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on. IEEE, 2009.
  * [11] Pearson, Karl. "Principal components analysis." The London, Edinburgh and Dublin Philosophical Magazine and Journal 6.2 (1901): 566.
  * [12] Press, William H., et al. "Singular value decomposition." Solution of Linear Algebraic Equation from Numerical Recipes in C (1992): 59-70.
  * [13] Rafael, C. "Gonzalez, and Richard E. Woods." Digital Image Processing (1992).
  * [14] Rodriguez, JR Tena. "3D face modelling for 2D+ 3D face recognition." (2007).
  * [15] Zhu, Xiangxin, and Deva Ramanan. "Face detection, pose estimation, and landmark localization in the wild." Computer Vision and Pattern Recognition (CVPR), 2012 IEEE Conference on. IEEE, 2012.
  * [16] Chao, Min-Wen, et al. "A high capacity 3D steganography algorithm." IEEE transactions on visualization and computer graphics 15.2 (2009): 274-284.