Efficient Image Splicing Localization via Contrastive Feature Extraction

TL;DR

This paper introduces cPCA++, a fast and efficient contrastive feature extraction technique, applied to image splicing localization, achieving comparable accuracy to deep learning methods without extensive training.

Contribution

The paper presents cPCA++, a novel contrastive PCA variant that is more efficient than cPCA and effective for image splicing localization without iterative training.

Findings

cPCA++ outperforms traditional PCA in highlighting dataset-specific structures.

The method is significantly faster than cPCA and deep learning approaches.

cPCA++ achieves comparable localization accuracy to state-of-the-art neural networks.

Abstract

In this work, we propose a new data visualization and clustering technique for discovering discriminative structures in high-dimensional data. This technique, referred to as cPCA++, utilizes the fact that the interesting features of a "target" dataset may be obscured by high variance components during traditional PCA. By analyzing what is referred to as a "background" dataset (i.e., one that exhibits the high variance principal components but not the interesting structures), our technique is capable of efficiently highlighting the structure that is unique to the "target" dataset. Similar to another recently proposed algorithm called "contrastive PCA" (cPCA), the proposed cPCA++ method identifies important dataset specific patterns that are not detected by traditional PCA in a wide variety of settings. However, the proposed cPCA++ method is significantly more efficient than cPCA, because it does not require the parameter sweep in the latter approach. We applied the cPCA++ method to the problem of image splicing localization. In this application, we utilize authentic edges as the background dataset and the spliced edges as the target dataset. The proposed method is significantly more efficient than state-of-the-art methods, as the former does not require iterative updates of filter weights via stochastic gradient descent and backpropagation, nor the training of a classifier. Furthermore, the cPCA++ method is shown to provide performance scores comparable to the state-of-the-art Multi-task Fully Convolutional Network (MFCN).

Figures

Tables

Table 1. TABLE I : Overview of the datasets used for comparing the PCA, cPCA, t-SNE, and cPCA++ dimensionality reduction methods. Note that N f subscript 𝑁 𝑓 N_{f} denotes the number of foreground samples, N b subscript 𝑁 𝑏 N_{b} denotes the number of background samples, and M 𝑀 M denotes the original feature dimension.

Example	$N_f$	$N_b$	$M$
Synthetic	$400$	$400$	$30$
MNIST over Grass	$5000$	$5000$	$784$
Mice Protein Expression	$270$	$135$	$77$
MHealth Measurements	$6451$	$3072$	$23$
Single Cell RNA-Seq of Leukemia Patient	$7898$	$1985$	$500$

Table 2. TABLE II : Time required for the different algorithms to perform the required dimensionality reduction for the various datasets studied in Sec.  III-B . All times listed in the table are in seconds. Boldface is used to indicate shortest runtimes and average cPCA++ speedup.

Example	cPCA	t-SNE	cPCA++
Synthetic	$0.062$	$6.62$	$0.0017$
MNIST over Grass	$26$	$811$	$1.00$
Mice Protein Expression	$0.13$	$2.23$	$0.0033$
MHealth Measurements	$0.091$	$1388$	$0.00065$
RNA-Seq of Leukemia Patient	$11.70$	$2155$	$0.86$
Average cPCA++ Speedup	$\mathrm{𝟓𝟏}$x	$\mathrm{𝟒𝟐𝟖𝟔𝟓𝟒}$x	$\mathrm{𝟏}$x

Table 3. TABLE III : Edge-based M ​ C ​ C 𝑀 𝐶 𝐶 MCC Scores on Columbia and Nimble WEB Datasets. Boldface is used to emphasize best performance.

Dataset	cPCA++	MFCN
Columbia	$0.385$	$0.329$
Nimble WEB	$0.388$	$0.297$

Table 4. TABLE IV : Edge-based F 1 subscript 𝐹 1 F_{1} Scores on Columbia and Nimble WEB Datasets. Boldface is used to emphasize best performance.

Dataset	cPCA++	MFCN
Columbia	$0.359$	$0.312$
Nimble WEB	$0.376$	$0.273$

Equations

$$\hat{\mathbf{\Sigma }}\triangleq \frac{1}{N}\mathbf{Z}{\mathbf{Z}}^{\mathsf{T}}$$

$${\mathbf{F}}_{\mathrm{PCA}}\triangleq {\mathrm{evecs}}_K\left(\hat{\mathbf{\Sigma }}\right)$$

$${\mathbf{Y}}_{\mathrm{PCA}}={\mathbf{F}}_{\mathrm{PCA}}^{\mathsf{T}}\mathbf{Z}.$$

$${\mathbf{Z}}_b\sim \mathcal{N}(0_{M\times N_b},{\mathbf{\Sigma }}_b),[{\mathbf{Z}}_b\text{ from background dataset}]$$

$$p(\mathbf{Z};\mathbf{A},\mathbf{\Sigma })\triangleq \frac{\exp \left(-\frac{1}{2}\mathsf{Tr}[{\mathbf{\Sigma }}^{-1}(\mathbf{Z}-\mathbf{A})(\mathbf{Z}-\mathbf{A}{)}^{\mathsf{T}}]\right)}{\sqrt{(2\pi {)}^{MN}|\mathbf{\Sigma }{|}^N}}.$$

$${\mathbf{Z}}_f\sim \mathcal{N}(\mathbf{W}{\mathbf{Y}}_f,{\mathbf{\Sigma }}_f),[{\mathbf{Z}}_f\text{ from foreground dataset}]$$

$$\displaystyle\boldsymbol{Z}\triangleq\left[\begin{array}[]{cc}\boldsymbol{Z}_{f}&\boldsymbol{Z}_{b}\end{array}\right]$$

$$\mathbb{E}\left[\mathbf{Z}\right]$$

$$\displaystyle=\left[\begin{array}[]{cc}\boldsymbol{W}\boldsymbol{Y}_{\!\!f}&\boldsymbol{0}_{M\times N_{b}}\end{array}\right]$$

$$=\mathbf{W}{\mathbf{Y}}_f\mathbf{T}$$

$$\displaystyle\boldsymbol{T}\triangleq\left[\begin{array}[]{cc}\boldsymbol{I}_{N_{f}\times N_{f}}&\boldsymbol{0}_{N_{f}\times N_{b}}\end{array}\right].$$

$$\mathbf{Z}\sim \{\begin{array}{ll}\mathcal{N}(0_{M\times N},{\mathbf{\Sigma }}_0),& {\mathcal{H}}_0\\ \mathcal{N}(\mathbf{W}{\mathbf{Y}}_f\mathbf{T},{\mathbf{\Sigma }}_1),& {\mathcal{H}}_1\end{array}.$$

$${\mathbf{Z}}_b$$

$${\mathbf{Z}}_f$$

$$p({\mathcal{H}}_1|\mathbf{Z})$$

$$p(\mathbf{Z}|{\mathcal{H}}_1)p({\mathcal{H}}_1)$$

$$\frac{p(\mathbf{Z}|{\mathcal{H}}_1)}{p(\mathbf{Z}|{\mathcal{H}}_0)}\ge \gamma \triangleq \frac{p({\mathcal{H}}_0)}{p({\mathcal{H}}_1)}.$$

$$\mathrm{GLRT}(\mathbf{W})\triangleq \frac{\underset{{\mathbf{\Sigma }}_1,{\mathbf{Y}}_f}{\max }{f}_1(\mathbf{Z};{\mathbf{\Sigma }}_1,{\mathbf{Y}}_f)}{\underset{{\mathbf{\Sigma }}_0}{\max }{f}_0(\mathbf{Z};{\mathbf{\Sigma }}_0)}\ge \gamma$$

$$f_1(\mathbf{Z};{\mathbf{\Sigma }}_1,{\mathbf{Y}}_f)$$

$$f_0(\mathbf{Z};{\mathbf{\Sigma }}_0)$$

$$\underset{{\mathbf{\Sigma }}_0}{\max }{f}_0(\mathbf{Z};{\mathbf{\Sigma }}_0)$$

$$={\left((2\pi e{)}^M\left|\frac{1}{N}\mathbf{Z}{\mathbf{Z}}^{\mathsf{T}}\right|\right)}^{-N/2}$$

$$\underset{{\mathbf{\Sigma }}_1}{\max }$$

$$=f_1\left(\mathbf{Z};\frac{1}{N}(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T})(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T}{)}^{\mathsf{T}},{\mathbf{Y}}_f\right)$$

$$={\left((2\pi e{)}^M\left|\frac{1}{N}(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T})(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T}{)}^{\mathsf{T}}\right|\right)}^{-N/2}.$$

$$\mathrm{GLRT}$$

$$=\frac{\underset{{\mathbf{Y}}_f}{\max }{\left((2\pi e{)}^M\left|\frac{1}{N}(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T})(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T}{)}^{\mathsf{T}}\right|\right)}^{-\frac{N}{2}}}{{\left((2\pi e{)}^M\left|\frac{1}{N}\mathbf{Z}{\mathbf{Z}}^{\mathsf{T}}\right|\right)}^{-\frac{N}{2}}}$$

$$\stackrel{(a)}{⟺}\frac{\left|\mathbf{Z}{\mathbf{Z}}^{\mathsf{T}}\right|}{\underset{{\mathbf{Y}}_f}{\min }\left|(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T})(\mathbf{Z}-\mathbf{W}{\mathbf{Y}}_f\mathbf{T}{)}^{\mathsf{T}}\right|}$$

$$\mathrm{GLRT}(\mathbf{W})$$

$$=\frac{\left|{\mathbf{W}}^{\mathsf{T}}{\mathbf{R}}_b^{-1}\mathbf{W}\right|}{\left|{\mathbf{W}}^{\mathsf{T}}(N_f/N_b\cdot {\mathbf{R}}_f+{\mathbf{R}}_b{)}^{-1}\mathbf{W}\right|}$$

$${\mathbf{R}}_b\triangleq \frac{1}{N_b}{\mathbf{Z}}_b{\mathbf{Z}}_b^{\mathsf{T}},{\mathbf{R}}_f\triangleq \frac{1}{N_f}{\mathbf{Z}}_f{\mathbf{Z}}_f^{\mathsf{T}}$$

$${\mathbf{Y}}_f={\mathbf{F}}^{\mathsf{T}}{\mathbf{Z}}_f$$

$$\mathbf{F}={\mathbf{R}}_b^{-1}\mathbf{W}({\mathbf{W}}^{\mathsf{T}}{\mathbf{R}}_b^{-1}\mathbf{W}{)}^{-1}.$$

$$\mathbf{Q}\triangleq {\mathbf{R}}_b^{-1}{\mathbf{R}}_f.$$

$${\mathbf{R}}_b$$

Full text

I Introduction

With the proliferation of low-cost and highly sophisticated image editing software, the generation of realistic-looking manipulated images has never been more easily accessible. At the same time, the prevalence of social media applications, such as Twitter, has made it very easy to quickly circulate these manipulated or fake images. Thus, there has been an increasing interest in developing forensic techniques to detect and localize forgeries (also referred to as manipulations or attacks) in images. One of the most common types of forgery is the image splicing attack. A splicing attack is a forgery in which a region from one image (i.e., the donor image) is copied and pasted onto another image (i.e., the host image). Forgers often use splicing to give a false impression that there is an additional object present in the image, or to remove an object from the image. Image splicing can be potentially used in generating false propaganda for political purposes. For example, the altered or manipulated image in the top row of Figure 1 shows Bill Clinton shaking hands with Saddam Hussein, although this event never occurred. The authentic images used to create the fake image can be seen in the bottom row of Figure 1.

An additional splicing example111https://www.nist.gov/itl/iad/mig/media-forensics-challenge is shown in Figure 2. The top row shows (from left to right): the manipulated image (also referred to as the probe or spliced image), the host image, and the donor image. In this example, the murky water from the donor image was copied and pasted on top of the clear water in the host image. This gives the false appearance that the water in the pool is not clean. The murky water in the manipulated image is referred to as the spliced surface or region. The bottom row in Figure 2 shows (from left to right): the spliced surface (or region) ground truth mask and the spliced boundary (or edge) ground truth mask. These two types of masks provide different ways of highlighting the splicing manipulation. The surface ground truth mask is a per-pixel binary mask which specifies whether each pixel in the given manipulated image is part of the spliced surface (or region). We use the color white to denote a pixel belonging to the spliced surface and the color black to denote a pixel not belonging to the spliced surface. The edge ground truth mask is a per-pixel binary mask which specifies whether each pixel in the given probe image is part of the boundary of the spliced surface. In this case, we use the color white to denote a pixel belonging to the spliced boundary and the color black to denote a pixel not belonging to the spliced boundary.

There are two main problems in the splicing literature: detection and localization. The detection problem refers to the problem of classifying an image as a whole as either spliced or authentic, without localizing the spliced region or boundary. Many of the current techniques only address the detection problem, and do not address the localization problem. However, recently, more effort has been made in the forensics community in addressing the localization problem, which aims to classify each pixel in an image as either belonging to the spliced surface/boundary or not. Most of the recent techniques utilize deep-learning-based techniques such as convolutional neural networks (CNNs) to localize splicing attacks. Although CNNs have yielded some improvement in the localization of splicing attacks, they rely on careful selection of hyperparameters, network architecture, and initial filter weights. Furthermore, CNNs require a very long training time because they utilize stochastic gradient descent and backpropagation to iteratively update the filter or kernel weights. It is desired to have an approach that relies more on fundamental mathematical theory, does not require a significant amount of experimental tuning, and does not require long training and testing times.

In this paper, we propose such an approach, the basis of which is Principal Component Analysis (PCA). As we will see, traditional PCA is not adequate for the task of image splicing localization, and so we will first present a new dimensionality-reduction technique for performing discriminative feature extraction, referred to as cPCA++, which is inspired by a recently proposed algorithm called “contrastive PCA” (cPCA) [1]. We then propose a new approach for image splicing localization based on cPCA++, which is significantly more efficient than state-of-the-art techniques such as the Multi-task Fully Convolutional Network (MFCN) [2], and still achieves comparable performance scores. Also, cPCA++ is derived via matrix factorization, which allows us to identify underlying bases and corresponding weightings in the data that can be used for denoising applications (see App. B for such an example).

The rest of the paper is organized as follows. Sec. II reviews existing splicing localization techniques. Sec. III presents cPCA++ as a new dimensionality-reduction technique, discusses how it differs from traditional PCA, and presents simulation results. The general framework or pipeline of the proposed cPCA++ approach for splicing localization is presented in Sec. IV. Experimental results on splicing datasets are presented in Sec. V. Finally, concluding remarks are given in Sec. VI.

II Related Work

Zampoglou et al. [3] conducted a comprehensive review of non-deep-learning-based techniques for image splicing localization, and provided a comparison of their performance. These techniques can be roughly grouped into the following three categories based on the type of feature or artifact they exploit: noise patterns [4, 5, 6, 7], Color Filter Array (CFA) interpolation patterns [8, 9], and JPEG-related traces [10, 11, 12, 13, 14, 15, 16, 17, 18]. Recently, there has been an increasing interest in the application of deep-learning-based techniques to general image forensics and splicing detection/localization [2, 19, 20, 21, 22, 23, 24, 25, 26]. Specifically, convolutional neural networks (CNNs) have attracted a significant amount of attention in the forensics community, due to the promising results they have yielded on a variety of image-based tasks such as object recognition and semantic segmentation [27, 28, 29]. One of the state-of-the-art CNN-based techniques for image splicing localization is the Multi-task Fully Convolutional Network (MFCN), proposed in [2]. In this work, one branch is used to learn the spliced surface ground truth mask (or label), while the other branch is used to learn the spliced boundary/edge ground truth mask. In addition to the surface ground truth mask, the boundaries between inserted regions and their host background can be an important indicator of a manipulated area. It is shown in [2] that simultaneously training on the surface and edge labels can yield finer localization of the splicing manipulation, as compared to training only on the surface label. It is also important to note that the MFCN-based method and the proposed cPCA++ method both output an edge-based probability map (as we will see in Sec. IV), which allows for a direct comparison.

Although CNN-based approaches have yielded promising results in the field of image forensics, they rely on careful selection of hyperparameters, network architecture, and initial filter weights. The behavior of the resulting CNN classifier, after learning of the coupled feature extractor and classifier blocks, is difficult to explain mathematically. Furthermore, CNNs require a long training time since the filter weights need to be iteratively updated via stochastic gradient descent and backpropagation. It is desired to have an approach that relies more on fundamental mathematical theory, does not require a significant amount of experimental tuning, and does not require long training/testing times.

In this paper, we propose such an approach, based on a new version of Principal Component Analysis (PCA). In the context of image splicing localization, the two classes of interest (i.e., spliced and authentic boundaries) are very similar in terms of their covariance matrices, and traditional PCA is not able to effectively discriminate between the two classes. Instead, we propose a new version of PCA, referred to as cPCA++, that is able to perform discriminative feature extraction when dealing with extremely similar classes. We then propose a new approach for image splicing localization based on cPCA++. The proposed approach is mathematically tractable and does not require a significant amount of experimental tuning. Unlike CNNs, the proposed approach does not require iterative updates of filter weights via stochastic gradient descent and backpropagation, and thus is much more efficient than CNN-based approaches. In addition, we will see that the cPCA++ approach does not require the training of a classifier (e.g., support vector machines or random forests), which greatly speeds up the method. Also, the proposed approach can be readily parallelized due to its lack of dependence on inherently serial methods such as stochastic gradient descent. In the next section, we will derive the cPCA++ method, and then apply it to the image splicing localization problem in Sec. IV.

III cPCA++

III-A The cPCA++ Method

As we will see in Sec. IV, one characteristic of the image splicing localization problem is the strong similarity between authentic and spliced edges, especially when spliced edges have been smoothed over with a low-pass filter. Thus, traditional Principal Component Analysis (PCA) is not able to effectively discriminate between spliced and authentic edges. In order to deal with this issue, we will first study the problem of dimensionality reduction for extremely similar classes. Doing so will allow us to develop a new and efficient dimensionality-reduction technique that we refer to as cPCA++, which is inspired by a recently proposed algorithm called “contrastive PCA” (cPCA) [1]. In order to derive the cPCA++ method, we shall first consider the traditional PCA method.

Suppose we are presented with a data matrix $\boldsymbol{Z}\in\mathbb{R}^{M\times N}$, where $M$ denotes the original feature dimension of the data and $N$ denotes the number of instances included in $\boldsymbol{Z}$. PCA first computes the empirical covariance matrix of $\boldsymbol{Z}$:

[TABLE]

where we assumed that $\boldsymbol{Z}$ has zero-mean. Next, PCA would compute the subspace spanned by the $K$ top or leading eigenvectors of $\widehat{\boldsymbol{\Sigma}}$ (i.e., those corresponding to the largest $K$ eigenvalues). The basis for this space would constitute the filters used to process the input data:

[TABLE]

where $K$ denotes the number of principal eigenvectors to return. Now, a low-dimensional version of the input data $\boldsymbol{Z}$ can be obtained as:

[TABLE]

It is known that the PCA filters $\boldsymbol{F}_{\mathrm{PCA}}$ preserve the most energy in $\boldsymbol{Z}$ after the transformation. Observe that this property may not yield separability of classes in $\boldsymbol{Z}$. That is, if the data matrix $\boldsymbol{Z}$ is composed of multiple classes (e.g., spliced and authentic edges in our case), performing traditional PCA will not necessarily allow us to find a representation in which we can separate the classes.

In this work, we will take a different approach, inspired by a recently proposed algorithm called the “contrastive PCA” (cPCA) method [1], which obtains discriminative filters. This approach focuses on finding directions that yield large variations for one dataset, referred to as the “target” or “foreground” dataset, while simultaneously yielding small variations for another dataset, referred to as the “background” dataset. There are two problem setups that benefit from this scheme, and we examine them in the following remark.

Remark 1**.**

In [1], the problem is to discover structures that are unique to one dataset relative to another. There are two ways of utilizing the solution to this problem. They are described as follows.

1)

Consider the example described more thoroughly in Sec. III-B2. We are provided a “target” or “foreground” dataset with digits superimposed atop of grass background images. Observe that the target dataset contains multiple classes (i.e., the two digits, [math] and $1$). Examples of instances of this dataset are illustrated in Fig. 4. We wish to learn filters that are tuned to the digits, as opposed to the relatively strong grass background images. In order to accomplish this task, consider having access to a “background” dataset that only contains instances of grass images (i.e., without digits). The task is to discover the structures that are unique to the target dataset and not present in the background dataset. If this is accomplished, the filters are hoped to be able to differentiate between the two digits in the “target” dataset.

2)

In this alternative setup, we would like to solve a binary classification problem with the caveat that the two classes are very similar, and are dominated by the structures that are present in both classes. This is the case in the image splicing localization problem, where one class is the spliced edges and the other class is the authentic edges. In this scenario, we consider one class (i.e., the spliced edges) to be the target dataset and the other class (i.e., the authentic edges) to be the background dataset. If we are able to learn filters that are tuned to the structures unique to the spliced edges, then these filters form a basis for an efficient classifier that is able to differentiate between spliced and authentic edges. We will see the form of this classifier in Sec. IV.

We will take a detection-based approach to dimensionality reduction. We will setup a detector whose task is to identify whether or not a presented data matrix contains the interesting or special structures present in the foreground dataset (observe that the data matrix may contain background instances as well). By following this route, we will see that the detector will tune itself to specifically look for these special structures in the presented data matrix. Fortunately, the detector will look for these structures by transforming the incoming data matrix with a set of linear filters. This means that we may simply utilize the output of these filters as a low-dimensionality representation of the input data matrix, tuned specifically to the interesting structures.

To setup the problem mathematically, we assume that the data matrix $\boldsymbol{Z}_{b}\in\mathbb{R}^{M\times N_{b}}$ is collected in an independently and identically distributed (i.i.d.) manner from the background dataset and has the following distribution:

[TABLE]

where $\boldsymbol{\Sigma}_{b}\in\mathbb{R}^{M\times M}$ denotes an unknown covariance matrix of the background dataset and $\boldsymbol{Z}\sim\mathcal{N}(\boldsymbol{A},\boldsymbol{\Sigma})$ with $\boldsymbol{Z},\boldsymbol{A}\in\mathbb{R}^{M\times N}$ and positive-definite $\boldsymbol{\Sigma}\in\mathbb{R}^{M\times M}$ means that $\boldsymbol{Z}$ has the following probability density function (pdf):

[TABLE]

On the other hand, we assume that when the instances of $\boldsymbol{Z}_{f}\in\mathbb{R}^{M\times N_{f}}$ are sampled independently from the foreground dataset, $\boldsymbol{Z}_{f}$ has the following distribution:

[TABLE]

where the mean of $\boldsymbol{Z}_{f}$ can be factored as the product of a basis matrix $\boldsymbol{W}\in\mathbb{R}^{M\times K}$ and a low-dimensionality representation matrix $\boldsymbol{Y}_{\!\!f}\in\mathbb{R}^{K\times N_{f}}$. The inner-dimension $K$ represents the underlying rank of the data matrix $\boldsymbol{Z}_{f}$, which is generally assumed to be much smaller than the initial dimension $M$.

Now, it is assumed that we are presented with a data matrix $\boldsymbol{Z}\in\mathbb{R}^{M\times N}$ that could be drawn completely from the background dataset (we refer to this case as the null hypothesis $\mathcal{H}_{0}$) or contain instances from the background dataset and some instances from the foreground dataset (we refer to this case as hypothesis $\mathcal{H}_{1}$). That is, the data matrix $\boldsymbol{Z}$ under $\mathcal{H}_{1}$ can be written as:

[TABLE]

where $\boldsymbol{Z}_{f}\in\mathbb{R}^{M\times N_{f}}$ and $\boldsymbol{Z}_{b}\in\mathbb{R}^{M\times N_{b}}$ denote data instances from the foreground and background datasets, respectively. When $\mathcal{H}_{1}$ is active, we assume that we know the value of $N_{f}$. Given (4) and (6) above, we have that the mean of $\boldsymbol{Z}$ under $\mathcal{H}_{1}$ can be written as:

[TABLE]

where we have introduced the matrix $\boldsymbol{T}\in\mathbb{R}^{N_{f}\times N}$, where $N\triangleq N_{f}+N_{b}\geq 2$, with the following structure:

[TABLE]

Therefore, the detection/classification problem we wish to examine is the determination of the active hypothesis (either $\mathcal{H}_{0}$ or $\mathcal{H}_{1}$), where the data matrix is described, under each hypothesis, by:

[TABLE]

Note that $\boldsymbol{\Sigma}_{0}$ and $\boldsymbol{\Sigma}_{1}$ denote the covariance matrix under $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$, respectively. The underlying assumption is that when the data matrix contains only instances from the background dataset, then the mean would be zero. We thus assume that all data instances of $\boldsymbol{Z}$ have the mean of the respective partitions subtracted off. For example, if the raw (unprocessed) data matrix is given by \widetilde{\boldsymbol{Z}}\triangleq\left[\begin{array}[]{cc}\widetilde{\boldsymbol{Z}}_{f}&\widetilde{\boldsymbol{Z}}_{b}\end{array}\right], then the partitions of the data matrix \boldsymbol{Z}=\left[\begin{array}[]{cc}\boldsymbol{Z}_{f}&\boldsymbol{Z}_{b}\end{array}\right] are given by:

[TABLE]

In solving this detection problem, we will be able to obtain the matrix $\boldsymbol{Y}_{\!\!f}\in\mathbb{R}^{K\times N_{f}}$, which is the desired low-dimensionality representation of the foreground dataset. While it may initially appear that the detection problem is relatively simple—especially when the mean under the foreground dataset is relatively large in magnitude compared to the power in the covariance matrices, the detection problem is difficult when the mean is small compared to the covariance power (i.e., it is masked by large variation that naturally occurs in both the background and foreground datasets), as outlined earlier and in [1]. This is true in the context of image splicing localization, so it is assumed that the mean under the foreground dataset is small compared to the covariance power.

Now, we would like to determine if a given data matrix $\boldsymbol{Z}$ belongs to the distribution in the null hypothesis ($\mathcal{H}_{0}$) or the alternative hypothesis ($\mathcal{H}_{1}$). A classical approach to such problems is the likelihood ratio test, which is derived via the test:

[TABLE]

where step $(a)$ is due to Bayes’ rule, and $p(\mathcal{H}_{0})$ and $p(\mathcal{H}_{1})$ denote the prior probabilities of hypotheses $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$, respectively. The above can readily be re-arranged as:

[TABLE]

This is referred to as the “likelihood ratio test.” It evaluates the likelihood ratio of the data matrix $\boldsymbol{Z}$ under each hypothesis and compares it to a threshold given by $\gamma$. The likelihood ratio on the left-hand-side of (20) is desired to have a large value when the hypothesis $\mathcal{H}_{1}$ is active, but is desired to take on a small value when the hypothesis $\mathcal{H}_{0}$ is active. The threshold on the right-hand-side of (20) is generally swept to achieve different performance that trades the probability of detection to false-alarm. Observe that the likelihood ratio test is only useful when the probabilities $p(\boldsymbol{Z}|\mathcal{H}_{1})$ and $p(\boldsymbol{Z}|\mathcal{H}_{0})$ can be computed (i.e., the probability density function is completely known or all of its parameters can be estimated). In our case, however, we do not have knowledge of many parameters, including the mean components $\boldsymbol{W}$ and $\boldsymbol{Y}_{\!\!f}$, as well as the covariance matrices $\boldsymbol{\Sigma}_{0}$ and $\boldsymbol{\Sigma}_{1}$. In order to deal with this situation, a common approach is to utilize what is referred to as the generalized likelihood ratio test (GLRT), which optimizes over parameters that are unknown from within the ratio. In our case, the GLRT is given by [30, p. 200]:

[TABLE]

where $\gamma$ denotes the threshold parameter, $f_{1}(\boldsymbol{Z};\boldsymbol{\Sigma}_{1},\boldsymbol{Y}_{\!\!f})$ denotes the Gaussian pdf of $\boldsymbol{Z}$ under the hypothesis $\mathcal{H}_{1}$, and $f_{0}(\boldsymbol{Z};\boldsymbol{\Sigma}_{0})$ denotes the Gaussian pdf of $\boldsymbol{Z}$ under the null-hypothesis $\mathcal{H}_{0}$. That is,

[TABLE]

In (21), we clearly indicate that the resulting GLRT statistic will be a function of the matrix $\boldsymbol{W}$, which is the basis for the feature mean under $\mathcal{H}_{1}$. We will return to this point later in the derivation. Now, we can optimize each pdf over the unknown covariance matrix, independently, as prescribed by the GLRT in (21) to obtain that the optimal covariance matrices for each case is given by the empirical covariance matrix (i.e., the maximum-likelihood covariance matrix estimator for a Gaussian distribution [31, p. 89]), so that:

[TABLE]

and the pdf for $\mathcal{H}_{1}$ is given by

[TABLE]

Thus, the GLRT simplifies to:

[TABLE]

where step $(a)$ states that the GLRT statistic is equivalent to (26). Now, using the derivation in [32], we have that the GLRT statistic can be written as (see equation (8) and related derivations in App. A from [32]):

[TABLE]

We will now optimize over our free variable $\boldsymbol{W}$. When the null-hypothesis is active, it can be shown that the GLRT statistic is independent of $\boldsymbol{W}$. On the other hand, when $\mathcal{H}_{1}$ is active, we have that the GLRT is dependent on $\boldsymbol{W}$ and is given by:

[TABLE]

where $(a)$ is obtained via (8) and we have defined the following quantities:

[TABLE]

where $\boldsymbol{R}_{b}$ denotes the second order statistic associated with the background dataset only, while $\boldsymbol{R}_{f}$ denotes the second order statistic associated with the foreground dataset only. Now, to ensure that the gap in the GLRT statistic between the two hypotheses is as large as possible, we must maximize the GLRT value over $\boldsymbol{W}$ when the hypothesis $\mathcal{H}_{1}$ is active (since the GLRT value is independent of $\boldsymbol{W}$ when the null-hypothesis is active). It turns out that the value of $\boldsymbol{W}$ that maximizes the GLRT is obtained by solving a generalized eigenvalue problem [33] [34, pp. 454–455]. The optimal choice for $\boldsymbol{W}$ turns out to be the principal eigenvectors of the matrix $\boldsymbol{I}_{M\times M}+N_{f}/N_{b}\cdot\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$, or, equivalently222Assuming a square matrix $\boldsymbol{C}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{U}^{-1}$ is diagonalizable with non-negative eigenvalues, the matrix $\boldsymbol{I}+\alpha\boldsymbol{C}$ for $\alpha>0$ can be written as $\boldsymbol{I}+\alpha\boldsymbol{C}=\boldsymbol{U}\boldsymbol{U}^{-1}+\alpha\boldsymbol{U}\boldsymbol{D}\boldsymbol{U}^{-1}=\boldsymbol{U}(\boldsymbol{I}+\alpha\boldsymbol{D})\boldsymbol{U}^{-1}$, which modifies all eigenvalues of $\boldsymbol{C}$ by adding one to a positive-scaled version of them. This does not change the order of the eigenvalues of $\boldsymbol{C}$ since they are assumed to be non-negative and thus $\boldsymbol{I}+\alpha\boldsymbol{C}$ and $\boldsymbol{C}$ share the same principal eigenvectors when $\alpha>0$ and the eigenvalues of $\boldsymbol{C}$ are non-negative., the principal eigenvectors of the matrix $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$ when the eigenvalues of $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$ are non-negative (we will show later that this is indeed the case333We will later show that the eigenvalues of $\boldsymbol{R}_{b}^{-1}\boldsymbol{R}_{f}$ are non-negative. Observe, however, that $\boldsymbol{R}_{b}^{-1}\boldsymbol{R}_{f}$ and $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$ are similar since $\boldsymbol{R}_{b}$ is non-singular [35, p. 53], which means that the eigenvalues of $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$ are also non-negative.). While the solution to $\boldsymbol{W}$ itself is not directly interesting for the problem of image splicing localization (it is relevant for applications related to matrix factorization, denoising, and dictionary learning [36]—which we explore in App. B), the reduced-dimensionality matrix $\boldsymbol{Y}_{\!\!f}$ is much more relevant for dimensionality reduction purposes.

The maximum-likelihood estimator for $\boldsymbol{Y}_{\!\!f}$ (which was used to obtain (27)) is written as

[TABLE]

where $\boldsymbol{F}\in\mathbb{R}^{M\times K}$ is a bank of filters yet to be determined, but has the following dependence on $\boldsymbol{W}$ [32]:

[TABLE]

In fact, the filters in (31) can be shown to be the principal eigenvectors of the matrix:

[TABLE]

This can be seen by noting that the optimal $\boldsymbol{W}$ (being the eigenvectors of $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$) must satisfy $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}\boldsymbol{W}=\boldsymbol{W}\boldsymbol{D}$ for some diagonal matrix $\boldsymbol{D}$ that contains the eigenvalues. From this, we have that $\boldsymbol{R}_{b}^{-1}\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}\boldsymbol{W}=\boldsymbol{R}_{b}^{-1}\boldsymbol{W}\boldsymbol{D}$. Denoting $\boldsymbol{P}=\boldsymbol{R}_{b}^{-1}\boldsymbol{W}\boldsymbol{\Lambda}$ for some diagonal invertible scaling matrix $\boldsymbol{\Lambda}$ that we are free to choose, we have that $\boldsymbol{P}$ denotes the eigenvectors of the matrix $\boldsymbol{Q}=\boldsymbol{R}_{b}^{-1}\boldsymbol{R}_{f}$ (since $\boldsymbol{Q}\boldsymbol{P}=\boldsymbol{P}\boldsymbol{D}$). Substituting this into (31), we have that $\boldsymbol{F}=\boldsymbol{P}\boldsymbol{\Lambda}^{-1}(\boldsymbol{W}^{\mathsf{T}}\boldsymbol{R}_{b}^{-1}\boldsymbol{W})^{-1}=\boldsymbol{P}$, with the choice of the matrix $\boldsymbol{\Lambda}=(\boldsymbol{W}^{\mathsf{T}}\boldsymbol{R}_{b}^{-1}\boldsymbol{W})^{-1}$. Now, we would like to show that the matrix $\boldsymbol{\Lambda}$ is diagonal. First, observe that $\boldsymbol{\Lambda}=\boldsymbol{\Lambda}\boldsymbol{W}^{\mathsf{T}}\boldsymbol{R}_{b}^{-1}\boldsymbol{W}\boldsymbol{\Lambda}=\boldsymbol{P}^{\mathsf{T}}\boldsymbol{R}_{b}\boldsymbol{P}$. The matrix $\boldsymbol{P}^{\mathsf{T}}\boldsymbol{R}_{b}\boldsymbol{P}$ can be shown to be diagonal since the columns of $\boldsymbol{P}$ are $\boldsymbol{R}_{b}$-orthogonal for distinct eigenvalues and the columns of $\boldsymbol{P}$ that correspond to the same eigenvalue can be chosen to be $\boldsymbol{R}_{b}$-orthogonal [37, p. 345]. Thus, $\boldsymbol{\Lambda}=\boldsymbol{P}^{\mathsf{T}}\boldsymbol{R}_{b}\boldsymbol{P}$ is a diagonal matrix. It can also be shown $\boldsymbol{\Lambda}$ is positive definite (and thus invertible) as long as $\boldsymbol{W}$ has full column rank. Furthermore, the matrix $\boldsymbol{W}$ has full column rank because the matrix $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$ can be shown to be diagonalizable. This is because $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$ can be shown to be similar to $\boldsymbol{R}_{b}^{-1/2}\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1/2}$ [38, p. 81], which is real and symmetric and thus diagonalizable [38, p. 97].

The obtained filters attempt to examine the specific structures that differentiate samples from $\mathcal{H}_{0}$ and $\mathcal{H}_{1}$. This leads to a relatively straightforward and efficient method for dimensionality reduction, which is listed in Alg. 1. Also, note that it is assumed that the matrix $\boldsymbol{R}_{b}$ is positive-definite. If it is not positive-definite, it is customary to apply diagonal-loading (addition of $\rho\boldsymbol{I}_{M\times M}$ to the covariance matrix for a small $\rho>0$) to force the covariance matrix to be positive-definite when it is rank-deficient (although, if enough diverse data samples are available, it is rare that the covariance matrix is rank-deficient).

It is not immediately clear that $\boldsymbol{F}$ is a real matrix, since $\boldsymbol{Q}$ may not be symmetric. However, since the matrix $\boldsymbol{Q}$ is similar to the real and symmetric matrix $\boldsymbol{R}\triangleq\boldsymbol{R}_{b}^{-1/2}\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1/2}$ [38, p. 81] (this is because $\boldsymbol{Q}$ is by definition similar to $\boldsymbol{R}_{b}^{1/2}\boldsymbol{Q}\boldsymbol{R}_{b}^{-1/2}=\boldsymbol{R}_{b}^{1/2}\boldsymbol{R}_{b}^{-1}\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1/2}=\boldsymbol{R}_{b}^{-1/2}\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1/2}=\boldsymbol{R}$) where $\boldsymbol{R}_{b}^{1/2}$ denotes the matrix square root of $\boldsymbol{R}_{b}$, we have that $\boldsymbol{R}$ and $\boldsymbol{Q}$ share a spectrum, which is real due to $\boldsymbol{R}$ being a real symmetric matrix [38, p. 78]. It is important to note that the eigenvalues of $\boldsymbol{Q}$ are also non-negative due to it being similar to $\boldsymbol{R}$, and $\boldsymbol{R}$ being a positive-semi-definite matrix444The definition of a positive semi-definite matrix $\boldsymbol{R}$ is that $\boldsymbol{R}$ is symmetric and satisfies $\boldsymbol{v}^{\mathsf{T}}\boldsymbol{R}\boldsymbol{v}\geq 0$ for all $\boldsymbol{v}\neq\mathbb{0}$. In our case, we have $\boldsymbol{v}^{\mathsf{T}}\boldsymbol{R}\boldsymbol{v}=\boldsymbol{v}^{\mathsf{T}}\boldsymbol{R}_{b}^{-1/2}\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1/2}\boldsymbol{v}=\boldsymbol{y}^{\mathsf{T}}\boldsymbol{R}_{f}\boldsymbol{y}=\frac{1}{N_{f}}\boldsymbol{y}^{\mathsf{T}}\boldsymbol{Z}_{f}\boldsymbol{Z}_{f}^{\mathsf{T}}\boldsymbol{y}=\boldsymbol{x}^{\mathsf{T}}\boldsymbol{x}=\|\boldsymbol{x}\|^{2}\geq 0$, where $\boldsymbol{y}\triangleq\boldsymbol{R}_{b}^{-1/2}\boldsymbol{v}$ and we defined $\boldsymbol{x}\triangleq 1/\sqrt{N_{f}}\cdot\boldsymbol{Z}_{f}^{\mathsf{T}}\boldsymbol{y}$ as a transformation of $\boldsymbol{y}$.. That is, since $\boldsymbol{R}$ is positive-semi-definite, its eigenvalues are non-negative, and since $\boldsymbol{Q}$ is similar to $\boldsymbol{R}$, its eigenvalues are also non-negative. This was a requirement for our earlier derivation to show that the principal eigenvectors of $\boldsymbol{I}_{M\times M}+N_{f}/N_{b}\cdot\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$ are the same as the principal eigenvectors of $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$. We also know that the eigenvectors of the matrix $\boldsymbol{R}$ are real since $\boldsymbol{R}$ is real and symmetric [38, p. 97]. Let us now denote one of the eigenvectors of $\boldsymbol{R}$, associated with eigenvalue $\lambda\in\mathbb{R}$ as $\boldsymbol{x}\in\mathbb{R}^{M\times 1}$. Then, we have that $\boldsymbol{R}\boldsymbol{x}=\lambda\boldsymbol{x}$ and $\boldsymbol{Q}\boldsymbol{v}=\lambda\boldsymbol{v}$ for some eigenvector $\boldsymbol{v}\in\mathbb{C}^{M\times 1}$ of $\boldsymbol{Q}$ (since $\lambda$ is an eigenvalue of $\boldsymbol{Q}$ as we have shown above—we will see, however, that we can always find eigenvectors in $\mathbb{R}^{M\times 1}$, as opposed to $\mathbb{C}^{M\times 1}$). Multiplying both sides on the left by $\boldsymbol{R}_{b}^{-1/2}$, we obtain:

[TABLE]

Then, we immediately obtain that:

[TABLE]

where

[TABLE]

denotes the eigenvector of $\boldsymbol{Q}$, which is real due to $\boldsymbol{x}\in\mathbb{R}^{M\times 1}$ and $\boldsymbol{R}_{b}^{-1/2}\in\mathbb{R}^{M\times M}$. Thus, we see that $\boldsymbol{F}$ is indeed a real matrix. Note that $\boldsymbol{F}$ (which contains the principal $K$ eigenvectors of $\boldsymbol{Q}$ as its columns) may not satisfy $\boldsymbol{F}^{\mathsf{T}}\boldsymbol{F}=\boldsymbol{I}_{K\times K}$ since the eigenvectors of $\boldsymbol{Q}$ are not guaranteed to be orthonormal.

To see the advantage of the cPCA++ method, it is useful to compare it with the “contrastive PCA” (cPCA) algorithm from [1]. The cPCA algorithm from [1] performs an eigen-decomposition of the following matrix:

[TABLE]

where $\alpha$ is a contrast parameter that is swept to achieve different dimensionality reduction output. We observe that the cPCA algorithm is inherently sensitive to the relative scales of the covariance matrices $\boldsymbol{R}_{b}$ and $\boldsymbol{R}_{f}$. For this reason, a parameter sweep for $\alpha$ (which adjusts the scaling between the two covariance matrices) was used in [1]. In other words, multiple eigen-decompositions are required in the cPCA algorithm. On the other hand, the cPCA++ algorithm does not require this parameter sweep, as the relative scale between the two covariance matrices will be resolved when performing the eigen-decomposition of the matrix $\boldsymbol{Q}$. To see this property analytically, consider Example 1.

Example 1** (A needle in a haystack).**

We consider the following $M$-feature example in order to better understand the traditional PCA, cPCA, and cPCA++ methods, and how they react to a small amount of information hidden in a significant amount of “noise.” Let the covariance matrices of the background and target datasets be respectively given by (where we note that $\boldsymbol{R}_{b}$ is positive-definite):

[TABLE]

where $\boldsymbol{a},\boldsymbol{c}\in\mathbb{R}^{M\times 1}$ and it is assumed that $\|\boldsymbol{a}\|_{2}=\|\boldsymbol{c}\|_{2}=1$ and $\boldsymbol{a}^{\mathsf{T}}\boldsymbol{c}=0$, making the vectors $\boldsymbol{a}$ and $\boldsymbol{c}$ orthonormal. It is further assumed that $\epsilon\ll\beta$, which means that the vector that is not common to both background and target has very small power in comparison to what is common between them. It is also assumed that $\beta\cdot\rho/\gamma<\epsilon$ and $\rho\ll\gamma$, which means that the diagonal-loading should not affect the result significantly. It is assumed that $\gamma,\beta,\epsilon,\rho>0$. In this example, we would like to obtain a filter to reduce the dimension to $K=1$ from our original $M$ features. Note that (41) describes a matrix with eigenvalues $\beta$ and $\epsilon$ and associated eigenvectors $\boldsymbol{a}$ and $\boldsymbol{c}$. The reason why we are interested in the above problem is that by examining the covariance structure above, we see that the data is inherently noise-dominated. This noise (i.e., the $\boldsymbol{a}\boldsymbol{a}^{\mathsf{T}}$ term) appears in both the background dataset and the target dataset, but the target dataset may contain some class-specific information (i.e., the $\boldsymbol{c}\boldsymbol{c}^{\mathsf{T}}$ term) that is not present in the background dataset. However, this “interesting” information is dominated by the noise due to the eigenvector $\boldsymbol{a}$. We wish to perform dimensionality reduction that rejects variance due to the eigenvector $\boldsymbol{a}$ but preserves information due to the eigenvector $\boldsymbol{c}$ of the target dataset. The obvious solution is to choose the dimensionality reduction filter to be $\boldsymbol{c}$. We will examine how each algorithm performs.

​​​​​​​PCA

When traditional PCA is executed on the target dataset, the principal eigenvector (and thus $\boldsymbol{F}_{\mathrm{PCA}}$) will be found to be $\boldsymbol{a}$ since $\epsilon\ll\beta$. In addition, since $\rho\ll\gamma$, $\boldsymbol{a}$ is also the principal eigenvector of the background dataset, which is not useful in extracting the interesting structure from the foreground dataset. This is because $\boldsymbol{F}_{\mathrm{PCA}}^{\mathsf{T}}\boldsymbol{c}=0$, which means that this filter will in fact null-out the “interesting” eigenvector.

​​​​​​​cPCA

When cPCA is executed on the background and target datasets, we have that $\boldsymbol{Q}_{\mathrm{cPCA}}=\boldsymbol{R}_{f}-\alpha\boldsymbol{R}_{b}$ will be given by:

[TABLE]

Now, (42) allows us to see exactly the reasoning behind sweeping $\alpha$. When $\alpha=0$, we obtain the PCA method operating on the target dataset (which would pick the filter defined by $\boldsymbol{a}$). When $\alpha\rightarrow\infty$, we obtain the PCA method operating on the background dataset (which still would pick the filter defined by $\boldsymbol{a}$). Instead, the desired value of $\alpha$ is one that nulls out the eigenvector $\boldsymbol{a}$—e.g., $\alpha=\beta/\gamma$. Observe that when this choice of $\alpha$ is made (or chosen through a sweep), we have that the principal eigenvector (and thus $\boldsymbol{F}_{\mathrm{cPCA}}$) of $\boldsymbol{Q}_{\mathrm{cPCA}}$ will be $\boldsymbol{c}$. This is precisely what we want the filter to be in order to differentiate between the two datasets (or extract the useful information from the target dataset relative to the background). Although the optimal $\alpha$ value is evident in this simple analytic example (due to the fact that we know the underlying structure indicated in (40) and (41)), this is not the case in most experiments and thus it is necessary to sweep $\alpha$ over the unbounded range $\alpha\in[0,\infty)$.

​​​​​​​cPCA++

For the analysis of the cPCA++ method, let us first consider the inverse of $\boldsymbol{R}_{b}$, via the matrix inversion lemma:

[TABLE]

where step $(a)$ is due to the assumption that $0<\rho\ll\gamma$ and $\boldsymbol{P}_{\boldsymbol{a}}^{\perp}=\boldsymbol{I}_{M\times M}-\boldsymbol{a}\boldsymbol{a}^{\mathsf{T}}$ denotes the projection onto the left null-space of $\boldsymbol{a}$ **[38, p. 22]** (since $\boldsymbol{P}_{\boldsymbol{a}}^{\perp}\boldsymbol{a}=\mathbb{0}_{M}$). Now, multiplying $\boldsymbol{R}_{b}^{-1}$ by $\boldsymbol{R}_{f}$, we obtain:

[TABLE]

where step $(a)$ is due to the fact that $\boldsymbol{a}$ and $\boldsymbol{c}$ are orthonormal. This means that the $\boldsymbol{Q}$ matrix in the case of the cPCA++ method is approximately a rank-$1$* matrix with principal eigenvector $\boldsymbol{c}$. Clearly then, the filter chosen by cPCA++ will be $\boldsymbol{F}=\boldsymbol{c}$, which is the desired result. Observe that this result was obtained without the need for a hyper-parameter sweep (i.e., like the sweep over the $\alpha$ parameter in cPCA).*

Note that another recent publication [39], also based on the general principle of [1], was brought to our attention. However, careful review of this article shows that the approach in [39] is different from ours. In [39], the algorithm was obtained in a similar manner as the cPCA algorithm (i.e., it was formulated as a slightly modified maximization problem). On the other hand, our proposed cPCA++ approach is fundamentally different because it directly addresses the matrix factorization and dimensionality reduction nature of our problem, and therefore, gives rise to an intuitive per-instance classification algorithm (see Alg. 2 further ahead) that can be applied to image splicing localization, and also addresses image denoising problems (see App. B).

In the next subsection, we will simulate both the cPCA method and the cPCA++ method for a collection of experiments from [1] (we will also simulate the popular t-SNE dimensionality reduction method [40]). We will observe that the cPCA++ algorithm will be able to discover features from within the target dataset without the need to perform a parameter sweep, which will improve the running time of the method significantly (see Sec. III-C).

III-B Performance Comparison of Feature Extraction Methods

In this subsection, we will examine many of the same experiments that the cPCA algorithm555https://github.com/abidlabs/contrastive was applied to in [1]. Table I lists the parameters of the various datasets that will be examined in this section. In all cases, the desired dimensions from the methods will be set to $K=2$ (i.e., after feature reduction). We also incorporated the t-Distributed Stochastic Neighbor Embedding (t-SNE) algorithm [40] in the results666https://lvdmaaten.github.io/tsne/##implementations, even though t-SNE does not provide an explicit function for the discovered feature reduction. For example, after learning the filters $\boldsymbol{F}$ through Alg. 1, it is possible to efficiently apply them to new data matrices. In the case of t-SNE, however, the dimensionality reduction is done in relationship to the other samples in the dataset, so when a new data matrix is to be reduced in dimensions, the t-SNE algorithm must be re-executed. However, we include it in this exposition since it is a common algorithm for modern feature reduction. The main take-away from these experiments is the fact that the cPCA++ algorithm can obtain an explicit filter matrix, $\boldsymbol{F}$, and does not require a parameter sweep, which makes the algorithm much more efficient. At the same time, it is still able to discover the interesting structures in the target dataset.

III-B1 Synthetic Data

We first consider a synthetic data example to help illustrate the differences between the traditional PCA, cPCA, and cPCA++ algorithms. Consider the following original feature structure (i.e., prior to feature reduction) for an arbitrary $i$-th sample from a foreground or target data matrix $\widetilde{\boldsymbol{Z}_{f}}$:

[TABLE]

in which $k\in\{1,2,3,4\}$ denotes the class index, with sub-vectors $\boldsymbol{h}_{i,1}^{k},\boldsymbol{h}_{i,2}^{k},\boldsymbol{h}_{i,3}\in\mathbb{R}^{10\times 1}$, where $\boldsymbol{h}_{i,1}^{k}$, $\boldsymbol{h}_{i,2}^{k}$, $\boldsymbol{h}_{i,3}$ are independent. Now, it is assumed that $\boldsymbol{h}_{i,3}\sim\mathcal{N}(\mathbb{0}_{10},10\boldsymbol{I}_{10\times 10})$. On the other hand, the other two sub-vectors are more useful to differentiate across the different classes $k$:

[TABLE]

while $\boldsymbol{h}_{i,2}^{k}$ is distributed according:

[TABLE]

Knowing the above structure, we can deduce that to differentiate between the different classes $1\leq k\leq 4$, we need only use $\boldsymbol{h}_{i,1}^{k}$ and $\boldsymbol{h}_{i,2}^{k}$, while completely ignoring $\boldsymbol{h}_{i,3}$ as it does not help differentiate between the different classes. On the other hand, traditional PCA executed on the target data matrix $\widetilde{\boldsymbol{Z}_{f}}$ would readily pick features from $\boldsymbol{h}_{i,3}$ since they have the largest variance—even though they do not help in distinguishing between the different classes.

Instead, it is assumed that there is a background dataset that, while not containing the differentiating features, follows the following distribution:

[TABLE]

where the sub-vectors are independent. We observe that $\tilde{\boldsymbol{z}}_{b,i}$ only informs us as to the relative variance across the dimensions, but not the explicit mean values (i.e., feature structure) within the dimensions. It can be shown that the filters obtained by the cPCA++ algorithm are given by (please see App. A):

[TABLE]

where $\boldsymbol{e}_{n}$ denotes the standard basis vector of dimension $3$ (i.e., $\boldsymbol{e}_{n}=[\mathbb{0}_{n-1}^{\mathsf{T}},1,\mathbb{0}_{3-n}^{\mathsf{T}}]^{\mathsf{T}}$). Note that these are the ideal filters to extract the signal structure from the target dataset described by (46). We ran the cPCA++ method on this synthetic dataset, and compared the results with the other feature reduction methods in Fig. 3. Note that the t-SNE method operated over the target or foreground dataset. For all methods examined, we extract the dimensionality reduced dataset for $K=2$ dimensions. All classes in the figure are equally likely, and a total of $400$ samples were generated. In the top row of plots in Fig. 3, we ran the cPCA algorithm for three positive choices of the contrast parameter $\alpha$. We see that a contrast factor of $\alpha=2.7$ yields a correct separation among the classes, while the other values of the contrast factor fail to do so. In the bottom-left plot, we have the performance of the traditional PCA algorithm on the target dataset. It is in fact a special case of the cPCA method outlined in [1] for a contrast parameter of $\alpha=0$. We observe that traditional PCA fails to cluster the data. This is expected by our investigation above, since the traditional PCA method is designed to pick the directions that best explain the most variance of the dataset, which are directions that explain $\boldsymbol{h}_{3}$. In the bottom-center plot, we executed the t-SNE algorithm [40] on the foreground dataset for $250$ iterations (no significant improvement was obtained by running t-SNE longer). We see that t-SNE also fails to select the correct dimensions with which to visualize the data. Finally, in the bottom-right plot, we show the data plotted against the feature directions obtained by the cPCA++ algorithm. We observe that our feature extractor was able to perform the correct feature selection in one step, without the need for a parameter sweep. We note that while this is a synthetic example that can be solved analytically, it helps in illustrating how the method actually works. In the next few simulations, we compare and contrast the performance of the cPCA++, cPCA, traditional PCA, and t-SNE algorithms on more complex datasets.

III-B2 MNIST over High Variance Backdrop

In this experiment (also conducted in [1]), the MNIST dataset [41] is used and the $28\times 28$ images corresponding to digits [math] and $1$ are extracted. Utilizing traditional PCA processing, it is possible to separate these two particular digits.

However, the authors of [1] (as we do here), super-imposed these digits atop of high-variance grass backgrounds in order to mask the fine details that make the digit recognizable and distinct. Figure 4 shows six examples of these “target” images. We see that the digits themselves are difficult to see atop the high-variance background images. The task now is to apply the feature reduction techniques to still separate these digits.

In order to obtain a “background” dataset for the cPCA and cPCA++ methods, grass images (without the MNIST digits) randomly sampled from the ImageNet dataset [42] were used to estimate $\boldsymbol{R}_{b}$. Note that the same grass images that were used to generate the target images need not be utilized, but grass images in general will do.

In Fig. 5, we evaluate the performance of the dimensionality reduction techniques on the target images illustrated in Fig. 4. It can be seen that both traditional PCA and t-SNE have difficulty clustering the two digits, due to the high variance background masking the fine structure of the digits. On the other hand, the cPCA and cPCA++ methods are able to achieve some separation between the two digits, although in both cases there is some overlap. Again, note that cPCA++ does not require a hyperparameter sweep.

III-B3 Mice Protein Expression

This next example utilizes the public dataset from [43]. This dataset was also used in [1]. In this dataset, we have access to protein expression measurements from mice that have received shock therapy and mice that have not. We would like to analyze the foreground dataset (mice with shock therapy applied) to attempt to cluster mice that have developed Down Syndrome from those that have not. We also have access to a background dataset, which consists of protein expression measurements from mice that did not experience shock therapy. In Fig. 6 we illustrate the results of the different feature reduction techniques. It can be seen that traditional PCA is unable to effectively cluster the data. On the other hand, the t-SNE, cPCA, and cPCA++ methods are all able to achieve some separation of the data. Although t-SNE was able to achieve good separation in this particular example, we re-iterate the fact that t-SNE is not an appropriate method to be used to extract filters or function transformations from the high-dimensional space to the lower-dimensional space as the algorithm must be re-executed with new samples in order to map previously unmapped data. There is one other significant observation about this example. As indicated in the main paper, there is a very small number of samples in the background dataset. This caused the background covariance matrix, $\boldsymbol{R}_{b}$, to be rank-deficient. For this reason, diagonal loading was employed to ensure that the background covariance matrix is invertible.

III-B4 MHealth Measurements

In the next experiment (also conducted in [1]), we use the IMU data from [44]. In this dataset, a variety of sensors are used to monitor subjects that are performing squats, jogging, cycling, or lying still. The sensors include gyroscopes, accelerometers, and electrocardiograms (EKG). In the top-left subplot of Fig. 7, we see that performing traditional PCA on a dataset of subjects that are performing squats or cycling does not yield a reasonable clustering for those two activities. On the other hand, the optimal contrast parameter was used for cPCA, and the optimal cPCA result is illustrated in the top-right subplot of Fig. 7. For both the cPCA and cPCA++ methods, the background covariance matrix $\boldsymbol{R}_{b}$ was obtained from subjects lying still. We observe that the cPCA algorithm was able to cluster the two activities effectively. The t-SNE output is then shown in the bottom-left subplot of the figure. We see that t-SNE was unable to distinguish the two activities. Finally, the output of the cPCA++ method is shown in the bottom-right subplot of the figure. We see that cPCA++ is able to effectively cluster the two activities, without the need for a parameter sweep.

III-B5 Single Cell RNA-Seq of Leukemia Patient

Our final dataset in this section was obtained from [45], and this analysis was also conducted in [1]. The target or foreground dataset contains single-cell RNA expression levels of a mixture of bone marrow mononuclear cells (BMMCs) from a leukemia patient before and after stem cell transplant, and the goal is to separate the pre-transplant samples from the post-transplant samples. The original features were reduced to only 500 genes based on the value of the variance divided by the mean. Once these features are obtained, we execute the cPCA algorithm for different values of $\alpha$ and plot the result in the top row of Fig. 8. In this example, the background covariance matrix $\boldsymbol{R}_{b}$ was obtained from RNA-Seq measurements from a healthy individual’s BMMCs, for both the cPCA and cPCA++ methods. As noted in [1], this should allow the dimensionality reduction methods to focus on the novel features of the target dataset as opposed to being overwhelmed by the variance due to the heterogeneous population of cells as well as variations in experimental conditions. We observe that the data is close to being separable for the contrast value of $\alpha=3.5$. Next, we plot the traditional PCA output in the bottom-left subplot of Fig. 8, and observe that it is incapable of clustering the data. Next, in the bottom-center plot, we illustrate the output of the t-SNE algorithm, and see that while it does achieve some separation of the data, it is not on par with the cPCA algorithm. Finally, the cPCA++ output is illustrated in the bottom-right subplot of Fig. 8, and we see that it is able to obtain a similar clustering performance as compared to the optimal cPCA clustering result. Note that cPCA++ does not require a parameter sweep. The covariance matrix $\boldsymbol{R}_{b}$ was rank deficient in this example, so we utilized diagonal loading in order to make it invertible.

III-C Computational Time Performance Comparison

In this section, we will examine the computational time requirements for the various datasets studied in Sec. III-B. We neglect measuring the computation time required by the traditional PCA algorithm since it will be very similar to the cPCA++ algorithm and since traditional PCA is not an effective algorithm for discriminative dimensionality reduction, as we saw in the results of Sec. III-B. As shown in Table II, the proposed cPCA++ method is significantly faster than both the cPCA method and the t-SNE method. The reason why the cPCA++ method is on average (i.e., averaged over the experiments outlined in Sec. III-B) $51$ times faster than the cPCA method is that while both methods require the eigen-decomposition of similarly sized covariance matrices, the cPCA algorithm requires a sweep over the contrast hyper-parameter $\alpha$ (i.e., multiple eigen-decompositions). On the other hand, the cPCA++ method does not require this hyper-parameter tuning since the scale discrepancy is automatically resolved. It is also notable to see that cPCA++ is on average $428,654$ times faster than t-SNE. The t-SNE method has a very high running time due to the fact that it works by finding relationships between all the points in a dataset, and not through the determination of compact parameters. In addition, it is not very well suited for dimensionality reduction when future data samples will be obtained (since it has to be re-executed on the entire new dataset).

IV The cPCA++ Framework for Image Splicing Localization

In this section, we describe the general framework of the cPCA++ approach for image splicing localization. In this approach, we focus on detecting the spliced boundary or edge.

We take this approach of determining the spliced boundary rather than the spliced region or surface because there is an ambiguity in the determination of which image is the donor and the host. For example, consider Figure 9, which shows a probe image and the corresponding spliced boundary ground truth mask. Please note that we use the color white to denote a pixel belonging to the spliced boundary, and the color black to denote a pixel not belonging to the spliced boundary. It is possible to consider region A as the spliced region and region B as the authentic region, or vice versa. Thus, there is an ambiguity in how one labels the spliced and authentic regions in a given manipulated image. However, this ambiguity is resolved if we instead consider the spliced boundary (rather than the spliced surface/region).

The “training” phase consists of two main steps: 1) collect labeled or “training” samples/patches, 2) build a feature extractor based on the training samples. As we will see later, the cPCA++ approach does not require the training of a classifier, such as a support vector machine or random forest. However, we will still refer to patches used in this phase as “training” samples, since we are using their labeled information (i.e., spliced vs. authentic edges). In the testing phase, we are given a new image (i.e., not seen during the training phase), and the goal is to output a probability map that indicates the probability that each pixel belongs to the spliced boundary. The testing phase consists of the following main steps: 1) divide a given test image into overlapping patches, 2) extract the feature vector for each test patch, 3) obtain the probability output for each test patch, and 4) reconstruct the final output map for the given test image. The training and testing phases are elaborated below.

In the training phase, the first step is to collect labeled or “training” samples/patches from a given set of images. In order to do this, we need access to two masks for each training image: 1) the spliced surface ground truth mask and 2) a mask which highlights the superset of all edges/boundaries in the image (i.e., spliced boundaries as well as natural/authentic boundaries), referred to as the edge detection output mask. Note that although we utilize the surface ground truth mask in selecting the training patches, the output of the cPCA++ method will be edge-based (i.e., the output will represent an estimate of the true spliced edge/boundary). This will be explained in more detail in the next paragraph. We utilized the CASIA v2.0 dataset for training purposes. For the CASIA v2.0 dataset, the surface ground truth masks are not provided. We generated the ground truth masks using the provided reference information. In particular, for a given spliced image in the CASIA dataset, the corresponding donor and host images are provided, and we used this information to generate the ground truth mask. The second mask that is needed for collecting the training patches is the edge detection output mask. We utilized structured edge detection for detecting spliced and natural edges [46].

Once we have these two masks, we can proceed to collect the training patches. Each training image is divided into overlapping patches, and we then select target/foreground and background samples in the following manner. We would like target/foreground samples to be those that contain a boundary between a spliced region and an authentic region (referred to as a spliced boundary). In order to achieve this, we utilize the surface ground truth masks, and we select foreground samples to be those that have a splicing fraction that lies in a certain range (e.g., 30-70 percent of the patch is a spliced area/surface). On the other hand, background samples do not contain a boundary between a spliced and authentic region. We impose an additional constraint on the background samples so that they contain a minimum amount of authentic edges (as specified by the structured edge detector). Figure 10 shows an example of a foreground patch and a background patch. Note that although we are using the surface ground truth masks in selecting foreground and background samples, what is actually important is the spliced boundary/edge (or lack thereof).

The next step is to build a feature extractor based on the covariance matrices of the foreground and background training samples. Suppose we have collected $N_{b}$ background patches and $N_{f}$ foreground patches, and let the size of each patch be $n\times n\times c$, where $c$ is the number of channels. Since we are working with RGB images, $c$ will be equal to $3$. Each training patch can be flattened so that it is a $M\times 1$ vector, where $M$ is equal to $cn^{2}$ (or 3$n^{2}$ if $c=3$). We can represent the flattened background and foreground patches as data matrices $\widetilde{\boldsymbol{Z}}_{b}\in\mathbb{R}^{M\times N_{b}}$ and $\widetilde{\boldsymbol{Z}}_{f}\in\mathbb{R}^{M\times N_{f}}$, respectively, such that different columns of $\widetilde{\boldsymbol{Z}}_{b}$ (or $\widetilde{\boldsymbol{Z}}_{f}$) correspond to different training samples.

To see why it is necessary to use the cPCA++ method as opposed to traditional PCA for the image splicing localization task, consider the top 50 eigenvectors of the covariance matrix of $\widetilde{\boldsymbol{Z}}_{b}$ (i.e., authentic edges or background) and the top 50 eigenvectors of the covariance matrix of $\widetilde{\boldsymbol{Z}}_{f}$ (i.e., spliced edges or target/foreground). We compute the power of the foreground eigenvectors ($\boldsymbol{U}_{f}$) after projecting them onto the background subspace (spanned by $\boldsymbol{U}_{b}$) to quantify the similarity in the subspaces spanned by $\widetilde{\boldsymbol{Z}}_{b}$ and $\widetilde{\boldsymbol{Z}}_{f}$: $\|\boldsymbol{P}_{\boldsymbol{U}_{b}}\boldsymbol{U}_{f}\|_{F}^{2}/\|\boldsymbol{U}_{f}\|_{F}^{2}=93.25\%$. What this indicates is that the subspaces spanned by the top $50$ principal components of each dataset ($\boldsymbol{U}_{f}$ and $\boldsymbol{U}_{b}$) are mostly overlapping—hinting at the fact that it is likely that any classifier will be overwhelmed by the similarity of the features of these two datasets. For this reason, we will utilize the cPCA++ method instead to alleviate this situation and perform discriminative feature extraction.

We thus utilize Alg. 1, described in Sec. III, to build the feature extractor. We first center the data matrices $\widetilde{\boldsymbol{Z}}_{b}$ and $\widetilde{\boldsymbol{Z}}_{f}$ by subtracting their respective means. We then calculate the second order statistics $\boldsymbol{R}_{b}=\frac{1}{N_{b}}\boldsymbol{Z}_{b}\boldsymbol{Z}_{b}^{\mathsf{T}}$ and $\boldsymbol{R}_{f}=\frac{1}{N_{f}}\boldsymbol{Z}_{f}\boldsymbol{Z}_{f}^{\mathsf{T}}$, and then perform eigenvalue decomposition on $\boldsymbol{Q}=\boldsymbol{R}_{b}^{-1}\boldsymbol{R}_{f}$. If the matrix $\boldsymbol{R}_{b}$ is rank-deficient, we utilize diagonal loading in order to make it invertible. We then compute the top $K$ right-eigenvectors $\boldsymbol{F}$ of $\boldsymbol{Q}$. The matrix $\boldsymbol{F}\in\mathbb{R}^{M\times K}$ is used to extract features during the testing phase, and is referred to as the transform matrix.

In the testing phase, we divide each test image into overlapping patches. Similar to what was done during the training phase, we can flatten the test patches and represent them as a matrix $\widetilde{\boldsymbol{Z}}_{t}\in\mathbb{R}^{M\times N_{t}}$. We then utilize Algorithm 2 to obtain a probability output for each patch. The motivation behind Algorithm 2 is explained as follows. The cPCA++ method will naturally find filters that will yield a small norm for the reduced dimensionality features of background samples and a larger norm for foreground samples. This means that once the filters are obtained, they are expected to yield small valued output when applied to background patches and larger valued output when applied to foreground patches. This presents us with a very simple and efficient technique for obtaining the output: simply measure the output power after dimensionality reduction and then convert the raw value to a probability-based one. The full output map for a given test image can then be reconstructed by averaging the contributions from overlapping regions. After reconstruction, we perform an element-wise multiplication of the output map and the structural edge detection output mask. The reason for doing this is explained as follows. In this task, we are attempting to classify the spliced edge. Therefore, if the structural edge detector labeled a given pixel as non-edge (i.e., neither an authentic edge nor a spliced edge), then the corresponding pixel value in the cPCA++ output map should automatically be set to zero.

V Experimental Analysis

V-A Evaluation/Scoring Procedure

In this subsection, we explain the procedure used to score the output maps from the cPCA++ method, as well as from the Multi-task Fully Convolutional Network (MFCN)-based method [2] we compared against. The MFCN-based method is a state-of-the-art approach that also outputs an edge-based probability map, and thus can be directly compared with the proposed cPCA++ method. The MFCN achieved the highest score in the splicing localization task in the 2017 Nimble Challenge, and the second highest score in the 2018 Media Forensics Challenge (which are part of the DARPA MediFor program777https://www.darpa.mil/program/media-forensics). We evaluated the performance of the methods using the $F_{1}$ and Matthews Correlation Coefficient ($MCC$) metrics, which are popular per-pixel localization metrics. As noted previously, the edge-based approach that we propose avoids the ambiguity in the labeling of the spliced and authentic surfaces/regions. Because of the ambiguity in the surface-based labeling, most surface-based scoring procedures score the original output map as well as the inverted version of the output map, and select the one that yields the best score. However, in edge-based scoring procedures, there is no need to invert the edge-based output map, and the output can be scored directly. For a given spliced image, the $F_{1}$ metric is defined as

[TABLE]

where $TP$ represents the number of pixels classified as true positive, $FN$ represents the number of pixels classified as false negative, and $FP$ represents the number of pixels classified as false positive. In this case, the target/foreground patches can be viewed as the positive samples and the background patches can be viewed as the negative samples. The $F_{1}$ metric ranges in value from [math] to $1$, with a value of $1$ being the best. For a given spliced image, the $MCC$ metric is defined as

[TABLE]

The $MCC$ metric ranges in value from $-1$ to $1$, with a value of $1$ being the best.

V-B Experimental Results

We first compared the training times of the proposed cPCA++ method and the MFCN-based method, and found the proposed method to be significantly more efficient. The MFCN-based method requires a training time of approximately 11.5 hours on a NVIDIA GeForce GTX Titan X GPU, while the proposed cPCA++ method only requires a training time of approximately 1.5 hours on CPU (Intel Xeon Gold 6126). Please note that the cPCA++ code has not yet been optimized to make full utilization of the multi-core processor. Next, we evaluated the proposed cPCA++ method on the Columbia and Nimble WEB datasets and compared its performance with the MFCN-based method. Tables III and IV show the Matthews Correlation Coefficient ($MCC$) and $F_{1}$ scores, respectively. It can be seen that the cPCA++ method yields higher scores (in terms of both $MCC$ and $F_{1}$), as compared to the MFCN-based method.

Figures 11 and 12 show multiple examples of localization output from the Columbia and Nimble WEB datasets, respectively. Each row shows (from left to right) a manipulated or probe image with the spliced edge highlighted in pink, the structural edge detection output mask highlighting both spliced and authentic edges, the cPCA++ raw probability output map, and the MFCN-based raw probability output map. In these figures, it can be seen that the cPCA++ method yields a finer localization output than the MFCN-based method.

VI Conclusion

In conclusion, we proposed cPCA++, which is a new technique for discovering discriminative features in high-dimensional data. The proposed approach is able to discover structures that are unique to a target dataset, while at the same time suppressing uninteresting high-variance structures present in both the target dataset and a background dataset. The proposed cPCA++ approach is compared with a recently proposed algorithm, called contrastive PCA (cPCA), and we show that cPCA++ achieves similar discriminative performance in a wide variety of settings, even though it eliminates the need for the hyperparameter sweep required by cPCA. Following this discussion, the cPCA++ approach was applied to the problem of image splicing localization, in which the two classes of interest (i.e., spliced and authentic edges) are extremely similar in nature. In the context of this problem, the target dataset contains the spliced edges and the background dataset contains the authentic edges. We show that the cPCA++ approach is able to effectively discriminate between the spliced and authentic edges. The resulting method was evaluated on the Columbia and Nimble WEB splicing datasets. The proposed method achieves scores comparable to the state-of-the-art Multi-task Fully Convolutional Network (MFCN), and it does so very efficiently––without the need to iteratively update filter weights via stochastic gradient descent and backpropagation, and without the need to train a classifier.

Appendix A Derivation of Filters for Synthetic Example

The sample covariance of the background dataset specified in (52) will be:

[TABLE]

where the $\boldsymbol{A}\otimes\boldsymbol{B}$ operator denotes the Kronecker product of the matrices $\boldsymbol{A}$ and $\boldsymbol{B}$. On the other hand, the sample covariance of the target/foreground dataset is given by (when the classes are equally likely):

[TABLE]

Then, the matrix $\boldsymbol{Q}$ for the cPCA++ algorithm becomes:

[TABLE]

The matrix $\boldsymbol{Q}$ turns out to be symmetric (and thus is diagonalizable) and has the following block-diagonal structure:

[TABLE]

where

[TABLE]

The eigenvalues of the block diagonal matrix $\boldsymbol{Q}$ are given by the eigenvalues of its diagonal blocks888This is because an eigenvalue $\lambda$ of $\boldsymbol{Q}$ by definition satisfies $|\boldsymbol{Q}-\lambda\boldsymbol{I}_{30\times 30}|=0$. However, since $\boldsymbol{Q}$ is block diagonal, we have that $|\boldsymbol{Q}-\lambda\boldsymbol{I}_{30\times 30}|=|\boldsymbol{A}-\lambda\boldsymbol{I}_{10\times 10}|\cdot|\boldsymbol{B}-\lambda\boldsymbol{I}_{10\times 10}|\cdot|\boldsymbol{C}-\lambda\boldsymbol{I}_{10\times 10}|$ [38, p. 5]. Clearly, if $\lambda$ is an eigenvalue of $\boldsymbol{A}$, $\boldsymbol{B}$, or $\boldsymbol{C}$, it is also an eigenvalue of $\boldsymbol{Q}$ since it will also satisfy $|\boldsymbol{Q}-\lambda\boldsymbol{I}_{30\times 30}|=0$. There are a total of $30$ such values ($10$ per block), some with multiplicity greater than one., i.e., if the eigen-decomposition of $\boldsymbol{Q}$ is given by $\boldsymbol{Q}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{U}^{\mathsf{T}}$, then

[TABLE]

where the matrices $\boldsymbol{D}_{A}$, $\boldsymbol{D}_{B}$, and $\boldsymbol{D}_{C}$ contain the eigenvalues of the matrices $\boldsymbol{A}$, $\boldsymbol{B}$, and $\boldsymbol{C}$, respectively, along their diagonal and the eigenvalues are sorted in decreasing order in step $(a)$. It is then easy to verify that:

[TABLE]

where

[TABLE]

indicating that $\boldsymbol{u}_{1}$ and $\boldsymbol{u}_{2}$ are the principal eigenvectors of $\boldsymbol{Q}$, yielding (54).

Appendix B The cPCA++ Approach For Matrix Factorization and Image Denoising

In this appendix, we explore the use of the cPCA++ method for matrix factorization and image denoising. This denoising example was performed in [1] for the MNIST over grass dataset examined in Sec. III-B2. In this exercise, we are given a single foreground image (i.e., an image containing a digit overlayed on top of grass imagery), flattened into a vector $\boldsymbol{z}_{n}\in\mathbb{R}^{M\times 1}$, where $M=784$ and the subscript $n$ is used to indicate that $\boldsymbol{z}_{n}$ is noisy. Since $\boldsymbol{z}_{n}$ is a foreground image, it follows the distribution from (6), and thus we may seek the factorization of the expected value $\mathbb{E}[\boldsymbol{z}_{n}]=\boldsymbol{W}\boldsymbol{y}_{n}$ where $\boldsymbol{W}\in\mathbb{R}^{M\times K}$ and $\boldsymbol{y}_{n}\in\mathbb{R}^{K\times 1}$, and $K$ is the effective rank of the denoised image. This is a typical setup of matrix factorization methods. In our case, the matrix $\boldsymbol{W}$ contains dictionary atoms along its columns. These contain the underlying bases that make up the labeled foreground dataset $\boldsymbol{Z}_{f}$. On the other hand, the vector $\boldsymbol{y}_{n}$ contains the weighting of each dictionary atom (i.e., how much each dictionary atom contributes to the final image) and is a function of the filters $\boldsymbol{F}$ and the original foreground image $\boldsymbol{z}_{n}$. Traditionally, it is assumed that the weighting vector is sparse. We do not directly impose this, but when $K$ is chosen to be small, this is implicitly imposed. That is, we are approximating the foreground image (which is corrupted by the grass background in addition to the actual digit) as:

[TABLE]

where the atom $\boldsymbol{w}_{k}$ denotes the $k$-th column of the matrix $\boldsymbol{W}$, the weighting $y_{n,k}$ denotes the $k$-th element of the vector $\boldsymbol{y}_{n}$, and $\hat{\boldsymbol{z}}_{n}$ denotes the denoised version of the image $\boldsymbol{z}_{n}$. When $K\ll M$, we expect that $\boldsymbol{W}\boldsymbol{y}_{n}$ will yield a denoised signal as it will only utilize the few strong atoms forming $\boldsymbol{W}$ (which best describe the bases present in the foreground dataset). We also saw in Sec. III that the matrices $\boldsymbol{F}\in\mathbb{R}^{M\times K}$ and $\boldsymbol{W}\in\mathbb{R}^{M\times K}$ contain the principal eigenvectors of the matrices $\boldsymbol{R}_{b}^{-1}\boldsymbol{R}_{f}$ and $\boldsymbol{R}_{f}\boldsymbol{R}_{b}^{-1}$, respectively. In addition, the vector $\boldsymbol{y}_{n}$ is obtained via $\boldsymbol{y}_{n}=\boldsymbol{F}^{\mathsf{T}}\boldsymbol{z}_{n}$. This means that the denoised version of the image is given by $\hat{\boldsymbol{z}}_{n}=\boldsymbol{W}\boldsymbol{F}^{\mathsf{T}}\boldsymbol{z}_{n}$. We now attempt to denoise an image containing the digit [math] over grass background. This situation is shown in Fig. 13. Please note that we use the color white to denote pixels corresponding to the digit [math]. In the top-left plot, we show the original noisy digit. It is relatively difficult to see the digit over the background, so a significant amount of denoising is necessary. In all of the denoising methods, we set the number of components $K=3$. In the top-right plot, we show the denoising achieved by traditional PCA. In the bottom-left plot, we show the denoising achieved by the cPCA algorithm. Finally, in the bottom-right plot, we show the denoising performance of the cPCA++ method, which is achieved as the low-rank approximation $\boldsymbol{W}\boldsymbol{y}_{n}$. We observe that the output of cPCA++ is far less noisy than that of the other methods.

Acknowledgment

This material is based on research sponsored by DARPA and Air Force Research Laboratory (AFRL) under agreement number FA8750-16-2-0173. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of DARPA and Air Force Research Laboratory (AFRL) or the U.S. Government. Credits for the use of the CASIA Image Tampering Detection Evaluation Database (CASIA TIDE) v1.0 and v2.0 are given to the National Laboratory of Pattern Recognition, Institute of Automation, Chinese Academy of Science, Corel Image Database and the photographers http://forensics.idealtest.org.