Automation of reversible steganographic coding with nonlinear discrete optimisation

Abstract

Authentication mechanisms are at the forefront of defending the world from various types of cybercrime. Steganography can serve as an authentication solution through the use of a digital signature embedded in a carrier object to ensure the integrity of the object and simultaneously lighten the burden of metadata management. Nevertheless, despite being generally imperceptible to human sensory systems, any degree of steganographic distortion might be inadmissible in fidelity-sensitive situations such as forensic science, legal proceedings, medical diagnosis and military reconnaissance. This has led to the development of reversible steganography. A fundamental element of reversible steganography is predictive analytics, for which powerful neural network models have been effectively deployed. Another core element is reversible steganographic coding. Contemporary coding is based primarily on heuristics, which offers a shortcut towards sufficient, but not necessarily optimal, capacity–distortion performance. While attempts have been made to realise automatic coding with neural networks, perfect reversibility is unattainable via such learning machinery. Instead of relying on heuristics and machine learning, we aim to derive optimal coding by means of mathematical optimisation. In this study, we formulate reversible steganographic coding as a nonlinear discrete optimisation problem with a logarithmic capacity constraint and a quadratic distortion objective. Linearisation techniques are developed to enable iterative mixed-integer linear programming. Experimental results validate the near-optimality of the proposed optimisation algorithm when benchmarked against a brute-force method.

Keywords:

• Automatic coding
• mathematical optimisation
• reversible steganography

1. Introduction

Steganography is the art and science of concealing information within a carrier object (Anderson & Petitcolas, 1998). The term encompasses a wide range of techniques and applications, including but not limited to covert communications (Fridrich et al., 2005), ownership identification (Cox et al., 1997), copyright protection (Barni et al., 1998), broadcast monitoring (Depovere et al., 1999) and traitor tracing (He & Wu, 2006). An important application of steganography is data authentication, which plays a vital role in cybersecurity. The advent of data-centric artificial intelligence has been accompanied by cybersecurity concerns (Boden et al., 2017). It has been reported that intelligent systems are vulnerable to adversarial attacks such as poisonous data collected for re-training during deployment (Muñoz-González et al., 2017), malware codes hidden in neural network parameters (Liu et al., 2020) and invisible perturbations crafted to cause erroneous decisions (Goodfellow et al., 2015). A proper authentication mechanism must ensure that the integrity of data has not been undermined and that the identity of users has not been forged, and thereby protect against these insidious threats.

Digital signatures are a type of authentication technology that is based upon modern cryptography (Rivest et al., 1978). This technology can be incorporated into a trustworthy surveillance camera in such a way that photographs are taken and stored along with digital signatures (Friedman, 1993). However, storing such auxiliary metadata as a separate file entails the risk of accidental loss and mismanagement during the data lifecycle. Steganography can allow auxiliary information about the data to be embedded invisibly within the data itself. Nevertheless, although generally imperceptible to human sensory systems, any degree of steganographic distortion might not be admissible in some fidelity-sensitive situations such as forensic science, legal proceedings, medical diagnosis and military reconnaissance. This is where the notion of reversible computing comes into play (Alattar, 2004; Chang et al., 2018; Coatrieux et al., 2013; De Vleeschouwer et al., 2003; Fridrich et al., 2001; Lee et al., 2007; Wu & Zhang, 2020).

A fundamental element of reversible steganography, in common with lossless compression, is predictive modelling (Rissanen, 1984; Shannon, 1948; Weinberger & Seroussi, 1997). Prediction error modulation is a cutting-edge reversible steganographic technique composed of a predictive analytics module and a reversible coding module (Celik et al., 2005; Dragoi & Coltuc, 2014; Fallahpour, 2008; Hwang et al., 2016; Li et al., 2011; Sachnev et al., 2009; Thodi & Rodriguez, 2007). The recent development of deep learning has advanced the frontier of reversible steganography. It has been shown that deep neural networks can be applied as powerful predictive models (Chang, 2020, 2021, 2022; Hu & Xiang, 2021). Despite inspiring progress in the analytics module, the design of the coding module is still based largely on heuristics. While there are studies on end-to-end deep learning that use neural networks for automatic reversible computing, perfect reversibility cannot be guaranteed (Duan et al., 2019; Lu et al., 2021; Z. Zhang et al., 2019). From a certain point of view, it is hard for a neural network, as a monolithic black box, to follow the intricate procedures of reversible computing (Castelvecchi, 2016). While deep learning is adept at handling the complex nature of the real world (LeCun et al., 2015), reversible computing is more of a mechanical process in which procedures have to be conducted in accordance with rigorous algorithms. Therefore, at the time of writing, it seems advisable to follow a modular framework.

The essence of reversible steganographic coding is determining how values change to represent different message digits. Different solutions can lead to different trade-offs between capacity and distortion. Instead of relying on heuristics, this study pursues the development of optimal coding for reversible steganography in order to attain optimal capacity–distortion performance. We model reversible steganographic coding as a mathematical optimisation problem and propose an optimisation algorithm for addressing the nonlinearity of this problem. In particular, the task is to minimise steganographic distortion subject to a capacity constraint, where both objective and constraint are nonlinear functions. We propose linearisation techniques for addressing this nonlinear discrete optimisation problem. The remainder of this paper is organised as follows. Section 2 outlines the background regarding reversible steganography. Section 3 formulates the nonlinear discrete optimisation problem and discusses the complexity of a brute-force search algorithm. Section 4 presents linearisation techniques for tackling the nonlinear discrete optimisation problem. Section 5 analyses the optimality of solutions through simulation experiments. Section 6 provides concluding remarks.

2. Background

Prediction error modulation is a reversible steganographic technique that consists of an analytics module and a coding module. The analytics module begins by splitting a cover image into context and query sets, denoted by $\mathit{c}$ and $\mathit{q}$ respectively. A conventional method is to arrange pixels in two groups according to a chequered pattern. Then a predictive model is applied to predict the intensities of the query pixels from the intensities of the context pixels. A contemporary practice of predictive modelling is to employ an artificial neural network originally designed for computer vision tasks. The coding module embeds a message $\mathit{\omega }$ into the cover image by modulating the prediction errors $\mathit{ϵ}=\mathit{q}-\tilde{\mathit{q}}$. The modulated errors ${\mathit{ϵ}}^{\mathrm{\prime }}$ are then added to the predicted intensities, causing distortion of the query pixels. The stego image is created by merging the context set $\mathit{c}$ and the modulated query set ${\mathit{q}}^{\mathrm{\prime }}$. The decoding procedure is similar to the encoding procedure. It begins by predicting the query pixel intensities. Since the context set is kept unchanged, the prediction in the decoding phase is guaranteed to be identical to that in the encoding phase given the same predictive model. The message is extracted and the query set is recovered by demodulating the prediction errors. The image is reversed to its original state by merging the context and recovered query sets. The procedures for encoding and decoding are depicted schematically in Figure 1 and also provided in Algorithms 1 and 2. We would like to note that the message may contain certain auxiliary information for handling pixel intensity overflow. This paper does not go into detail about every aspect of the stego-system; instead, our study focuses on the mathematical optimisation of reversible steganographic coding.

Figure 1. Workflow of reversible steganography with prediction error modulation.

3. Nonlinear discrete optimisation

The essence of reversible steganographic coding is designating one or multiple error values as the carrier and determining how these values change to represent different message digits. A rule of thumb for reversible steganographic coding is to choose the prediction errors of the peak frequency as the carrier. While the peak frequency implies the highest capacity, this capacity-greedy strategy is not necessarily optimal in terms of minimising distortion.

3.1. Problem definition

According to the typical law of error, the frequency of an error can be expressed as an exponential function of its numerical magnitude, disregarding sign (Wilson, 1923). In other words, the frequency distribution of prediction errors is expected to centre around zero. In general, a smaller absolute error tends to have a higher occurrence. A special exception is that the occurrence of zero might be lower than the occurrence of a certain small absolute error considering that the latter is the sum of both positive and negative error occurrences. Consider an absolute error histogram as shown in Figure 2. The problem of reversible steganographic coding is to establish a mapping between the values in $[0,n]$ and the values in $[0,n+\vartheta ]$, where ϑ denotes the extra quota and is typically defined as less than or equal to the number of successive empty bins in the absolute error histogram. Encoding is a one-to-many mapping that links a cover value to one or more stego values. A message digit can only be represented if the connections are greater than one. Different cover values can never yield the same stego value. This is done in order to avoid an overlap between values (i.e. an ambiguity in decoding). Therefore, a cover value may be changed to a different stego value even if it does not represent any message digit. We impose a constraint that each cover value can only be mapped to the nearest available stego values since a non-cross mapping drastically reduces the problem dimension. An example of a cover/stego mapping is illustrated in Figure 3.

Figure 2. Example of absolute error distribution with highlighted zero occurrences.

Figure 3. Example of reversible steganographic coding.

3.2. Model formulation

Let us denote by $a_i$ the frequency of the value i and by $x_i$ the number of extra cover-to-stego links for the value i. The total number of links for the value i equals $x_i+1$. The number of bits that can be represented by modifying the value i is ${\log }_2(x_i+1)$ and thus the capacity is computed by (1) $\begin{array}{r}\mathfrak{C}=\sum _{i=0}^n{a}_i{\log }_2(x_i+1).\end{array}$(1) In fact, the number of bits that can be represented by modifying 0 equals ${\log }_2(2x_0+1)$ because 0 can be mapped to both positive and negative values. For example, there are three different states 0 and $\pm 1$ when $x_0=1$, and five different states 0, $\pm 1$ and $\pm 2$ when $x_0=2$. To be concise, we simplify the case by mapping 0 randomly to a positive or negative value so that the capacity computation for 0 is identical to that for other values at the cost of slightly underestimating the capacity offered by the former. The probability of changing a cover value to each stego value is $1/(x_i+1)$. The deviations of the first to the last stego value are $0+y_i$ to $x_i+y_i$ respectively, where $y_i$ denotes the sum of all the previous extra links (i.e. the cumulative deviation). Hence, the expected distortion in terms of the squared deviations is computed by (2) $\begin{array}{r}\mathfrak{D}=\sum _{i=0}^n{a}_i\left(\frac{(0+y_i{)}^2+(1+y_i{)}^2+\cdots +(x_i+y_i{)}^2}{x_i+1}\right),\end{array}$(2) where (3) $\begin{array}{r}{y}_i=\sum _{j=0}^{i-1}{x}_j.\end{array}$(3) We can simplify the algebraic expression by (4) $\begin{array}{rl}& \frac{(0+y_i{)}^2+(1+y_i{)}^2+\cdots +(x_i+y_i{)}^2}{x_i+1}\\ & =\frac{(0^2+2y_i\cdot 0+y_i^2)+(1^2+2y_i\cdot 1+y_i^2)+\cdots +(x_i^2+2y_i\cdot x_i+y_i^2)}{x_i+1}\\ & =\frac{(0^2+1^2+\cdots +x_i^2)+2y_i(0+1+\cdots +x_i)+y_i^2(x_i+1)}{x_i+1}\\ & =\frac{x_i(x_i+1)(2x_i+1)}{6(x_i+1)}+\frac{2y_i{x}_i(x_i+1)}{2(x_i+1)}+\frac{y_i^2(x_i+1)}{x_i+1}\\ & =\frac{1}{3}{x}_i^2+\frac{1}{6}{x}_i+x_i{y}_i+y_i^2.\end{array}$(4) The reason for computing squared deviations rather than absolute deviations is that image quality is often measured by the peak signal-to-noise ratio (PSNR), which is defined via the mean squared error (MSE). Our goal is to solve for the decision variables $x_i\in \{0,\dots \vartheta \}$ which minimise the distortion objective subject to the capacity constraint. The sum of all the extra cover-to-stego links is not allowed to exceed the quota ϑ. To summarise, the mathematical optimisation problem for reversible steganographic coding is $\begin{array}{rl}min& \mathfrak{D}=\sum _{i=0}^n{a}_i(\frac{1}{3}{x}_i^2+\frac{1}{6}{x}_i+x_i{y}_i+y_i^2),\\ \mathrm{s}.\mathrm{t}.& \mathfrak{C}=\sum _{i=0}^n{a}_i{\log }_2(x_i+1)\ge \mathrm{p}\mathrm{a}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d},\\ & \sum _{i=0}^n{x}_i\le \vartheta ,\\ \mathrm{v}\mathrm{a}\mathrm{r}.& x_i\in \{0,\dots ,\vartheta \},\mathrm{\forall }i=0,\dots ,n.\end{array}$

3.3. Brute-fORCE search

Brute-force search is a baseline method for benchmarking optimisation algorithms. The solution space that exhausts all possible combinations of the decision variables is equal to $(\vartheta +1{)}^{n+1}\in \mathcal{O}(c^n)$. By taking account of the quota constraint, we can reduce the solution space from the number of possible combinations to the number of feasible combinations. In number theory and combinatorics, the partition function $\mathrm{part}\left(t\right)$ computes the number of ways of writing t as a sum of the positive integers in $[1,t]$. Let ${\mathbf{\Lambda }}_t$ denote a matrix of $\mathrm{part}\left(t\right)$ rows and t columns which enumerates all possible partitions: (5) $\begin{array}{r}{\mathbf{\Lambda }}_t=\left[\begin{array}{c}{\mathit{\lambda }}_1\\ ⋮\\ {\mathit{\lambda }}_{\mathrm{part}\left(t\right)}\end{array}\right]=\left[\begin{array}{ccc}{\lambda }_{1,1}& \cdots & {\lambda }_{1,t}\\ ⋮& \ddots & ⋮\\ {\lambda }_{\mathrm{part}\left(t\right),1}& \cdots & {\lambda }_{\mathrm{part}\left(t\right),t}\end{array}\right].\end{array}$(5) Each vector ${\mathit{\lambda }}_{\ell }$ represents a possible partition in which each element is the quantity of a candidate integer (i.e. the summand). For example, ${\mathbf{\Lambda }}_2$, ${\mathbf{\Lambda }}_3$ and ${\mathbf{\Lambda }}_4$ are $\begin{array}{rl}ccc\text{{1}}& \text{{2}}\\ \left[cc\right]c2& 0& {\mathit{\lambda }}_1\\ 0& 1& {\mathit{\lambda }}_2\\ ,cccc\text{{1}}& \text{{2}}& \text{{3}}\\ \left[ccc\right]c3& 0& 0& {\mathit{\lambda }}_1\\ 1& 1& 0& {\mathit{\lambda }}_2\\ 0& 0& 1& {\mathit{\lambda }}_3\\ ,ccccc\text{{1}}& \text{{2}}& \text{{3}}& \text{{4}}\\ \left[cccc\right]c4& 0& 0& 0& {\mathit{\lambda }}_1\\ 2& 1& 0& 0& {\mathit{\lambda }}_2\\ 0& 2& 0& 0& {\mathit{\lambda }}_3\\ 1& 0& 1& 0& {\mathit{\lambda }}_4\\ 0& 0& 0& 1& {\mathit{\lambda }}_5\end{array}.$The total number of feasible solutions can be calculated by adding up the number of feasible solutions given by each individual partition matrix from ${\mathbf{\Lambda }}_1$ to ${\mathbf{\Lambda }}_{\vartheta }$; that is, (6) $\begin{array}{r}\sum _{t=1}^{\vartheta }\mathrm{feasible}({\mathbf{\Lambda }}_t,n^{\ast }),\end{array}$(6) where $n^{\ast }=n+1$ denotes the number of integers in $[0,n]$. For each matrix ${\mathbf{\Lambda }}_t$, the number of feasible solutions is computed by summing the number of possible combinations given by each partition vector ${\mathit{\lambda }}_{\ell }$, denoted by (7) $\begin{array}{r}\mathrm{feasible}({\mathbf{\Lambda }}_t,n^{\ast })=\sum _{\ell =1}^{\mathrm{part}\left(t\right)}\mathrm{comb}({\mathit{\lambda }}_{\ell },n^{\ast }).\end{array}$(7) A combination is a selection of values from a set of $n^{\ast }$ values based on a given partition vector and hence the number of combinations is computed by (8) $\begin{array}{r}\mathrm{comb}({\mathit{\lambda }}_{\ell },n^{\ast })=\prod _{i=1}^t\left(\begin{array}{c}{\displaystyle n^{\ast }-\sum _{j=1}^{i-1}{\lambda }_j^{\ast }}\\ {\lambda }_i^{\ast }\end{array}\right),\end{array}$(8) where ${\lambda }_i^{\ast }={\lambda }_{\ell ,i}$ represents a convenient notation without explicitly writing out the index of the partition vector (for reducing the verbosity). The number of combinations is a product of t binomial coefficients and each term is meant to choose (and remove) an unordered subset of ${\lambda }_i^{\ast }$ values from the remaining values in the set of $n^{\ast }$ values. Let us take ${\mathbf{\Lambda }}_3$ for example. The number of combinations for partition vectors ${\mathit{\lambda }}_1$, ${\mathit{\lambda }}_2$ and ${\mathit{\lambda }}_3$ are computed as follows: $\begin{array}{rl}\mathrm{comb}({\mathit{\lambda }}_1,n^{\ast })& =(\genfrac{}{}{0ex}{}{n^{\ast }}{3})(\genfrac{}{}{0ex}{}{n^{\ast }-3}{0})(\genfrac{}{}{0ex}{}{n^{\ast }-3-0}{0}),\\ \mathrm{comb}({\mathit{\lambda }}_2,n^{\ast })& =(\genfrac{}{}{0ex}{}{n^{\ast }}{1})(\genfrac{}{}{0ex}{}{n^{\ast }-1}{1})(\genfrac{}{}{0ex}{}{n^{\ast }-1-1}{0}),\\ \mathrm{comb}({\mathit{\lambda }}_3,n^{\ast })& =(\genfrac{}{}{0ex}{}{n^{\ast }}{0})(\genfrac{}{}{0ex}{}{n^{\ast }-0}{0})(\genfrac{}{}{0ex}{}{n^{\ast }-0-0}{1}).\end{array}$The number of combinations can be approximated by (9) $\begin{array}{rl}& \prod _{i=1}^t\left(\begin{array}{c}{\displaystyle n^{\ast }-\sum _{j=1}^{i-1}{\lambda }_j^{\ast }}\\ {\lambda }_i^{\ast }\end{array}\right)\\ & =(\genfrac{}{}{0ex}{}{n^{\ast }}{{\lambda }_1^{\ast }})(\genfrac{}{}{0ex}{}{n^{\ast }-{\lambda }_1^{\ast }}{{\lambda }_2^{\ast }})\cdots (\genfrac{}{}{0ex}{}{n^{\ast }-{\lambda }_1^{\ast }-{\lambda }_2^{\ast }-\cdots -{\lambda }_{t-1}^{\ast }}{{\lambda }_t^{\ast }})\\ & =\frac{n^{\ast }!}{{\lambda }_1^{\ast }!(n^{\ast }-{\lambda }_1^{\ast })!}\times \frac{(n^{\ast }-{\lambda }_1^{\ast })!}{{\lambda }_2^{\ast }!(n^{\ast }-{\lambda }_1^{\ast }-{\lambda }_2^{\ast })!}\times \cdots \times \frac{(n-\sum _{j=1}^{t-1}{\lambda }_j^{\ast })!}{{\lambda }_t^{\ast }!(n-\sum _{j=1}^t{\lambda }_j^{\ast })!}\\ & =\frac{n^{\ast }!}{{\lambda }_1^{\ast }!{\lambda }_2^{\ast }!\dots {\lambda }_t^{\ast }!(n^{\ast }-\sum _{j=1}^t{\lambda }_j^{\ast })!}\\ & =\frac{n^{\ast }(n^{\ast }-1)(n^{\ast }-2)\dots (n^{\ast }-(\sum _{j=1}^t{\lambda }_j^{\ast }-1\left)\right)}{{\lambda }_1^{\ast }!{\lambda }_2^{\ast }!\dots {\lambda }_t^{\ast }!}\\ & \le \frac{n^{\ast }(n^{\ast }-1)(n^{\ast }-2)\dots (n^{\ast }-(t-1\left)\right)}{{\lambda }_1^{\ast }!{\lambda }_2^{\ast }!\dots {\lambda }_t^{\ast }!}\approx n^t.\end{array}$(9) Hence, the complexity of this brute-force algorithm is approximately equal to (10) $\begin{array}{r}\sum _{t=1}^{\vartheta }\sum _{\ell =1}^{\mathrm{part}\left(t\right)}\mathrm{comb}({\mathit{\lambda }}_{\ell },n^{\ast })\approx \sum _{t=1}^{\vartheta }\mathrm{part}\left(t\right)\cdot n^t\in \mathcal{O}\left(n^c\right).\end{array}$(10)

4. Linearisation

The difficulty of our optimisation problem lies in the nonlinear nature of the capacity constraint and the distortion objective. To apply off-the-shelf optimisation tools, we have to tackle these nonlinearities.

4.1. Logarithmic capacity constraint

The capacity constraint involves the calculation of logarithm of variables ${\log }_2(x_i+1)$. The logarithmic function is nonlinear. A useful linearisation trick is to re-model the problem with binary-integer variables. We binarise each decision variable $x_i$ with the domain $[0,\vartheta ]$ into a 0/1 vector (or a one-hot vector) of length $\vartheta +1$, as illustrated in Figure 4. The vector consists of 0s with the exception of a single 1 whose position indicates the value of $x_i$; that is, (11) $\begin{array}{r}{\mathbf{x}}_i=[{\mathrm{x}}_i^0,\dots ,{\mathrm{x}}_i^{\vartheta }]\in \{0,1{\}}^{\vartheta +1},\end{array}$(11) such that (12) $\begin{array}{r}\text{𝟙}\cdot {\mathbf{x}}_i^{⊺}=1,\mathrm{\forall }i=0,\dots ,n.\end{array}$(12) We can retrieve $x_i$ using the dot product of vectors (13) $\begin{array}{r}{x}_i=[0,\dots ,\vartheta ]\cdot {\mathbf{x}}_i^{⊺}=\mathbf{v}{\mathbf{x}}_i^{⊺}.\end{array}$(13) Accordingly, the quota constraint becomes (14) $\begin{array}{r}\sum _{i=0}^n\mathbf{v}{\mathbf{x}}_i^{⊺}\le \vartheta .\end{array}$(14) In a similar manner, the logarithm can be derived using the dot product of vectors (15) $\begin{array}{rl}{\log }_2(x_i+1)& =[{\log }_2(0+1),\dots ,{\log }_2(\vartheta +1\left)\right]\cdot {\mathbf{x}}_i^{⊺}\\ & ={\mathbf{v}}_{\log }{\mathbf{x}}_i^{⊺}.\end{array}$(15) Hence, we rewrite the capacity constraint as (16) $\begin{array}{r}\mathfrak{C}=\sum _{i=0}^n{a}_i{\mathbf{v}}_{\log }{\mathbf{x}}_i^{⊺}.\end{array}$(16)

Figure 4. Example of binary-integer decision variable.

4.2. Quadratic distortion objective

The distortion objective involves three nonlinear terms $x_i^2$, $y_i^2$ and $x_i{y}_i$. These terms are quadratic functions of variables. The first term can be approached using the dot product as before; that is (17) $\begin{array}{r}{x}_i^2=[0^2,\dots ,{\vartheta }^2]\cdot {\mathbf{x}}_i^{⊺}={\mathbf{v}}_{\mathrm{sq}}{\mathbf{x}}_i^{⊺}.\end{array}$(17) The remaining two terms contain the partial sum of variables $y_i$, which is computed by (18) $\begin{array}{r}{y}_i=\sum _{j=0}^{i-1}\mathbf{v}{\mathbf{x}}_j^{⊺}.\end{array}$(18) To linearise the univariate quadratic term $y_i^2$ and the bivariate quadratic term $x_i{y}_i$, we introduce two non-negative continuous slack variables $z_{y_i^2}\ge 0$ and $z_{x_i{y}_i}\ge 0$. Replacing the quadratic terms with the dot product and the slack variables results in a linear distortion objective (19) $\begin{array}{r}\mathfrak{D}=\sum _{i=0}^n{a}_i(\frac{1}{3}{\mathbf{v}}_{\mathrm{sq}}{\mathbf{x}}_i^{⊺}+\frac{1}{6}\mathbf{v}{\mathbf{x}}_i^{⊺}+z_{x_i{y}_i}+z_{y_i^2}).\end{array}$(19) We begin by solving this mixed-integer linear programming problem, which does not yet reflect the quadratic terms regarding cumulative distortion, and obtain an initial solution comprising ${\tilde{\mathbf{x}}}_i$, ${\tilde{z}}_{y_i^2}$, and ${\tilde{z}}_{x_i{y}_i}$. The initial slack variables would be zeros because the objective is to minimise distortion. To make the slack variables reflect the quadratic terms properly, we add the following constraints (20) $\begin{array}{r}\begin{array}{rl}{z}_{y_i^2}& \ge y_i^2,\\ z_{x_i{y}_i}& \ge x_i{y}_i.\end{array}\end{array}$(20) In this way, we reformulate a problem with a nonlinear objective into a problem with a linear objective and nonlinear constraints. We make use of the solution obtained previously to linearise these nonlinear constraints and solve the mixed-integer linear programming problem iteratively. To begin with, we express the variables in terms of the previous solution: (21) $\begin{array}{r}\begin{array}{rl}{x}_i& ={\tilde{x}}_i+{\delta }_{x_i},\\ y_i& ={\tilde{y}}_i+{\delta }_{y_i},\end{array}\end{array}$(21) where ${\tilde{x}}_i$ and ${\tilde{y}}_i$ are treated as constants. Then, we apply the first-order Taylor series to approximate the univariate quadratic term as (22) $\begin{array}{rl}f\left(y_i\right)& =f({\tilde{y}}_i+{\delta }_{y_i})\\ & =f\left({\tilde{y}}_i\right)+f^{\mathrm{\prime }}\left({\tilde{y}}_i\right){\delta }_{y_i}+\cdots \\ & ={\tilde{y}}_i^2+2{\tilde{y}}_i{\delta }_{y_i}+\cdots \\ & \approx {\tilde{y}}_i^2+2{\tilde{y}}_i(y_i-{\tilde{y}}_i)\\ & =2{\tilde{y}}_i{y}_i-{\tilde{y}}_i^2,\end{array}$(22) and similarly the bivariate quadratic term as (23) $\begin{array}{rl}f(x_i,y_i)& =f({\tilde{x}}_i+{\delta }_{x_i},{\tilde{y}}_i+{\delta }_{y_i})\\ & =f({\tilde{x}}_i,{\tilde{y}}_i)+\frac{\mathrm{\partial }f({\tilde{x}}_i,{\tilde{y}}_i)}{\mathrm{\partial }{x}_i}{\delta }_{x_i}+\frac{\mathrm{\partial }f({\tilde{x}}_i,{\tilde{y}}_i)}{\mathrm{\partial }{y}_i}{\delta }_{y_i}+\cdots \\ & ={\tilde{x}}_i{\tilde{y}}_i+{\tilde{y}}_i{\delta }_{x_i}+{\tilde{x}}_i{\delta }_{y_i}+\cdots \\ & \approx {\tilde{x}}_i{\tilde{y}}_i+{\tilde{y}}_i(x_i-{\tilde{x}}_i)+{\tilde{x}}_i(y_i-{\tilde{y}}_i)\\ & ={\tilde{x}}_i{y}_i+{\tilde{y}}_i{x}_i-{\tilde{x}}_i{\tilde{y}}_i.\end{array}$(23) As a result, the nonlinear constraints are transformed into linear constraints (24) $\begin{array}{r}\begin{array}{rl}2{\tilde{y}}_i{y}_i-z_{y_i^2}& \le {\tilde{y}}_i^2,\\ {\tilde{x}}_i{y}_i+{\tilde{y}}_i{x}_i-z_{x_i{y}_i}& \le {\tilde{x}}_i{\tilde{y}}_i.\end{array}\end{array}$(24) To recapitulate, the nonlinear discrete optimisation problem is approached by means of an iterative method that solves a mixed-integer linear programming problem with binary-integer variables and non-negative continuous slack variables: $\begin{array}{rl}min& \mathfrak{D}=\sum _{i=0}^n{a}_i(\frac{1}{3}{\mathbf{v}}_{\mathrm{sq}}{\mathbf{x}}_i^{⊺}+\frac{1}{6}\mathbf{v}{\mathbf{x}}_i^{⊺}+z_{x_i{y}_i}+z_{y_i^2}),\\ \mathrm{s}.\mathrm{t}.& \mathfrak{C}=\sum _{i=0}^n{a}_i{\mathbf{v}}_{\log }{\mathbf{x}}_i^{⊺}\ge \mathrm{p}\mathrm{a}\mathrm{y}\mathrm{l}\mathrm{o}\mathrm{a}\mathrm{d},\\ & \sum _{i=0}^n\mathbf{v}{\mathbf{x}}_i^{⊺}\le \vartheta ,\\ & \text{𝟙}\cdot {\mathbf{x}}_i^{⊺}=1,\mathrm{\forall }i=0,\dots ,n,\\ & 2{\tilde{y}}_i{y}_i-z_{y_i^2}\le {\tilde{y}}_i^2,\mathrm{\forall }i=0,\dots ,n,\\ & {\tilde{x}}_i{y}_i+{\tilde{y}}_i{x}_i-z_{x_i{y}_i}\le {\tilde{x}}_i{\tilde{y}}_i,\mathrm{\forall }i=0,\dots ,n,\\ \mathrm{v}\mathrm{a}\mathrm{r}.& {\mathbf{x}}_i\in \{0,1{\}}^{\vartheta +1},\mathrm{\forall }i=0,\dots ,n,\\ & z_{y_i^2}\ge 0,\mathrm{\forall }i=0,\dots ,n,\\ & z_{x_i{y}_i}\ge 0,\mathrm{\forall }i=0,\dots ,n,\\ \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}& x_i=\mathbf{v}{\mathbf{x}}_i^{⊺}\mathrm{\&}{y}_i=\sum _{j=0}^{i-1}\mathbf{v}{\mathbf{x}}_j^{⊺}.\end{array}$

5. Simulation

We carry out experimental analysis on the optimality of the proposed method benchmarked against the brute-force method. The experimental setup is described as follows. For the predictive model, we use the residual dense network (RDN), which has its origins in low-level computer vision tasks such as super-resolution imaging (Y. Zhang et al., 2018) and image restoration (Y. Zhang et al., 2021). This neural network model is characterised by a tangled labyrinth of residual and dense connections. It is trained on the BOSSbase dataset (Bas et al., 2011), which originated from an academic competition for digital steganography, and comprises a large collection of greyscale photographs covering a wide variety of subjects and scenes. The algorithms are tested on selected images from the USC-SIPI dataset (Weber, 2006). All the images are resized to a resolution of $256\times 256$ pixels via Lanczos resampling (Duchon, 1979). The border pixels along with half of the rest of the pixels are designated as the context. Accordingly, the number of query pixels equals $(254\times 254)/2$. We display both distortion and capacity as divided by the number of query pixels.

Figure 5 shows the absolute error distribution for each test image. It is observed that most of the error values are below around 30 to 50, depending on the image. We set n = 55 conservatively in the sense that nearly every value of non-zero occurrence is included. We implement the algorithms with respect to different quota settings ($\vartheta =1,2,3,4$). Figures 6–9 show performance evaluations of the proposed optimisation algorithm. Each point of the curve indicates the minimum distortion of a solution under a specific capacity constraint. In the vast majority of cases, the solutions found by the proposed method are identical to those given by the brute-force method. When failing to find the optimal solutions, the objective values reached are within a small distance from the optimal ones. Hence, even though optimal solutions cannot always be guaranteed, the results suggest that the proposed method can attain near-optimal performance.

Figure 5. Absolute error histograms with highlighted empty bins. (a) Aeroplane. (b) Lena. (c) Mandrill and (d) Peppers.

Figure 6. Payload–distortion curves for optimality analysis against brute-force search with highlighted discrepancies ($\vartheta =1$).

Figure 7. Payload–distortion curves for optimality analysis against brute-force search with highlighted discrepancies ($\vartheta =2$).

Figure 8. Payload–distortion curves for optimality analysis against brute-force search with highlighted discrepancies ($\vartheta =3$).

Figure 9. Payload–distortion curves for optimality analysis against brute-force search with highlighted discrepancies ($\vartheta =4$).

6. Conclusion

This paper studies a mathematical optimisation problem applied to reversible steganography. We formulate automatic coding in prediction error modulation as a nonlinear discrete optimisation problem. The objective is to minimise distortion under a constraint on capacity. We discuss the complexity of a brute-force search algorithm and the linearisation techniques for the logarithmic capacity constraint and the quadratic distortion objective. The problem is transformed into an iterative mixed-integer linear programming problem with binary-integer variables and slack variables. Our simulation results validate the near-optimality of the proposed algorithm.

Disclosure statement

No potential conflict of interest is reported by the author.