Compressed Smooth Sparse Decomposition

1. Introduction

High-quality image sensing systems are widely used in manufacturing processes for product quality monitoring and fault diagnosis. The resolution and acquisition rate of such systems increase significantly benefiting from the rapid development of image sensing technology. For example, in a hot rolling process, an in situ image-based sensor can detect a micrometer-sized seam on a rolling bar at a speed of up to 225 miles per hour (Yan et al. 2018). For another example, to monitor solar activity, satellites can capture high-resolution solar images with a high acquisition rate, producing terabytes of data per day (Wang et al. 2018). To achieve real-time inspection, a large volume of high-resolution images needs to be transmitted and processed in real time. Such a large volume of high-resolution image data poses a big challenge not only on the speed of image processing algorithms but also, for the storage and transmission of the data itself.

Matrix decomposition-based image processing techniques are widely used in image-based process monitoring and anomaly detection. They achieve the goal by integrating the prior for background and anomaly components into the optimization problem. In terms of utilizing the low-rank and sparse property, robust principal component analysis was first proposed by Candès et al. (2011) to decompose a data matrix into low-rank and element-wise sparse components. One of its famous applications is dynamic foreground and static background separation (Bouwmans and Zahzah 2014). Following this approach, numerous algorithm variants have been proposed, including outlier pursuit (Xu et al. 2012), which aims to decompose the data matrix into a low-rank component and a column-wise sparse component, and low-rank plus compressed sparse decomposition (Mardani et al. 2013), which aims to decompose the data matrix into low-rank and compressed sparse components and so on. For utilizing smooth and sparse properties, smooth and sparse decomposition (SSD) methods (Minaee et al. 2015, Yan et al. 2017) are proposed for anomaly detection in images with smooth backgrounds. Following this approach, several explorations have been conducted, including spatiotemporal smooth sparse decomposition (ST-SSD) (Yan et al. 2018), additive tensor decomposition (Mou et al. 2021), and so on. By adopting the matrix decomposition approach, both the background and anomaly can be captured without detection time delay. However, because of the requirement of storage, transmission, and processing of the whole image signal, it cannot be applied in the scenario with low-transmission bandwidth but high processing speed requirements: for example, the solar flare detection application (Augusto et al. 2011).

To mitigate the data storage and transmission burden and improve sensing efficiency, compressive sensing (CS) (Candes et al. 2006) has been proposed, in which the data are directly collected in a compressed form and then, reconstructed accurately with high probability. More specifically, suppose that the original signal is a sparse vector $\mathit{y}\in {\mathbb{R}}^n$, and the main idea is to store and transmit a small set of compressive measurements $\mathit{y\text{'}}=\mathit{Ay}\in {\mathbb{R}}^p$, where $\mathit{A}\in {\mathbb{R}}^{p\times n}$ is an underdetermined sensing matrix ($p\ll n$) satisfying specific properties. Then, the original signal can be reconstructed from its compressed form $\mathit{y\text{'}}$, on which assorted image processing algorithms can be applied for defect detection and so on. For a comprehensive review of CS, please refer to Marques et al. (2018) and Rani et al. (2018). Even though promising, the naïve approach that tries to first reconstruct the image from the compressed measurement and then, apply matrix decomposition algorithms for anomaly detection has two issues.

1. For smooth plus sparse signals, the existence of such a sensing matrix $\mathit{A}$ satisfying specific properties is unknown.

2. The reconstruction process is usually computationally intensive (Marques et al. 2018), which on the other hand, slows down the overall computational speed of the defect detection algorithms.

Recently, to integrate the CS with matrix decomposition algorithms, Waters et al. (2011) proposed an SpaRCS method to recover low-rank and sparse matrices directly from compressive measurements, and Tanner and Vary (2020) gave a rigorous performance discussion. However, those methods do not consider the smooth plus sparse decomposition problems and are not efficient in dealing with high-order data.

In this paper, we discuss the possibility of adopting compressive data acquisition systems for image-based quality monitoring and fault detection in applications where the background is smooth, and anomalies are sparse. To achieve so, we propose a compressed smooth sparse decomposition (CSSD) framework. In this framework, the signal processing algorithms are directly applied to the compressed data, and no reconstruction step is needed. By doing so, a significant cost reduction in sensing, storage, and transmission as well as a boost in the speed of image processing algorithms can be achieved with negligible performance loss. We also established the theoretical foundation of adopting such a compressive data acquisition system for smooth plus sparse signals as well as the performance guarantee of the proposed algorithm. To further improve its performance in high-order scenarios, a Kronecker compressed smooth sparse sensing (KronCSSD) is proposed.

The remainder of this paper is organized as follows. In Section 2, we present the CSSD framework. In Section 3, we use simulation studies to validate the proposed framework. In Section 4, we demonstrate the proposed framework using several case studies. Finally, Section 5 concludes the paper.

2. Compressed Smooth Sparse Decomposition Framework

In this section, we will present the proposed CSSD framework. As mentioned in Section 1, we aim to design a fast and data-efficient method for sparse anomaly (sparse signal component) detection in signals with a smooth background (smooth signal component). For simplicity, we first discuss the methodology for one-dimensional (1D) signals and then, generalize it to n-dimensional images. Figure 1 provides an overview of the proposed methodology. There are three stages in the proposed methodology. (i) The signal is acquired in its compressed form through compressive measurement. (ii) Then, the compressed data are transmitted to the server. (iii) Finally, a decomposition algorithm will be applied directly to the compressed signal to decompose it into its corresponding smooth and sparse signal components.

Mathematically, let $\mathit{y}\in {\mathbb{R}}^n$ be a smooth plus sparse signal (which will be defined formally in Section 2.1). We aim to store and transmit a small set of compressive measurements $\mathit{y\text{'}}$ (i.e., $\mathit{y\text{'}}=\mathit{Ay}$) and then, reconstruct the smooth and sparse signal components from $\mathit{y\text{'}}$. To achieve so, there are three questions to be addressed.

1. What is a smooth plus sparse signal?

2. How do we compress such a signal?

3. How do we reconstruct such a signal from compressive measurements?

The remainder of this section is organized as follows to answer those three questions. We start with a formal definition of 1D smooth plus sparse signals in Section 2.1. Based on that, we introduce the compressive measurement method and discuss its theoretical properties for such signals in Section 2.2. In Section 2.3, we present the proposed CSSD framework that can directly reconstruct the smooth and sparse signal components from the compressive measurement by solving the following optimization problem:

where $\mathit{B}$ and ${\mathit{B}}_a$ are bases, $\mathit{\theta }$ and ${\mathit{\theta }}_a$ are corresponding coefficients, and ${ϵ}_1$ is the bound for measurement error. The reconstruction accuracy is also characterized theoretically. Then, we generalize the CSSD algorithm to n dimensions using Kronecker compressive sensing (KCS) (Duarte and Baraniuk 2011) and propose a KronCSSD formulation in Section 2.4. Finally, in Section 2.5, we present the advantage of the proposed framework and give the strategy of selecting the compressive ratio, tuning parameters, and bases in practice.

2.1. The Set of Smooth Plus Sparse Signals

In this section, we define the set of smooth plus sparse signals mathematically. The smooth signals originate from the spline smoothing (De Boor and De Boor 1978, Eilers and Marx 1996), where the raw signal is approximated by a linear combination of a set of spline basis functions for smooth interpolation and denoising, and a spline regression technique is usually utilized. To improve the regression robustness with respect to outliers, outliers are explicitly accounted for in the regression model as a sparse component of the raw signal (Giannakis et al. 2011, Mateos and Giannakis 2011). This idea was further generalized to incorporate the special structure of outliers (Yan et al. 2017, 2018).

As a summary, a 1D signal $\mathit{y}\in {\mathbb{R}}^n$ is defined as a smooth plus sparse signal if it can be decomposed into two signal components: (i) a smooth signal $\mathit{m}\in {\mathbb{R}}^n$ in a low-dimensional subspace spanned by a set of smooth bases (i.e., $\mathit{m}=\mathit{B}\mathit{\theta }$, where $\mathit{B}\in {\mathbb{R}}^{n\times r}$ is a basis matrix with $r\ll n$) and (ii) a sparse signal $\mathit{a}\in {\mathbb{R}}^n$ in a relatively high-dimensional subspace spanned by a set of predefined bases, of which the coefficients admit sparse property (i.e., $\mathit{a}={\mathit{B}}_a{\mathit{\theta }}_a$, where ${\mathit{B}}_a\in {\mathbb{R}}^{n\times q}$ is a basis matrix with $q\le n$ and ${\mathit{\theta }}_a\in {\mathbb{R}}^q$ is an $s$-sparse vector; i.e., $\parallel {\mathit{\theta }}_a{\parallel }_0\le s$). Given a smooth plus sparse signal $\mathit{y}$, we define the aforementioned decomposition as SSD (i.e., $\mathit{y}=\mathit{m}+\mathit{a}$).

To ensure the uniqueness of SSD in a nontrivial case when $n\le r+q$, the following definition is introduced.

Definition 2.1.

The local support property of ${\mathit{B}}_a$, ${\mathit{B}}_a\in {\mathbb{R}}^{n\times q}$. ${\mathit{B}}_a$ only has local support such that each column of ${\mathit{B}}_a$ only has nonzero values inside a specific interval. The length of this interval is defined as $l\left({\mathit{B}}_a\right)$.

Notice that $l\left({\mathit{B}}_a\right)\in \left\{1,\dots ,n\right\}$. The local support property with a small $l$ ensures the sparsity of $\mathit{a}$. For example, the B-spline basis has a local support property (Unser 1999).

Definition 2.2.

The incoherence condition of $\mathit{B}, \mathit{B}\in {\mathbb{R}}^{n\times r}$. Let $\mathit{B}=\mathit{U\Sigma }{\mathit{V}}^T$ be the reduced singular value decomposition (SVD) of $\mathit{B}$, where $\mathit{U}\in {\mathbb{R}}^{n\times r}, \mathit{\Sigma }\in {\mathbb{R}}^{r\times r}$, and $\mathit{V}\in {\mathbb{R}}^{r\times r}$. Its incoherence condition parameter $\mu \left(\mathit{B}\right)$ is defined as the smallest value such that

$$\underset{i\in \left\{1,\mathrm{..},r\right\}}{\max }\parallel {\mathit{U}}^T{\mathit{e}}_i{\parallel }_2\le \sqrt{\frac{\mu \left(\mathit{B}\right)r}{n}},$$

where ${\mathit{e}}_i$ is the $i$ th standard basis vector in ${\mathbb{R}}^n$.

Notice that $\mu \left(B\right)\in [1,\sqrt{n/r}]$. The incoherent condition with a small $\mu$ ensures that the $\mathit{m}$ is not sparse (Candès et al. 2011).

The following theorem ensures the uniqueness of the SSD decomposition.

Theorem 2.1.

If $\mu \left(\mathit{B}\right)<n{\left(2rsl\right)}^{-1}$, then the SSD decomposition is unique with respect to $\mathit{m}$ and $\mathit{a}$.

Theorem 2.1 gives the condition that the smooth plus sparse signal can be uniquely decomposed into a smooth part and a sparse part. The proof of Theorem 2.1 is in Appendix A.

Formally, we define the set of smooth plus sparse signals as follows.

Definition 2.3.

The set of smooth and sparse signals is defined as $M{S}_{r,s,\mu ,l}:$

$$\begin{array}{l}M{S}_{r,s,\mu ,l}=\{\mathit{y}\in {\mathbb{R}}^n|\mathit{y}=\mathit{B}\mathit{\theta }+{\mathit{B}}_a{\mathit{\theta }}_a, {\mathit{B}}_a\in {\mathbb{R}}^{n\times q}, l\left({\mathit{B}}_a\right)=l, {\mathit{\theta }}_a\in {\mathbb{R}}^q,\parallel {\mathit{\theta }}_a{\parallel }_0\le s, \mathit{B}\in {\mathbb{R}}^{n\times r},\hfill \\ \mu \left(\mathit{B}\right)=\mu <n{\left(2rsl\right)}^{-1},\mathit{\theta }\in {\mathbb{R}}^r\}.\hfill \end{array}$$

For such a smooth plus sparse signal, how to compress it while ensuring the reconstruction performance will be discussed in the next section.

2.2. Compressive Sensing for Smooth Plus Sparse Signals

As mentioned in Section 1, to ensure the reconstruction performance, the sensing matrix $\mathit{A}$ has to satisfy the so-called restricted isometry property (RIP) (Candes 2008). It has been proven that a random matrix can satisfy the RIP property for the sparse signal (Candes 2008), the low-rank signal (Recht et al. 2010), and the rank plus sparse signal (Tanner and Vary 2020). However, the existence of such a matrix for the smooth plus sparse signal, which is the foundation of adopting compressed data acquisition techniques in applications with smooth background and sparse anomalies, is still unknown. In this section, we will discuss the existence of such a sensing matrix. Before stating the result, we will first present the relevant definitions that are necessary to derive the main result.

Definition 2.4.

RIP for $M{S}_{r,s,\mu ,l}$. Let $\mathit{A}\in {\mathbb{R}}^{p\times n}$ be a linear measurement matrix. For every quadruple $\left(r,s,\mu ,l\right)$, define the restricted isometry constant (RIC) ${\delta }_{r,s,\mu ,l}$ to be the smallest positive constant such that

$$\left(1-{\delta }_{r,s,\mu ,l}\right)\parallel \mathit{y}{\parallel }_2\le \parallel \mathit{Ay}{\parallel }_2\le \left(1+{\delta }_{r,s,\mu ,l}\right)\parallel \mathit{y}{\parallel }_2, \forall \mathit{y}\in M{S}_{r,s,\mu ,l}.$$

If such a ${\delta }_{r,s,\mu ,l}\in \left(0,1\right)$ exists, we say that $\mathit{A}$ satisfies the RIP.

Theorem 2.2.

Suppose that ${\delta }_{r,2s,\mu ,l}<1$ for some integer $r,s,l\ge 1$ and positive numbers $\mu <n{\left(2rsl\right)}^{-1}$; then, there is a ${\mathit{y}}_0$ in the set $M{S}_{r,s,\mu ,l}$, which is the only solution for $\mathit{A}{\mathit{y}}_0=\mathit{b}.$

Theorem 2.2 guarantees the uniqueness of the smooth plus sparse signal that satisfies the sensing equation when $\mathit{A}$ satisfies the RIP. The proof of Theorem 2.2 is in Appendix B.

Next, we prove that for the set of smooth plus sparse signals, $M{S}_{r,s,\mu ,l}$, there exists such a matrix $\mathit{A}$ satisfying the RIP property with RIC = ${\delta }_{r,s,\mu ,l}$ with high probability.

Notice that the RIP for a matrix is difficult to verify. A suitable set of random matrices that obey the RIP for the set of sparse vectors with high probability (Recht et al. 2010, Tanner and Vary 2020) is defined as follows.

Definition 2.5.

Nearly isometric matrices (Baraniuk et al. 2008). Let $\mathit{A}\in {\mathbb{R}}^{p\times n}$ be a random variable that takes values in linear maps from ${\mathbb{R}}^n$ to ${\mathbb{R}}^p$; then, for any $\mathit{y}\in {\mathbb{R}}^n$, $\mathit{A}$ is nearly isometrically distributed if

1. $\mathbb{E}[\parallel \mathit{Ay}{\parallel }_2^2]=\parallel \mathit{y}{\parallel }_2^2$ and

2. $Pr(|\parallel \mathit{Ay}{\parallel }_2^2-\parallel \mathit{y}{\parallel }_2^2|\ge ϵ\parallel \mathit{y}{\parallel }_2^2)\le 2e^{-p{c}_0\left(ϵ\right)}, 0<ϵ<1,$

   where $c_0\left(ϵ\right)$ is a constant that only depends on $ϵ$.

The $p\times n$ matrix with independent, identically distributed (i.i.d.) Gaussian entries satisfies those two properties (Baraniuk et al. 2008) $($i.e., ${\mathit{A}}_{ij}\sim \mathcal{N}\left(0,\frac{1}{p}\right),$ with $c_0\left(ϵ\right)={ϵ}^2/4-{ϵ}^3/6$$)$. There are also other distributions satisfying the nearly isometric property, such as the $p\times n$ matrix with i.i.d. Bernoulli entries and their related distribution (Baraniuk et al. 2008).

Then, the following theorem states that the nearly isometric matrices can also serve as the sensing matrix for smooth plus sparse signals and gives the magnitude of the number of linear measurements.

Theorem 2.3.

Let $\mathit{A}\in {\mathbb{R}}^{p\times n}$ be a matrix from the families described in Definition 2.5. Furthermore, assume that $\mu <n{\left(2rsl\right)}^{-1}$ and the basis matrix ${\mathit{B}}_a$ for the sparse signal component satisfies the RIP with RIC ${\delta }_{{\mathit{B}}_a,s}\in \left(0,1\right)$: that is, ${\delta }_{{\mathit{B}}_a,s}$ to be the smallest positive constant such that

$$\left(1-{\delta }_{{\mathit{B}}_a,s}\right)\parallel {\mathit{\theta }}_a{\parallel }_2\le \parallel {\mathit{B}}_a{\mathit{\theta }}_a{\parallel }_2\le \left(1+{\delta }_{{\mathit{B}}_a,s}\right)\parallel {\mathit{\theta }}_a{\parallel }_2, \forall {\mathit{\theta }}_a\in \left\{{\mathit{\theta }}_a\in {\mathbb{R}}^q,\parallel {\mathit{\theta }}_a{\parallel }_0\le s\right\}.$$

For a given $\delta \in \left(0,1\right)$, there exists constants $c_1,c_2>0$ depending only on $\delta$, such that the RIC for $M{S}_{r,s,\mu ,l}$ is upper bounded by $\delta ,$ with the probability of at least $1-\exp \left(-c_{1}p\right)$, whenever

$$p\ge c_2\left(\ln 2+r\ln \frac{24}{\delta }{\tau }_1+s\left(1+\ln \frac{24}{\delta }{\tau }_0+\ln \frac{n}{s}\right)\right),$$

where $\eta =\sqrt{\frac{\mu \mathit{rsl}}{n}}$, ${\tau }_0=\frac{1}{\sqrt{\left(1-{\delta }_{{\mathit{B}}_a,s}\right)\left(1-{\eta }^2\right)}}$, ${\tau }_1=\parallel {\mathit{B}}^{†}{\parallel }_2\left(1+\frac{1}{\sqrt{1-{\eta }^2}}\right)$, and ${\mathit{B}}^{†}={\left({\mathit{B}}^{\mathit{T}}\mathit{B}\right)}^{-1}{\mathit{B}}^{\mathit{T}}$.

Theorem 2.3 states that if the nearly isometric matrix is selected as the sensing matrix, the RIC for $M{S}_{r,s,\mu ,l}$ is upper bounded with high probability. In practice, $\mathit{B}$ and ${\mathit{B}}_a$ are prespecified based on the understanding/engineering knowledge of the process (please refer to Section 2.5.4 for more detail). Theorem 2.3 also provides a guidance on the magnitude of linear measurements $p$, which determines the compressive measurement matrix $\mathit{A}$.

The proof of Theorem 2.3 is in Appendix C.

In this section, we answered two fundamental questions. (i) How do we compress the smooth plus sparse signal (design the compressive measurement matrix $\mathit{A}$)? (ii) How many linear measurements are needed to preserve the information in smooth plus sparse signals with high probability?

In the next section, we will discuss the problem formulation to recover the smooth and sparse signal components simultaneously from the compressed signal using the CSSD framework.

2.3. Compressed Smooth Sparse Decomposition

As mentioned in Section 1, one way to recover the smooth component and sparse component is first reconstructing the compressed image and then, using the SSD algorithm. However, this will slow down the speed of the defect detection algorithm. Instead, we propose to solve the one-step convex relaxation Problem (1). Natural questions are if we can recover the smooth and sparse signals from the compressed measurement ${\mathit{y}}^{\prime }$ by solving Problem (1) and what the accuracy is. The following theorem guarantees the recovery performance of the proposed convex relaxation Problem (1).

Theorem 2.4.

Let $\mathit{A}\in {\mathbb{R}}^{p\times n}$ be a matrix from the families described inDefinition 2.5. Let the signal $\mathit{y}=$${\mathit{m}}_0+{\mathit{a}}_0\in M{S}_{r,s,\mu ,l}$, where ${\mathit{m}}_0=\mathit{B}{\mathit{\theta }}_0$ and ${\mathit{a}}_0={\mathit{B}}_a{\mathit{\theta }}_{a0}$. Assume that Problem (1) is feasible, and let the optimal solution be ${\mathit{m}}^{*}=\mathit{B}{\mathit{\theta }}^{*}$ and ${\mathit{a}}^{*}={\mathit{B}}_a{\mathit{\theta }}_a^{*}$. Assume that the basis matrix $B_a$ for sparse signal component satisfies the RIP with RIC ${\delta }_{B_a,2s}\in \left(0,1\right)$. Let $a=\left(1+{\alpha }_1+{\alpha }_2\right){\gamma }^2+2{\alpha }_2+2$ and c $=1-{\gamma }^2{\alpha }_1{\alpha }_2-{\alpha }_2{}^2$, where ${\alpha }_1=\frac{\eta }{1-{\eta }^2}$, ${\alpha }_2=\frac{\sqrt{2}\eta }{1-2{\eta }^2}$, $\gamma =\sqrt{\frac{1+{\delta }_{B_a},2s}{1-{\delta }_{B_a,2s}}}$, and $\eta =\sqrt{\frac{\mu \mathit{rsl}}{n}}$. Suppose that $r,s,l\in \mathbb{N}$ and $\mu <n{\left(2rsl\right)}^{-1}$, such that $c>0$ and ${\delta }_{r,3s,\mu ,l}\in \left(0,c/a\right)$; then,

$${\parallel {\mathit{a}}_0-{\mathit{a}}^{*}\parallel }_2={\parallel {\mathit{B}}_a{\mathit{\theta }}_{a0}-{\mathit{B}}_a{\mathit{\theta }}_a^{*}\parallel }_2\le C_a{ϵ}_1$$

and

$${\parallel {\mathit{m}}_0-{\mathit{m}}^{*}\parallel }_2={\parallel \mathit{B}{\mathit{\theta }}_0-\mathit{B}{\mathit{\theta }}^{*}\parallel }_2\le C_m{ϵ}_1,$$

where

$$C_a=\frac{\left(1+{\gamma }^2\right)\left(1+{\alpha }_2\right)\sqrt{1+{\delta }_{r,3s,\mu ,l}}}{c-a{\delta }_{r,3s,\mu ,l}}$$

and

$$\begin{array}{c}{C}_m=\frac{\sqrt{1+{\delta }_{r,3s,\mu ,l}}+\left({\delta }_{r,3s,\mu ,l}+\frac{{\gamma }^2}{1+{\gamma }^2}{\alpha }_1+\frac{1}{1+{\gamma }^2}{\alpha }_2\right)C_a}{\left(1-{\delta }_{r,3s,\mu ,l}\right)}.\hfill \end{array}$$

Theorem 2.4 gives the conditions that the proposed convex relaxation Problem (1) can recover the true smooth and sparse signal up to a constant times the noise bound.

The proof of Theorem 2.4 is in Appendix D.

Notice that one advantage of the proposed CSSD framework is that it is compatible with existing decomposition algorithms. For example, for the SSD algorithm proposed by Yan et al. (2017), the problem formulation becomes

which can be solved efficiently by using the algorithm proposed by Yan et al. (2017).

2.4. Kronecker Compressed Smooth Sparse Decomposition

In the previous section, we discussed the CSSD framework for the 1D signal. In this section, we will generalize the proposed CSSD to KronCSSD framework for high-order tensor data.

Let $\mathcal{Y}\in \hspace{0.17em}{\mathbb{R}}^{n_1\times \cdots \times n_d}$ be the original signal and $\mathit{y}\in {\mathbb{R}}^N$ be its corresponding vectorized signal (i.e.,$\mathit{y}=\text{vec}\left(\mathcal{Y}\right)$ and $N={\prod }_{i=1}^d{n}_i$). Let $\mathit{A}\in {\mathbb{R}}^{p\times N}$ be a measurement matrix satisfying the RIP. Let $\mathit{y}\prime \in {\mathbb{R}}^p$ be the compressed data, such that $\mathit{y}\prime =\mathit{Ay}$. $\mathit{B}={\otimes }_{i=1}^d{\mathit{B}}_i$ and ${\mathit{B}}_a={\otimes }_{i=1}^d{\mathit{B}}_{ai}$ are the known bases for smooth and sparse components, respectively. Notice that problem formulation (1) can still be used for recovering the high-order smooth and sparse signal components from the compressive measurement. However, there are two issues. (i) In practice, the global CS measurements matrix $\mathit{A}$ is hard to realize using the CS device (Duarte and Baraniuk 2011). (ii) The resulting bases $\mathit{B}$ and ${\mathit{B}}_a$ can be extremely large as the dimension of the data increases, which will not only cause a big challenge in the storage of such a large matrix but also, result in the computational issue in handling such large matrices.

In reality, the high-dimensional tensor data usually have low-rank properties along each mode, which have been extensively exploited in tensor low-rank modeling techniques, such as CANDECOMP/PARAFAC (CP)/Tucker decompositions (Kolda and Bader 2009). This makes it possible for designing a sensing matrix for each mode, which is called Kronecker CS (Duarte and Baraniuk 2011). Inspired by the KCS method, we propose the KronCSSD framework formulation as follows.

Let ${\mathit{A}}_i\in {\mathbb{R}}^{p_i\times n_i}$ be a measurement matrix with RIP for each mode of the tensor; then, we have the following formulation,

where ${\mathcal{Y}}^{\prime }=\mathcal{Y}{\times }_1{\mathit{A}}_1{\times }_2\cdots {\times }_d{\mathit{A}}_d$ is the compressive measurement. $\mathbf{\Theta }\in {\mathbb{R}}^{r_1\times \cdots \times r_d}$ and ${\mathbf{\Theta }}_a\in {\mathbb{R}}^{q_1\times \cdots \times q_d}$ are the basis coefficients for the smooth and sparse signal components, respectively.

The proposed KronCSSD framework is a nontrivial generalization of the CSSD framework for 1D signals. Its performance will be shown empirically in simulation and case studies. The theoretical discussion of the KronSSD framework is left for future work.

For example, for a two-dimensional (2D) image $\mathit{Y}\in {\mathbb{R}}^{n_1\times n_2}$, when adopting the SSD algorithm (Yan et al. 2017), the problem formulation becomes

where $\mathit{Y}\prime$ is the compressed image (i.e., ${\mathit{Y}}^{\prime }={\mathit{A}}_1\mathit{Y}{\mathit{A}}_2^T$). It can be solved efficiently by using the algorithm proposed by Yan et al. (2017).

2.5. Discussion

2.5.1. Advantages of the Proposed CSSD/KronCSSD Methods.

In this section, we will give a brief discussion about the advantages of the proposed CSSD/KronCSSD methods. We define the compressive ratio as $c={\prod }_{i=1}^d{p}_i/n_i$. The smaller the compressive ratio, the fewer data will be transmitted. The proposed CSSD/KronCSSD methods have the following characteristics.

1. We propose to directly acquire the compressed image, which not only reduces the sensing cost but also, reduces the data transmission and storage cost by ${\prod }_{i=1}^d{p}_i/n_i$ times.

2. The smooth and sparse signal components can be recovered by solving a smaller-scale convex optimization problem with input data ${\prod }_{i=1}^d{p}_i/n_i$ times smaller than that of the original problem, which significantly boosts the computation.

2.5.2. Compressive Ratio Selection in Practice.

Theorems 2.3 and 2.4 show that the 1D smooth plus sparse signal can be recovered with high probability from a compressive measurement if the compressive ratio is above a specific threshold. However, there are several parameters (such as $c_1,c_2$) that are difficult to obtain when calculating the threshold in practice. We propose a practical procedure of selecting the compressive ratio utilizing the historical data. Suppose some training signals with their real background and anomalies are available; the compressive ratio can be selected with the following guidelines.

1. If there is a requirement for reconstruction accuracy for smooth and sparse signal components, then $\widehat{p}$ is chosen as the smallest value that satisfies such a requirement.

2. If there is no such requirement, we recommend choosing the compressive ratio corresponding to the sharp change point of the slope of the loss function-compressive ratio curve. This point exists because of the existence of such a threshold, after which the reconstruction with high probability is guaranteed by Theorems 2.3 and 2.4. We will demonstrate this in the simulation study in Section 3.1.

For high-order tensor data, the selection of the $p_i$ for each mode can be challenging. We provide some empirical guidelines as follows.

1. If the smoothness of the background is similar along different modes, a unified compressive ratio is recommended.

2. If the smoothness of the background is different, more sensing budget should be allocated to the mode along which the background is less smooth.

3. We recommend fixing the ratio among $p_i$ along different modes and adopting steps (i) and (ii) for 1D signal to determine the compressive ratio.

2.5.3. Tuning Parameter Selection in Practice.

Notice that in problem formulations (1) and (3), there is a hyperparameter ${ϵ}_1$, indicating the bound for measurement noise. If the measurement error bound is known from the accuracy of the measurement device, it can be directly used here. Otherwise, a crossvalidation step using historical data is recommended. Similarly, there is a tuning parameter $\lambda$ in their corresponding Lagrangian form Equations (2) and (4), which controls the sparsity of the decomposed anomaly. Crossvalidation can be used to determine this parameter. For more detail, please refer to Yan et al. (2017).

2.5.4. Bases Selection in Practice.

The bases $\mathit{B}$ and ${\mathit{B}}_a$ are prespecified based on the understanding/engineering knowledge of the process. The selection of such bases is discussed in detail in section 3.4 of Yan et al. (2017). In general, any smooth basis, such as splines or kernels, can be used for the background. For sparse anomalies, such as small regions scattered over the background or in the form of thin lines, an identity basis is recommended. Linear (quadratic) B splines are recommended for anomalous regions with sharp corners (curved boundaries). However, to ensure the uniqueness of the smooth sparse decomposition, we do require the bases $\mathit{B}$ and ${\mathit{B}}_a$ to satisfy specific properties, such as those mentioned in Definitions 2.1 and 2.2 and Theorem 2.1, which should be checked for selected bases. We will demonstrate this in the simulation study in Section 3.1.

3. Simulation Study

In this section, we will demonstrate the proposed CSSD and KronCSSD framework with simulation studies. First, the CSSD method is applied on 1D signals in Section 3.1, and then, we apply the KronCSSD method on 2D images in Section 3.2.

3.1. CSSD on a 1D Signal

A 1D signal $\left(\mathit{y}\in {\mathbb{R}}^n\right)$ is assumed to be a superposition of the smooth signal component ($\mathit{m}=\mathit{B}\mathit{\theta }$), the sparse signal component $\left(\mathit{a}={\mathit{B}}_a{\mathit{\theta }}_a\right)$, and noise $\mathit{e}$ (i.e., $\mathit{y}=\mathit{m}+\mathit{a}+\mathit{e}$). To demonstrate the proposed CSSD framework, we will first conduct compressive data acquisition on the raw signal (i.e., $\mathit{y}\prime =\mathit{Ay}$). Then, the data reconstruction and decomposition are achieved in one step. The performance is evaluated by the relative error between the true signal $\mathit{a}$ and the reconstructed one $\widehat{\mathit{a}}$ (i.e., ${\parallel \mathit{a}-\widehat{\mathit{a}}{\parallel }_2/\parallel \mathit{a}\parallel }_2$ and $\parallel \mathit{m}-\widehat{\mathit{m}}{\parallel }_2/\parallel \mathit{m}{\parallel }_2$ for the smooth background).

In the simulation study, we generate a 1D signal from $M{S}_{r,s,\mu ,l}$ and let $n=\mathrm{1,000}$. The smooth background is generated from a random linear combination of B-spline bases with three knots $\left(r=10\right)$ (i.e., $\mathit{m}=\mathit{B}\mathit{\theta }$, where $\mathit{B}\in {\mathbb{R}}^{n\times r}$ and $\mathit{\theta }\in {\mathbb{R}}^r$ is a random vector such that ${\mathit{\theta }}_i\sim \mathit{N}\left(0,1\right)$, $i\in \left\{1,\dots ,4\right\}$). The incoherence condition parameter $\mu =\mu \left(\mathit{B}\right)=0.82$. The sparse signal component is generated by a sparse random linear combination of degree 2 B-spline bases with 500 knots $\left(q=500,l=4\right)$: that is, $\mathit{a}={\mathit{B}}_a{\mathit{\theta }}_a$, where ${\mathit{B}}_a\in {\mathbb{R}}^{n\times q}$ and ${\mathit{\theta }}_a\in {\mathbb{R}}^q$ is a four-sparse random vector ($s=4$) such that its nonzero elements follow i.i.d. standard normal distribution. Notice that the RIC, ${\delta }_{{\mathit{B}}_a,s}$ for matrix ${\mathit{B}}_a$, is hard to calculate in general. However, there is a loose upper bound that can be used, which is ${\delta }_{{\mathit{B}}_a,s}\le {\delta }_{{\mathit{B}}_a,p}=\max \left\{{\lambda }_{\mathit{max}}-1,1-{\lambda }_{\mathit{min}}\right\}=0.60$, where ${\lambda }_{\mathit{max}}$ and ${\lambda }_{\mathit{min}}$ are the maximum and minimum singular values of ${\mathit{B}}_a$, respectively. The noise signal generated as a random vector $\mathit{e}\in {\mathbb{R}}^n$ is a random vector such that ${\mathit{e}}_i\sim \mathit{N}\left(0, {0.001}^2\right)$, $i\in \left\{1,\dots ,n\right\}$. Figure 2 shows a sample raw signal and its smooth, sparse components.

Before we state the result, we first check the assumption of Theorems 2.3 and 2.4. For Theorem 2.3, it is easy to check that $0.82=\mu \le \frac{n}{2rsl}=\frac{1,000}{2\times 10\times 4\times 4}=6.25$ and also, that ${\delta }_{{\mathit{B}}_a,s}\le 0.60\in \left(0,1\right)$. For Theorem 2.4, ${\delta }_{{\mathit{B}}_a,2s}\le 0.60\in \left(0,1\right)$, $c=0.37>0$, and ${\delta }_{r,3s,\mu ,l}\in \left(0,0.04\right)$. Notice that the range for ${\delta }_{r,3s,\mu ,l}$ is small in this case because of the loose upper bound for ${\delta }_{{\mathit{B}}_a,2s}$, which can be further improved. This indicates that the smooth and sparse signal can be recovered by the proposed algorithm with high probability provided that the compressive ratio is above a specific threshold, which is demonstrated by the following observation.

The simulation is repeated 100 times, and the average log relative error for the background and sparse signal components with respect to compressive ratio are shown in Figure 3. We can observe a large error at the beginning, and it drops very fast with an increase of the compressive ratio. Then, above a threshold (0.1 approximately) of compressive ratio, the error becomes small (below 3%), and the decrease of error becomes less significant, which demonstrates the effectiveness of the reconstruction algorithm. The threshold of 0.1 can be chosen as the compressive ratio mentioned in Section 2.5.2.

The reconstructed signal components in 1 of the 100 simulations are shown in Figure 4. We can observe that when adopting the compressive ratio of 0.1, both the smooth and sparse signal components can be reconstructed with high accuracy.

We also vary the magnitude of sparse signals to examine the reconstruction performance of the proposed method. The result and analysis are provided in Appendix H.

3.2. CSSD on a 2D Image

In this simulation study, we aim to decompose an image into the smooth background, sparse anomalies, and noise. A $350\times 350$ image with smooth background and the sparse anomaly is generated similar to Yan et al. (2017) (i.e., $\mathit{Y}=\mathit{M}+\mathit{A}+\mathit{E}$, where $\mathit{M}$ is the smooth background, $\mathit{A}$ is the sparse anomalies, and $\mathit{E}$ is i.i.d. Gaussian noise such that ${\mathit{E}}_i\sim \mathit{NID}\left(0,{\sigma }^2\right)$). The smooth background is generated from a linear combination of B-spline bases with $3\times 3$ knots, and the anomalies are generated from a sparse linear combination of B-spline bases with $88\times 88$ knots. The background and anomalies are shown in Figure 5.

We can see that the anomaly size covers a large range. The mean absolute value of the background plus anomaly is $\mu =0.21$. In this simulation, we first study the reconstruction performance under a noise-free scenario. Then, we increase the noise level to test the robustness of the algorithm.

3.2.1. Noise-Free Case.

In this section, we vary the compressive ratio from 4% to 100%. For each compressive ratio, we simulate 100 times and the average false-negative rate (FNR) and false-positive rate (FPR) are reported in Figure 6. The FPR is defined as the portion of normal pixels predicted as anomaly,

and the FNR is defined as the portion of anomalous pixels predicted as normal background,where $\mathit{A}$ is the true anomaly, $\widehat{\mathit{A}}$ is the predicted anomaly, and $I_{\text{Ω}}(\cdot )$ is the indicator function: that is,

From Figure 6, we can see that both the FPR and FNR ratios decrease as the compressive ratio increases. The FPR drops so fast that when the compressive ratio achieves 8%, there is no false alarm, which is desired for the anomaly detection algorithm. The FNR drops slower, and less than 5% of the anomaly pixels are ignored when the compressive ratio achieves 8%. However, it does not miss any cluster of the anomalies, even though some of them are small.

The recovered sparse signal components are shown in Figure 7 for compressive ratios 4% (Figure 7(a)), 8% (Figure 7(b)), 33% (Figure 7(c)), and 73% (Figure 7(d)). For comparison, we apply the SSD algorithm (Yan et al. 2017), of which the FNR and the FPR are zero and computation time is 0.16 seconds.

We record the computation time in Figure 8. A significant boosting of the computation is observed when the compressive ratio is 8%, which speeds up the SSD algorithm by 4.3 times.

3.2.2. Noisy Case.

The signal-to-noise ratio $\mu /ϵ \in \left[4,40\right]$ is studied to evaluate the robustness of the algorithm. The three-dimensional plot of FPR, the signal-to-noise ratio, and the compressive ratio are shown in Figure 9.

The purple line indicates the equipotential line of FPR = 0.01. We can see that with the increase of the signal-to-noise ratio, we are allowed to use a less compressive ratio to achieve satisfactory decomposition results. The majority of the area lies below the equipotential line of FPR = 0.01, which means that the proposed algorithm is robust.

4. Case Study

In this section, we use three real cases to demonstrate the effectiveness of the proposed CSSD/KronCSSD framework. For comparison, we also apply the SSD method in each case study. The compressive ratio ($c$) and the average computational time ($t$) for a single image are reported in Table 1. A significant transmission bandwidth reduction and computation boost can be observed with negligible performance degradation (Figures 10–12), compared with the vanilla SSD algorithm.

Table 1. Comparison Between KronCSSD and SSD

4.1. Surface Defect Detection in a Steel Rolling Process

As mentioned in Section 1, vision sensors collect high-resolution images of the product surface with a high data acquisition rate in the rolling processes. This poses a challenge in data storage, transmission, and processing. One sample image of size $128\times 512$ with typical anomalies is shown in Figure 10(a). The black scratches shown in the red block are surface anomalies. For a detailed description, the readers are encouraged to refer to Yan et al. (2018). The detected anomalies are shown in Figure 10, (b) and (c) by using the proposed KronCSSD algorithm and the SSD algorithm, respectively. The data set has 100 images, and more example results can be found in Appendix I.

The KronCSSD method achieves a similar anomaly detection performance with 54% bandwidth and is 1.8 times faster. By adopting the KronCSSD method, (i) we can achieve a faster anomaly detection and thus, reduce the loss through a timely intervention of the manufacturing process, and (ii) we can keep the manufacturing inspection information for a longer time period with the same storage capability, which is important for root cause analysis.

4.2. Solar Flare Detection

Another important application is solar flare detection from satellite images. A solar flare is defined as a sudden, transient, and intense variation in brightness over the sun’s surface. It has a significant influence on radio communication on Earth. Each second, thousands of high-resolution images are captured by a satellite, which poses a big challenge to real-time data transmission and processing (Yan et al. 2018). The data set has 300 images, and more example images can be found in Appendix I.

One sample image of size $232\times 292$ with a typical solar flare is shown in Figure 11(a), where the yellowish bright region is the solar flare. The detected anomalies are shown in Figure 11, (b) and (c) by using the proposed KronCSSD and SSD algorithms, respectively. The KronCSSD method achieves a similar solar flare detection performance as the SSD method but with 22% bandwidth, and it is 2.5 times faster. Notice that the decomposed images can also be used for downstream tasks, such as control charts and so on, which are beyond the scope of this paper.

By adopting the KronCSSD method, we can improve the transmission rate under the same transmission bandwidth and thus, achieve almost five times faster solar flare detection, which is of vital importance for protecting radio communications, power grids, and navigation systems.

4.3. Silicon Surface Indentation Detection

The stress map of size $90\times 550$ of a silicon surface laminate with surface indentation is shown in Figure 12, where clusters of high-stress areas indicate the surface indentation (Yan et al. 2017). We aim to detect those high-stress areas. The detected anomalies are shown in Figure 12, (b) and (c). We can see that the KronCSSD method achieves a similar detection performance with the SSD method but with 40% bandwidth, and it is 1.8 times faster, which significantly reduces the storage and transmission cost and improves the processing speed.

5. Conclusion

In this paper, we proposed a CSSD framework for efficient data acquisition, transmission, and processing for sparse anomaly detection in smooth backgrounds. To further enhance its computational efficiency, a KronCSSD framework is proposed for tensor data.

The contributions of this work are twofold. (i) Theoretically, we showed the feasibility of combining compressive sensing and smooth sparse decomposition. This enables the adoption of a compressive data acquisition approach. (ii) Practically, the proposed framework is compatible with many existing decomposition-based anomaly detection algorithms, such as SSD, ST-SSD, and so on, which achieve both a significant cost reduction in sensing, storage, and transmission and a boost in speed but with negligible loss in their performance.

In this article, we use a simulation study to demonstrate the effectiveness and robustness of the proposed CSSD/KronCSSD framework. Three case studies across different applications demonstrate the versatility of the proposed framework. The authors believe that the CSSD/KronCSSD framework can be applied in a wider range of applications toward more efficient data acquisition, transmission, and processing.

Further studies on the theoretical properties of the proposed KronCSSD can be a future direction.