SMDS-Net: Model Guided Spectral-Spatial Network for Hyperspectral Image Denoising

Fengchao Xiong, Member, IEEE, Jun Zhou, Senior Member, IEEE,  Shuyin Tao, Jianfeng Lu, Member, IEEE, Jiantao Zhou, Senior Member, IEEE and Yuntao Qian, Senior Member, IEEE
This work was supported in part by Jiangsu Provincial Natural Science Foundation of China under Grant BK20200466, and the National Natural Science Foundation of China under Grant 62002169, 61905114, 61971476 and 62071421, the National Key Research and Development Program of China under Grant  2018YFB0505000 and 2018AAA0100500, Macao Young Scholars Program (No. AM2021020) and Macau Science and Technology Development Fund under 077/2018/A2. (Corresponding author: Fengchao Xiong.)
Fengchao Xiong is  with the School of Computer Science and Engineering, Nanjing University of Science and Technology, Nanjing 210094,  China and also  with the State Key Laboratory of Internet of Things for Smart City, Department of Computer and Information Science, University of Macau,  Macau 999078, China (e-mail: fcxiong@njust.edu.cn).
Jun Zhou is with the School of Information and Communication Technology, Griffith University, Nathan, Australia (e-mail: jun.zhou@griffith.edu.au).
Shuyin Tao and Jianfeng Lu are with the School of Computer Science and Engineering, Nanjing University of Science and Technology, Nanjing 210094, China (e-mail<EMAIL_ADDRESS>lujf@njust.edu.cn).
Jiantao Zhou is with the State Key Laboratory of Internet of Things for
Smart City, Department of Computer and Information Science, University of
Macau, Macau 999078, China (e-mail: jtzhou@um.edu.mo).

Yuntao Qian  is with the College of Computer Science, Zhejiang University, Hangzhou 310027, China (e-mail: ytqian@zju.edu.cn).

Deep learning (DL) based hyperspectral images (HSIs) denoising approaches directly learn the nonlinear mapping between  noisy  and  clean HSI pairs.  They normally do not consider the physical characteristics of HSIs. This drawback makes the models  lack  interpretability that is key to understand their denoising mechanism and also limits their denoising ability.   To tackle this problem, we introduce a novel model guided interpretable network for HSI denoising.  Fully considering the spatial redundancy, spectral low-rankness and spectral-spatial correlations of HSIs, we first establish a subspace based multidimensional sparse (SMDS) model under the umbrella of tensor notation.    After that, the model is unfolded into an end-to-end network named  as  SMDS-Net whose fundamental modules are seamlessly connected with the denoising procedure and optimization of  the SMDS model. This makes SMDS-Net convey clear physical meanings, i.e., learning the low-rankness and sparsity of HSIs.  Finally, all key variables are obtained by an end-to-end training. Extensive experiments and comprehensive analysis  on  synthetic and real-world HSIs confirm the strong denoising ability, learning capability and high interpretability of SMDS-Net   against the state-of-the-art HSI denoising methods.

Hyperspectral image  denoising, model-based neural network, low-rank representation, multidimensional sparse representation.

§ INTRODUCTION

Hyperspectral images (HSIs) have enabled many practical applications such as pedestrian detection <cit.>, object tracking <cit.>, medical diagnosis <cit.> and more <cit.> thanks to its material identification ability provided by numerous light-wavelength indexed bands. Because of narrow bands, limited sensitivity of sensors or their defects, and the interference of the imaging environment, the captured HSIs tend to be corrupted by noises. The degradation of HSIs brings many drawbacks in the exploitation of HSIs and significantly decreases their utility as well. Therefore, HSI denoising is often needed as a critical preprocessing step to enhance the quality of HSIs.

[NGMeet <cit.>]
[QRNN3D <cit.>]
Denoising results on Real-world Urban HSI with respect to bands 1, 109, and 208. Our SMDS-Net obtains the best qualitative  results. (Best view on screen with zoom)

HSI denoising  methods can be divided into two  categories, i.e., model-based and learning-based. Model-based methods make explicit hypotheses on the underlying physical characteristics of HSIs and introduce various prior structures such as spectral low-rankness <cit.>, spatial redundancy <cit.> and nonlocal similarity <cit.>.    These model-based methods are  physically interpretable   but heavily depend on the effectiveness of hand-crafted  priors, computational demanding iterative optimization, and exhaustive hyperparameter tuning.

Illustration of the proposed method.  The observed noisy HSI $\Ycal$ is first projected into a low-dimensional image $\Gcal$ with respect to learned subspace $\A$ to depict the spectral low-rankness.  A number of overlapping cubes are then extracted from $\Gcal$ and sparsely encoded by a number of multidimensional sparsity coding blocks and a cube reconstruction block to model the spatial redundancy. After that, All the resulted cubes are aggregated and reconstructed with $\A$ to yield  the denoised image $\widehat{\Xcal}$.

As an alternative paradigm, deep learning (DL)-based methods directly learn the nonlinear mapping between noisy-clean HSI pairs for  denoising <cit.> rather than making subjective physical assumptions on the underlying HSIs.  Attributing to their strong representation ability, DL-based methods provide competitive performance when trained with a large number of noisy-clean HSI pairs.  However, the high cost of acquiring hyperspectral data can dramatically limit their generalization ability and make the learning infeasible. Moreover,  complicated network architectures also hinder deep analysis of the functionality of different modules and  comprehensive understanding of the rationality behind their denoising mechanisms.

Given the above discussion, it is preferable to design a  HSI denoising method that combines the advantages of both model-based and  DL-based methods, i.e., of high interpretability while supporting discriminative learning.  The consideration of spectral-spatial correlation  facilities information sharing among bands and improves denoising performance <cit.>. Spectral low-rankness and spatial redundancy priors are beneficial to better model underlying clean HSI, helping more effectively remove noises <cit.>.  Therefore, to reach the above goal, the spectral-spatial correlation, spectral low-rankness, and spatial redundancy structure priors and end-to-end learning are very important factors to be considered.

In this paper, a model guided network is introduced for HSI denoising,  which simultaneously considers all the above factors in a unified network whose framework is shown in Fig. <ref>.  Specifically,  we first build a  subspace based multidimensional sparse (SMDS) model under the umbrella of multidimensional tensor sparse representation, which projects the high-dimensional HSIs into a low-dimensional spectral subspace and exploits the spatial redundancy property on the projected image.  SMDS has higher capability of removing noises because it considers both spectral and spatial structure powered by tensors. After that, an end-to-end deep neural network named SMDS-Net is obtained by unfolding the denoising procedure and optimization algorithm of SMDS. SMDS-Net includes five stages, i.e. spectral subspace projection, cube extraction, spatial multidimensional sparse modeling, cube aggregation, and image reconstruction. Each stage has a clear physical meaning, implying the interpretability of our network.

To the best of our knowledge, this is the first time that low-rankness, sparsity, and spectral-spatial structural priors were simultaneously taken into consideration in a unified end-to-end network for HSI denoising. The main contribution of this paper can be summarized as follows:

* We propose a novel subspace based multidimensional sparse  model which captures  the spectral-spatial correlations, spectral low-rankness and spatial sparsity  of HSIs under the umbrella of tensors, making the model highly capable of  removing noises.

* We construct an  end-to-end deep neural network by unfolding the denoising procedure and optimization algorithm of SMDS.  SMDS-Net combines the hybrid advantages of  model-based and DL-based methods,  such as clear interpretability, higher representation capability, and strong learning ability.

* Extensive experiments confirm the advantages of the SMDS-Net over other state-of-the-art model-based and DL-based methods.   Moreover, a comprehensive analysis of the network is carried out to fully explain its denoising mechanism and learning ability.

This paper is organized as follows: Section <ref> briefly reviews recent works on model-based and DL-based HSI denoising as well as deep unfolding.  In Section <ref>, we  present the subspace based multidimensional sparse  model and its optimization procedure. Section <ref> shows how to unfold the SMDS model into an end-to-end neural networks as well as the implementation details.   Section <ref> reports the experimental results on both synthetic and real-world data and  compare them against several other competing approaches.  This paper concludes in  Section <ref>.

§ RELATED WORK

This section briefly introduces the recent works concerning model-based and deep learning-based HSI denoising methods and also deep unfolding.

§.§ Model Based Approaches

HSIs can be considered as a stack of 2D images along the spectral dimension. Therefore, traditional 2D image denoising methods such as  Non-local means (NLM) <cit.>, BM3D <cit.> and low-rank/sparse representation <cit.> can be directly applied for HSI denoising in a band-wise manner. Band-wise methods ignore the spectral-spatial correlation in HSIs, leading to unsatisfactory denoising performance.  Alternatively, spectral-spatial methods  jointly leverage the information in different bands and simultaneously model HSIs in both spectral and spatial domains, thereby achieve more promising performance. Typical examples include BM4D <cit.> and 3D NLM <cit.>, which model HSIs at 3D cube level rather than patch level.

Following this line,  sparse representation based methods   encode the observed HSIs using a few atoms in predefined fixed  dictionaries <cit.> or data-driven adaptive dictionaries <cit.>.  Besides sparsity, the spectral correlation among bands and spatial correlation among  pixels also imply there exist low-dimensional structures  in HSIs that can be depicted by matrix rank minimization <cit.> or low-rank decomposition <cit.>.  Moreover, structural sparse/low-rank representation also demonstrate their effectiveness in HSI denoising by taking the local and nonlocal similarity structures of HSIs into sparse/low-rank models <cit.>. Matrix-based sparse/low-rank models inevitably destroy the spectral-spatial structure of HSIs when converting 3D HSIs into a 2D matrix and require more parameters.  For this reason, tensor-based sparse/low-rank representations are introduced and successfully employed in HSI denoising <cit.>.

Although these methods work well when the prior model assumptions match well with the HSIs to be processed, they normally require iterative numerical optimization which is computationally expensive. Exhaustive hyperparameter tuning also limits their  effectiveness in complex scenarios.

§.§ Deep Learning Based Approaches

Instead of subjective physical assumptions on the underlying HSIs, DL-based methods leverage their strong image representation capability and learn a direct mapping function between noisy and clean HSIs <cit.>.  Chang et al. <cit.> learned a series of multiple channels of 2D filters for the consideration of both spatial and spectral structures of HSIs. Yuan et al. <cit.> proposed multiscale spectral-spatial denoising through multi-scale feature extraction and multilevel feature representation using deep neural networks. To effectively extract the joint spectral-spatial information, Zhao et al. <cit.> and Maffei et al. <cit.> took the current band and its $K$ adjacent bands as input. Alternatively, structural spectral-spatial correlations in HSIs can also be extracted by 3D convolutions <cit.>. Instead of directly operating on the original HSIs, noisy HSIs can be converted into low-dimensional HSIs before applying CNNs on the reduced images <cit.>.

Although these methods are able to recover the underlying image via training deep models on a large number of noisy-clean HSI pairs, they often ignore very important intrinsic prior structures of  HSIs such as sparsity and low-rankness as mentioned previously.  Instead, they perform denoising in a “black box" mechanism whose architecture is empirically determined in a heuristic manner and exhausting trials, making further model improvement difficult. Therefore, the interpretability of the actual functionality is a matter of concern.

§.§ Deep Unfolding

Deep unfolding starts from the iterative algorithm deduced by the model-based methods and maps each iteration into a typical layer of a deep neural network. By stacking a predefined number of layers, a hierarchical deep network architecture can be obtained. Deep unfolding provides an attractive way to inherit the advantages of model-based methods, for example, high interpretability and good generalization capability and absorb the merits of learning-based methods such as strong learning ability and computational efficiency.  Attributing to this, deep unfolding networks have been successfully applied in many tasks, such as hyperspectral fusion <cit.>, rain removal <cit.>, image super-resolution <cit.> and  deblurring <cit.>. Our work falls in the line of deep unfolding and converts the subspace based multidimensional sparse model into an end-to-end network to achieve more  more effective  HSI denoising.

§ SMDS HSI DENOISING MODEL

In this section, we introduce our subspace based multidimensional sparse model and also its associated optimization procedure.

§.§ Notations

In this paper, we represent a scalar by a lowercase letter, e.g. $y$; a vector by a   lowercase letter, e.g. $\y \in \Rbb^{I}$; a matrix by  a  capital letter, e.g. $\Y \in \Rbb^{I_1 \times I_2}$; and a tensor by  a  calligraphic letter, e.g. $\Ycal \in \Rbb^{I_1 \times I_2 \times \ldots \times I_N}$.  The $L_0$ norm of a tensor $\Ycal$ is defined as the number of non-zero entries, i.e, $\|\Ycal\|_0=\#\{\Ycal_{i_1,i_2,\cdots, i_N}\neq 0\}$.  The $L_1$ norm and Frobenius norm of $\Ycal$ are respectively defined as $\|\Ycal\|_1=\sum_{i_1}\sum_{i_2}\cdots\sum_{i_N}|\Ycal_{i_1,i_2,\cdots, i_N}|$,  $\|\Ycal\|_F=\sqrt{\sum_{i_1}\sum_{i_2}\cdots\sum_{i_N}(\Ycal_{i_1,i_2,\cdots, i_N})^2}$.  The n-mode unfolding  is the process of  arranging all the n-mode vectors as columns of a matrix, denoted by $\Y_{(n)}$.  The n-mode product of a tensor $\Ycal$ and matrix $\U \in \Rbb^{J\times I_n}$ is written as $\Ycal\times_n\U$ and computed by $(\Ycal\times_n\U)_{i_1\cdots i_{n-1}j i_{n+1}\cdots i_N}=\sum_{i_n=1}^{I_n}y_{i_1i_2\cdots i_N} u_{ji_n}$.

§.§ Model Formulation

Given an HSI $\Ycal$ with $H \times W$ pixels and $B$ bands corrupted by additive zero-mean Gaussian noise $\Ncal$, the observation model can be represented as

\begin{equation}
\Ycal=\Xcal+\Ncal \label{eq:model}
\end{equation}

where $\Xcal$ is the clean HSI to be estimated. This  is an ill-posed inverse problem and the prior knowledge on $\Xcal$ is very important to solve the problem. Specifically, we consider spectral-spatial correlation,  spatial sparsity and spectral low-rankness priors to model $\Xcal$.

Spectral low-rankness:  As stated in <cit.>, the spectral correlations among bands make acquired spectra in fact lie in a low-dimensional subspace, i.e.,

\begin{equation}
\Xcal=\Gcal\times_3\A \label{eq:spectral}
\end{equation}

Here $\A=[\a_1,\cdots,\a_r] \in \Rbb^{B\times R}$ is the spectral subspace spanned by $R$  bases with $R$ far less than $B$ and $\Gcal \in H \times W \times R$ is the representation tensor with respect to $\A$. In general, $\A$ can be obtained by many approaches, for example, unmixing and singular value decomposition (SVD) on the unfolded $\Ycal$ along spectral models. As orthogonal subspace shares many merits such as disentangling the correlation among bands <cit.> and also simplifying the optimization, we choose a set of orthonormal bases, i.e.,  $\A^{\Trm}\A=\I$.

Spatial sparsity: As the low-rank projection is performed in the spectral domain, the projected image can be considers  as a linear combinations of the bands of $\Xcal$, meaning that the pixels nearby shares similarity.  On the other hand, the spectral low-rankness is a global  transformation for all the pixels, indicating that for pixels in a local region, there still exist dependencies.   Instead of separately encoding each mode (channel) of $\Gcal$, we  model the local 3D cubes within $\Gcal$ under the framework of tensor based multidimensional (MD) sparse  representation.  In this way,  all the spatial structures in each mode remain  and  can be shared  across channels with less parameters, facilitating better handling  of noises. This can be implemented  by dividing $\Gcal$ into a number of small  overlapping cubes, calculating their sparsity, and putting them back to the original position using their denoised estimation by averaging with their overlapped cubes. Specifically, let $\Rcal_i$ denote an operator extracting the $i$-th cube from $\Gcal$, under the framework of  MD sparse representation,  the sparse representation problem of each cube $\Gcal_i=\Rcal_i\Gcal$ is formulated as

\begin{equation}
\arg \min_{\Bcal_i}\frac{1}{2}\|\Gcal_i-\Bcal_i\times_1\D_1 \times_2\D_2 \times_3\D_3\|_F^2, \text{s.t.}\: \|\Bcal_i\|_0 \leq K \label{eq:spatial}
\end{equation}
where $\D_j \in \Rbb^{I_j \times M_j } (j=1,2,3)$ is the overcomplete dictionary matrix along $j-$th mode with $I_j \leq M_j$. Unlike Tucker decomposition which only models the space that the tensor itself belongs to, the model in Eq. (<ref>) can simultaneously and adaptively characterize redundant basis of the space where multiple tensors reside in a data-driven way.

Substituting Eq. (<ref>) and  Eq. (<ref>) into Eq. (<ref>), the proposed subspace based MD sparse representation denoising method is written as

\begin{equation}
\begin{split}
&\!\{\A^*, \Gcal^*, \Bcal_{i}^*\}=\arg\min_{\A,\Gcal, \Bcal_i} \frac{1}{2}\|\Gcal\times_3\A-\Ycal\|_F^2\\
&+\mu\sum_{i}\left(\frac{1}{2}\|\Gcal_i-\Bcal_i\times_1\D_1 \times_2\D_2 \times_3\D_3\|_F^2+\lambda\|\Bcal_i\|_0\right)\\
& \text{s.t.}\:\A^{\Trm}{\A}=\I
\end{split}\label{eq:obj}
\end{equation}

where $\D_j$ is the dictionary along each dimension. In Eq. (<ref>), the first term is the data fidelity term that demands proximity between the observed noisy $\Ycal$ and estimated image $\Gcal\times_3\A$. The second term stands for the prior assumption on $\Gcal$ that every cube $\Gcal_i$ has an MD sparse representation and $\lambda>0$ controls sparse regularization level. $\mu$ balances the data fidelity term and sparse representation term and can be implicitly set to 1 as in <cit.>.

§.§ Model Optimization

SMDS for HSI denoising.
[1]
Noisy HSI $\Ycal$.

Denoised HSI $\Xcal$.

Estimate $R$ using Hysime <cit.>.

Conduct SVD on  $\Y_{(3)}$, i.e., $\Y_{(3)}=\U\S\V^{\Trm}$ and set $\A=[\v_1,\cdots, \v_R]$.

Project $\Ycal$ to obtain $\Gcal$ with Eq. (<ref>).//spectral subspace projection

Extract overlapping cubes from $\Gcal$. //cube extraction
Compute $\Bcal_i$ with Eq. (<ref>). //spatial MD sparse coding
TISTA converges.
Obtain $\widehat{\Gcal}$ with Eq. (<ref>). //cube reconstruction and aggregation

Output the estimated $\widehat{\Xcal}$ with Eq. (<ref>). //image reconstruction

Eq. (<ref>) is difficult to solve because of the complicated constraint on $\A$ and $\Bcal_i$. Similar to <cit.>, we solve Eq. (<ref>) in a greedy way, which includes sequentially learning  the subspace $\A$ from $\Ycal$ and  MD sparse representation $\Bcal_{i}^*$ from $\Gcal_i$.

Learning $\A$: The optimization of $\A$ can be considered as a Tucker-1 decomposition problem. Under the framework of higher-order singular value decomposition (HOSVD), we can obtain $\A$ by performing singular decomposition on the unfolded matrix $\Y_{(3)}$ along the spectral dimension, i.e., $\Y_{(3)}=\U\S\V^{\Trm}$. Then, $\A$ is assigned with $\A=[\v_1,\cdots, \v_R]$ where $R$ is estimated by the Hysime algorithm <cit.>. In order to more accurately estimate $R$ and $\A$, some fast denoising methods such as NLM <cit.> can be applied. Once $\A$ is obtained, $\Gcal$ can be obtained by

\begin{equation}
\Gcal=\Ycal\times_3\A^{\Trm} \label{eq:projG}
\end{equation}

Learning $\Bcal_{i}$: The $L_0$ norm of $\Bcal$ is intractable because of its NP-hard property. It can be relaxed to the $L_1$ norm to result in  a convex optimization problem. Along this line, the optimization problem of $\Bcal_i$ can be formulated by

\begin{equation}
\arg \min_{\Bcal_i}\frac{1}{2}\|\Gcal_i-\Bcal_i\times_1\D_1 \times_2\D_2 \times_3\D_3\|_F^2+\lambda\|\Bcal_i\|_1
\end{equation}
This problem can be solved with Tensor-based Iterative Shrinkage Thresholding Algorithm (TISTA) <cit.> with the following solution:

\begin{equation}
\begin{split}
\Bcal_i^{(k+1)}&\!\!=\!P_{\lambda/L}(\Bcal_i^{(k)}\!\!-\!\!\frac{1}{L}(\Bcal_i^{(k)}\!\times_1\D_1^{\Trm}\D_1\!\times_2\D_2^{\Trm}\D_2\!\times_3\D_3^{\Trm}\D_3\\
&-\Gcal_i\times_1 \D_1^{\Trm}\times_2\D_2^{\Trm}\times_3\D_3^{\Trm}))
\end{split} \label{eq:tista}
\end{equation}
where $L$ is a Lipschitz constant, $P_{\lambda/L}=\text{sgn}(x)(|x|-\lambda/L)_+$ performs soft-thresholding operation to achieve the sparsity of $\Bcal_i^{(k+1)}$ and $k$ indexes the iteration number.

After  obtaining $\Bcal_i$,  $\widehat{\Gcal}$ can be obtained by averaging $m$ estimates for each entry in the cube, i.e.,

\begin{equation}
\widehat{\Gcal}=\frac{1}{m}\sum \limits_{i}\Rcal_{i}^{'}\Bcal_i\times_1\D_1 \times_2\D_2 \times_3\D_3 \label{eq:gupdate}
\end{equation}

where $\Rcal_i^{'}$ is a linear operator  placing back the cube  $\Bcal_i\times_1\D_1 \times_2\D_2 \times_3\D_3$ at the position centered on entry $i$.  Finally, the output denoised HSI $\widehat{\Xcal}$ can be obtained  by Eq. (<ref>).  Algorithm <ref> summarizes the denoising procedure of  the SMDS method.  SMDS simultaneously considers the spectral-spatial correlation, spectral low-rankness and spatial sparsity structures of HSIs and has strong denoising ability. However, SMDS formulates denoising as an optimization problem, which normally requires computationally inefficient iterations and exhausting hyperparameter tuning, making the denoising efficiency and effectiveness still limited. In the next section, we will introduce an unfolded SMDS network to solve this problem.

§ UNFOLDED SMDS NETWORK

In this section, we unfold all the steps of SMDS in Algorithm <ref> as network layers to yield an end-to-end neural network for discriminative training.

§.§ Network Architecture

Algorithm <ref> can be decomposed into five stages, which corresponds to five parts in the network shown in Fig. <ref>.  Here, we describe the details of key parts as follows:

Spectral subspace projection: The first stage of the network is mapping $\Ycal$ into a low-dimensional subspace held by $\A$ following Eq. (<ref>).  Since different HSIs may reside in different subspaces, it is unrealistic to encode $\A$ as a network parameter to learn the unique subspace for all the HSIs.    For this reason,  $\A$ is precomputed and taken as network input along with  $\Ycal$. Consequently, our network can not only  employ the internal data, i.e., the noisy HSI itself for adaptive learning of spectral subspace but also  the extern data, i.e., the training set to learn the shared dictionaries  for all the HSIs.  Moreover, with the spectral subspace as input, our network can be used for input HSI with an arbitrary number of bands, offering much flexibility.

MD sparse coding: After cube extraction, MD sparse coding stage learns the sparse coding of each cube with respect to multidimensional dictionaries, corresponding to Eq. (<ref>). Eq. (<ref>) can be  decomposed into the following three steps:

\begin{align}
&\Ecal_i^{(k+1)}=\Gcal_i-\Bcal_i^{(k)}\times_1\D_1 \times_2\D_2 \times_3\D_3\label{eq:objE}\\
\end{align}
where $\LAmbda^{(k)} \in \Rbb^{M_1\times M_2\times M_3}$ denotes the thresholding parameter tensor  for $\Hcal_i^{(k)}$ in the $k$-th layer of MD sparse coding module. $1/L$ is absorbed into $\LAmbda^{(k)}$ and $\D_j$. Following <cit.>, we  decouple $\C$ from  $\D$  and set different $\LAmbda$ in each layer for accelerated convergence.

Eq. (<ref>) and  Eq. (<ref>)  are similar in function, mapping the processed tensor to a new dimension along three modes. They can  be regarded as three linear transformation layers.  $P_{\LAmbda_{(k)}}$ in Eq. (<ref>) corresponds to the nonlinear function allowing for the sparsity of $\Hcal^{(k+1)}$, which can be implemented with $\text{Relu}(x-\theta)-\text{Relu}(-x-\theta)$. Eqs. (<ref>), (<ref>) and (<ref>) make up  a typical MD sparse coding block in  Fig. <ref>. Stacking this block  for $K$ times  is equivalent to running Eqs. (<ref>), (<ref>) and (<ref>) for $K$ times, resulting in a subnetwork for MD sparse coding.

Cube reconstruction and aggregation: Similarly, Eq. (<ref>) can  be interpreted as three linear transformation layers followed by an averaging layer.  As in <cit.>, we decouple $\W$ from $\D$ for more effective end-to-end training, i.e., $\widehat{\Gcal_i}=\Bcal_i^{(K)}\times_1\W_1 \times_2\W_2 \times_3\W_3$. This corresponds to the cube reconstruction block in Fig. <ref>. After that the overlapped cubes are aggregated, i.e., $\widehat{\Gcal}=\frac{1}{m}\sum \Rcal_i^{'}\widehat{\Gcal_i}$.

Image reconstruction:  The last image reconstruction stage uses Eq. (<ref>) to produce the denoised HSIs $\widehat{\Xcal}$.

Comparison of different methods on 50 testing HSIs of ICVL  Dataset. The top two values are marked Red and Blue.

4c|Sparse methods   5c|Low-rank  methods   5cDL  methods

2*$\sigma$   2*Index   Noisy   BM3D   BM4D   TDL   MTSNMF   PARAFAC   LLRT   NGMeet   LRMR   LRTDTV   Dn-CNN   HSI-SDe   HSID-   QRNN3D   SMDS-Net

<cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>   CNN <cit.>   CNN <cit.>    <cit.>   (Ours)

3*[0-15]   PSNR   33.47   45.68    45.12    38.84   45.62   39.09    46.17    39.42    40.48   42.17      41.92   41.06   38.08   43.36      47.81

SSIM   0.6295   0.9786   0.9753    0.8481    0.9562   0.9468   0.9649   0.8647   0.9601   0.9717      0.9596   0.9565   0.9669   0.9884   	0.9899

SAM   0.1622   0.0252    0.0254    0.0829    0.0409    0.0299   0.0319   0.0891   0.0177   0.0229      0.0274   0.0329   0.0511   0.0336   0.0139

3*[0-55]   PSNR   21.67   38.67   38.21   29.45   37.57   35.71   37.19   30.48   36.20   39.69   37.90   35.66   32.55   40.15   42.00

SSIM   0.2361   0.9363   0.9217   0.5238   0.8569   0.8808   0.8389   0.6748   0.9295   0.9625   0.9293   0.8736   0.8421   0.9729   0.9733

SAM   0.5014   0.0541   0.0552   0.2409   0.1372   0.0542   0.1020   0.2830   0.0334   0.0335   0.0504   0.0602   0.0918
0.0382   0.0242

3*[0-95]   PSNR   16.97   35.99    35.27    25.40   34.20   32.80    32.49    27.20   34.27    38.12    34.65   32.45   29.12   37.66   39.34

SSIM    0.1442   0.9080   0.8764     0.3735    0.7998   0.7810   0.7191    0.5506   0.9153   0.9539   0.8442   0.7875   0.7049   	0.9479   			0.9593

SAM   0.7199    0.0738   0.0799   0.3641   0.2128   0.0997   0.1687    0.4275   0.0484   0.0402   	0.1094   0.0876   	0.1316   	0.0468   0.0303

§.§ Network Training

Training loss:
The training loss for a given training set of noisy-clean pairs is defined as the Euclidean distance  between the output of SMDS-Net  and the ground truth $\Xcal$, i.e.,

\begin{equation}
L=\|\text{SMDS-Net}(\Ycal, \Theta)-\Xcal\|_F^2
\end{equation}

where $\Theta=\{(\C_j,\D_j,\W_j)_{j=1,2,3}, (\LAmbda^{(k)})_{k=1,\cdots,K}\}$ represents all parameters to be learned.

Implement details: SMDS-Net is implemented on the Pytorch platform and trained on  two  NVIDIA GeForce RTX 3090 GPUs with 300 epochs. We adopt the Adam optimizer with the batch size of 2 and  patch size of $56\times 56$.  The initial learning rate is $5\times10^{-3}$ and multiplied by 0.35 for every 80 epochs.  The cube size $[I_1, I_2, I_3]$ is set as $[9,9,9]$.   The dictionaries were initialized with discrete cosine transform (DCT) basis whose number  $[M_1, M_2, M_3] $ are given by [9, 9, 9], leading to dictionaries size of $9\times 9$ along three modes. The number of unfolding $K$ is set to 6.

Training dataset:  The same as <cit.>, 100 images containing $1392\times 1300$ pixels and 31 spectral bands are selected from ICVL hyperspectral dataset [http://icvl.cs.bgu.ac.il/hyperspectral/] to construct the training set.   Data augmentation is performed  to enlarge the training  set including random flipping, cropping and resizing. After that  images with a size of $56\times 56\times 31 $ are used to train the network.

Denoising results on the  gavyam_0823-0933 image with the noise variance in [0-95]. The false-color images were generated by combining bands 5, 18, 25. SMDS-Net achieves the best visual results with less artifacts. (Best view on screen with zoom)

Denoising results of pixel (451, 135) in  gavyam_0823-0933 HSI with the noise variance in [0-95].

§ EXPERIMENTAL RESULTS

In this section, we conduct experiments on both close-range HSIs and also remote sensing HSIs to show the denoising ability of our methods.  Moreover, comprehensive studies are conducted to illustrate its high learning ability.

§.§ Experimental Settings

Methods of Comparison:  The competing denoising methods include: four sparse representation based methods, i.e., BM3D <cit.> [https://webpages.tuni.fi/foi/GCF-BM3D/], BM4D <cit.> [http://www.cs.tut.fi/foi/], TDL <cit.> [http://gr.xjtu.edu.cn/web/dymeng/], MTSNMF <cit.>, five low-rank representation  based methods, PARAFAC <cit.>, LLRT <cit.> [https://owuchangyuo.github.io/], NGMeet <cit.> [github.com/quanmingyao/NGMeet], LRMR <cit.> [https://sites.google.com/site/rshewei/home], LRTDTV <cit.> [http://gr.xjtu.edu.cn/en/web/dymeng/3] and four deep learning based method, i.e., DnCNN <cit.> [https://github.com/cszn/DnCNN], HSI-SDeCNN <cit.> [https://github.com/mhaut/HSI-SDeCNN], HSID-CNN <cit.> [https://github.com/qzhang95/HSID-CNN], QRNN3D <cit.> [https://github.com/Vandermode/QRNN3D]. BM3D and DnCNN are executed in a band by band  manner.

Noise Setting: Due to the sensitivity of the sensor at different wavelengths, the illumination condition, and other factors, the noise levels in different bands are usually different <cit.>. That is to say, there exist nonindependent and identically distributed (non-i.i.d.) noises in the spectral domain,  see Fig. <ref>(a) for an example on the Urban HSI. The noise distribution in bands 1, 109, and 208 show obvious differences. For this reason, following <cit.>, the synthetic data are accordingly generated by manually adding different noise variances ($\sigma^2$) of additive white Gaussian noise (AWGN)  to each band of clean HSI.   Specifically $\sigma$ is chosen in three ranges, i.e, [0-15], [0-55] and  [0-95].  In these three cases, the accurate noise variances for all the bands  are provided for testing.  Except for BM3D and DnCNN, all other competing methods which require the noise level as input are provided with the averaged noise level.

Evaluation Measures: PSNR, SSIM, and SAM are used for quantitative assessment. SAM depicts the spectral differences between clean and denoised HSIs. Smaller SAM and larger PSNR and SSIM imply better denoising. Moreover, we compare the restored HSI with clean HSI and noisy HSI to illustrate their differences for qualitative analysis.

As in the training, the input HSIs are cropped into multiple overlapping full-band cubes containing $56\times 56$ pixels, with a stride of 12, to test our SMDS-Net. After that, the denoising results of all the cubes are summarized to produce estimated HSI.

§.§ Experiments with ICVL Dataset

§.§.§ Quantitative Comparison
As in <cit.>, we use 50 images in ICVL  whose main region size of  $512\times 512\times 31$ are cropped to quantitatively evaluate the denoising performance of all the methods.     The denoising results are shown in Table <ref>, in which the best two results are highlighted in red and blue.  As can be seen, the proposed SMDS-Net achieves the state-of-the-art results and provides the best denoising performance in all the cases, confirming its high denoising capacity.   Among model-based methods,  TDL, PARAFAC, LLRT, and NGMeet are less effective, especially in the case where the noise levels largely differ between bands, for example, $\sigma \in [0-95]$.  The main reason is that all of them inherently assume the noise distributions are identical over all bands, which is far from the requirements of real-world HSI denoising that the noises are spectrally non-i.i.d. Additionally, the spectral non-i.i.d noises make it more difficult to select parameters that are suitable for all the bands, even after intensive parameter tuning.  This however is not the case with the proposed SMDS-Net as it learns all the parameters from a large scale of data and treats the denoising as a feedforward process without parameter tuning, making it more capable of removing complex noises.

Compared with MTSNMF and LRMR,  LRTDTV and SMDS-Net achieve more satisfactory results as the former two methods treat HSIs as a matrix in the process of denoising, unavoidably losing the spatial structural information contained in HSIs.  QRNN3D and SMDS-Net dominate alternative DL-based methods, i.e., HSI-SDeCNN, HSID-CNN, and band-wise Dn-CNN because of better incorporation of spectral-spatial structural correlation and global spectral correlation.  The hybrid advantages of the low-rank multidimensional sparse model and the strong learning ability of DL help our SMDS-Net outperform QRNN3D, especially concerning PSNR and SAM index. This experiment evidently shows the benefit of developing DL-based denoising methods considering the physical model of HSIs.

Comparison of different methods on  32 testing HSIs of CAVE dataset. The top two values are marked Red and Blue.

4c|Sparse methods   5c|Low-rank  methods   5cDL  methods

2*$\sigma$   2*Index   Noisy   BM3D   BM4D   TDL   MTSNMF   PARAFAC   LLRT   NGMeet   LRMR   LRTDTV   Dn-CNN   HSI-SDe   HSID-   QRNN3D   SMDS-Net

<cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>    <cit.>   CNN <cit.>   CNN <cit.>    <cit.>   (Ours)

3*[0-15]   PSNR   33.42   45.01   44.78   37.84    42.79   34.91    45.02   37.81    35.95   36.88    40.10    41.23   35.97     40.10    45.40

SSIM   0.6256    0.9753   0.9744   0.8248   	0.9688   	0.8678   	0.9634   0.8619   0.9425   0.9268    0.9692   	0.9593   0.9186   	0.9692   0.9826

SAM   	0.6113    0.1143   0.1325   0.4348   0.1388   0.2174   0.1773   0.4259   0.1286   0.1835   	0.1767   0.1476   0.2646   0.1767   0.1023

3*[0-55]   PSNR   21.51    37.87    38.04   28.52   35.63   32.87        36.34   28.22   32.74    34.49    37.11    36.10    33.06     37.95    39.28

SSIM   0.2389   0.9149   0.9102   0.4994   0.8413   0.8124   0.8244   0.5788   0.9058   0.9015   0.8510   0.8830   0.8499   0.9481   0.9510

SAM   1.0008   0.2224   0.2936   0.7630   0.3764   0.3324   0.4446   0.8359   0.1975   0.2917   0.2374   0.2234   0.3951   0.2242   0.1671

3*[0-95]   PSNR      17.55   35.39   34.93   24.46   32.57   30.61   31.59   24.47   31.17   33.92   34.23   33.30   30.16   35.94    36.22

SSIM    0.1632   0.8707   0.8330   0.3546   0.7677   0.7124   0.6868   0.4367   0.8780   0.8939   0.8213   0.8061   0.7374   0.9153   0.9247

SAM   1.1417   0.2855   0.4305   0.9139   0.5191   0.4814   0.6103   0.9619   0.2724   0.3051   0.3097   0.2694   0.5437   0.2873   0.1977

Denoising results on the  chart and stuffed toy HSI with the noise variance in [0-95]. The false-color images were generated by combining bands 5, 18, 25. SMDS-Net achieves the best visual results with less artifacts. (Best view on screen with zoom)

Denoising results of pixel (200, 100) in  chart and stuffed toy HSI.

§.§.§ Visual Comparison

To thoroughly show the denoising effectiveness of our method, we present the visual denoising results of all the methods on gavyam_0823-0933 image in Fig. <ref>, where the false-color images are generated using bands 5, 18, 25.   Compared with other methods, band-wise BM3D and DnCNN lose some structures because of the ignorance of the spectral-spatial correlation among bands. The advantages of tensor in spectral-spatial modeling and the consideration of spectral low-rankness help PARAFAC and LRTDTV   provide more appealing visual results.  From highlighted areas, we can see LRMR, MTSNMF, BM4D, HSIDCNN, and HSI-SDeCNN can not remove all the noises. QRNN3D removes noises but introduces impulse noises in the denoised image.  Thanks to the strong physical model that considers the spatial and spectral properties of HSIs and the learning ability of the unfolded network,  SMDS-Net produces the best visual quality.

Fig. <ref> illustrates the spectral signature of pixel (453,135) in   gavyam_0823-0933  HSIs before and after denoising. As a intuitive comparison shown in  Fig. <ref> (a) and (b), there are obvious spectral distortion because of the
the contamination of heavy noises.  Among all the compared methods, our SMDS-Net can more capable of recovering the spectrum by best matching with the reference one.   Overall, the qualitive comparisons in Fig. <ref> and Fig. <ref> further verifies strong denoising capacity of proposed SMDS-Net.

§.§ Experiments with CAVE Dataset

Besides ICVL dataset,  we ran all the methods on the CAVE dataset. The CAVE dataset contains 32 HSIs with size of $512\times 512\times 31$. we directly used the trained models from ICVL dataset to test deep learning based methods.  As can be seen from Table <ref>, our method produces the best denoising results in all three cases and surpasses state-of-the-art  deep learning-based QRNN3D and model-based BM3D as well as BM4D.   Fig. <ref> and Fig. <ref> respective depicts false-color image and spectrum of chart and stuffed toy HSI before and after denoising.  It can be observed that our SMDS-Net produces more appealing denoising quality than competing methods which either produce oversmooth image/spectrum, see Fig. <ref> (c) and Fig. <ref> (k), fail to remove all the noises, see Fig. <ref> (d)-(i) or introduce annoying artefacts, see Fig. <ref> (o). This experiment further shows the generalization and strong denoising capability of our SMDS-Net.

Denoising results on real-world Urban dataset. The false-color images were generated by combining bands 109, 141, 209.  The proposed SMDS-Net provides the best visual results. (Best view on screen with zoom)

Denoising results of pixel (249, 216) in real-world  Urban data set.

§.§ Experiments with Real-world Remote Sensing HSI

Further, we ran all the competing methods on the widely used Urban dataset [https://www.agc.army.mil/] to evaluate their real-world denoising ability. Urban HSI is a remote sensing HSI containing $307\times 307$ pixels and 210 bands. It is contaminated with very complex noises, especially in bands 1, 109, 141, 208, 209, making it very suitable to test the denoising ability of methods. Some typical bands are shown in Fig. <ref>. Though there are large differences between close-range HSI and remote sensing HSIs in terms of spatial and spectral resolution, we directly use the model trained on ICVL HSIs rather than fine-tuning the model on remote sensing HSIs.  Fig. <ref> shows the band-wise denoising results with respect to bands 1, 109 and 208 and Fig. <ref> presents the false-color images  by compositing bands 109, 141, 209.   As can be seen, alternative methods either can not remove all the noises, for example, LLRT, MTSNMF, LRTDTV, HSI-SDeCNN, and HSID-CNN, or lose some important textures such as QRNN3D.  In contrast, our method achieves the best denoising results by removing most of the noises while retaining important structures. Fig. <ref> presents the spectral signature of pixel (249, 216) before and after denoising.  Our SMDS-Net also gives very competitive performance.  This experiment further verifies the strong denoising capability and generalization ability of the proposed SMDS-Net.

§.§ Network Analysis

Here, we conduct a throughout study on our SMDS-Net to show its virtue of strong learning ability,  attractive interpretability,  less requirement of training samples and also less number of parameters. All the experiments were conducted on ICVL HSI dataset.

The impact of the number of training samples.

§.§.§ Study of the Number of Training Samples
To further examine the learning capacity of our proposed network, we change the number of training samples from 10 to 100.  As shown in Fig. <ref>, more training samples contribute to network learning, leading to better denoising performance. However, the denoising performance changes are not very significant when the number of training samples increases from 20 to 100, indicating the strong learning ability of SMDS-Net. This makes our method very practical for real-world HSI denoising where training samples are very rare and and expensive to acquire.

The impact of dictionary learning(with/without).

[Initial dictionary]

Visualization of  learned dictionaries for SMDS-Net  when $\{\D_j, \C_j, \W_j\}_{j=1,2,3}\in \Rbb^{9\times 9}$.

§.§.§ Study of Dictionary Learning
Here, we conduct an ablation study to show the importance of dictionary learning of SMDS-Net. As for the case without dictionary learning, we direct set $\D_j$, $\W_j$, and $\C_j$ as DCT basis rather than parameters to be learned. Fig. <ref> presents the effect of dictionary learning. It can be seen that there is a significant improvement in using dictionary learning because it can increase the ability of the network to model and adapt the data. Fig. <ref> visualizes the learned dictionaries which are constructed by the Kronecker product of $\D$, $\W$ and $\C$ along arbitrary dimensions. It can be seen that, on the one hand, there are  significant differences between the initialized DCT basis and the learned DCT basis, which means that the network learns to fit the data. On the other hand, there are differences between dictionaries, see (b)-(d), (e)-(g), (h)-(j), which proves that the learned dictionaries can capture the structural features of all dimensions.

§.§.§ Study of Key Parameters

In general, the number of unfolding ($K$), the dictionary number ($[M_1,M_2,M_3]$) and the cube size $([I_1, I_2, I_3])$ are three factors influencing the denoising effectiveness  of SMDS-Net.  Here, we conduct a comprehensive study on these key parameters.  We add  additive white Gaussian  noise (AWGN) to every band whose  stand deviation $\sigma$ is set between 0 and 95. The averaged denoising performance  is reported.

Impact of Dictionary Number  We measure the influence of the dictionary number   in Table <ref>.  For simplicity, we set $M_1=M_2=M_3$ and choose the value from the set of $\{5,7,9,11,13\}$. As shown in  Table <ref>, using larger dictionaries allows for better denoising performance at the cost of more computation and parameters to be learned. Moreover, increasing the dictionary number from  $[5,5,5]$ to  $[13,13,13]$, we only  measure a difference of 0.15dB. Therefore, the dictionary number is suggested to set as $[9,9,9]$ considering a  trade-off between speed and the number of parameters.

The influence of the  dictionary number.

$[M_1,M_2,M_3]$   [5, 5, 5]   [7, 7, 7]   [9, 9, 9]   [11, 11, 11]   [13, 13, 13]

PSNR   39.21    39.31   39.34   39.43   39.36

SSIM   0.9583   0.9576   0.9593   0.9595   0.9586

SAM   0.0307   0.0306   0.0303    0.0305   0.0304

#Parameters   1155   2625   5103   8877   14235

The influence of unfolding iteration $K$.

$K$   3   6   9   12   15

PSNR    39.21   39.34    39.33   39.41    39.15

SSIM   0.9589   0.9593   0.9593   0.9587   0.9545

SAM   0.0305   0.0303   0.0306   0.0303   0.0304

#Parameters   2916   5103   7290   9477   11664

Impact of Unfolded Iteration $K$:  $K$ indexes the number of unfolding and also directly controls the depth of the network.   Table <ref> presents the denoising result with respect to different $K$.   SMDS-Net takes
the advantages of backpropagation and forward propagation to establish the relationship between current iteration and historical iterations which enforces the network to simultaneously takes historical information and current information into consideration for parameter update. This makes $K$ have insignificant influence on denoising performance as shown in Table <ref>, implying layer number of SMDS-Net is easy to choose. When setting $K=6$, the number of parameters is 5103 which is noticeably less than that of QRNN3D <cit.> and HSI-SDeCNN <cit.>, i.e., over 100 thousand. Small number of parameters also reduce the requirement of training samples size.  This attributes to the adoption of spectral low-rank projection and spatial multidimensional sparse modeling.

The influence of   cube size.

$[I_1,I_2,I_3]$   [3, 3, 3]   [5, 5, 5]   [7, 7, 7]   [9, 9, 9]   [11, 11, 11]

PSNR      37.05   38.78     39.13    39.34       39.34

SSIM   0.9220   0.9537   0.9575   0.9593    0.9595

SAM   0.0358   0.0316   0.0311   0.0303   0.0305

#Parameters   4617   4779   4941   5103   5265

Study of Cube Size:  Here, we analyze the influence of cube size $[I_1, I_2, I_3]$ on HSI denoising.  As the experiment in the study of dictionary number,  $I_1$, $I_2$ and $I_3$ are set the same.  By changing  the size from 3 to 11 with an interval of 2, Table <ref> presents the denoising ability of SMDS-Net. Larger cube size is able to include more spatial nearby pixels for representation, therefore leading to better denoising performance especially in highly noisy cases.   From this experiment, it is recommended that the cube size is set as $[9,9,9]$, balancing the denoising performance and the number of parameters.

§ CONCLUSION

In this paper, we introduce a model guided deep neural network, i.e., SMDS-Net for HSI denoising.
SMDS-Net simultaneously considers the spectral-spatial correlation, spectral low-rankness, and spatial sparsity priors in a unified end-to-end network.  Experimental results show  SMDS-Net provides the state-of-the-art denoising performance with strong learning capacity and higher interpretability, highlighting the benefit of equipping DL with model priors.   In the future, we will integrate the nonlocal similarity prior of HSIs into the network for more powerful denoising performance.


11institutetext: Department of Mathematics and Statistics, University of
Turku, FI-20014 Turku, Finland 11email<EMAIL_ADDRESS>22institutetext:  LIP,
ENS de Lyon – CNRS – INRIA – UCBL – Université de Lyon , École Normale
Supérieure de Lyon, 69007 Lyon, France 22email<EMAIL_ADDRESS>33institutetext: Faculty of Computing, University of Latvia, Raina bulv. 19,
Rīga, LV-1586, Latvia 33email<EMAIL_ADDRESS>44institutetext: Turku Centre for
Computer Science (TUCS)

Mika Hirvensalo, Etienne Moutot, and Abuzer Yakaryılmaz

# On the computational power of affine automata

Mika Hirvensalo,, Partially supported by Väisälä Foundation1144 Etienne
Moutot, Partially supported by TUCS COM3-project and ANR project CoCoGro
(ANR-16-CE40-0005)22 Abuzer Yakaryılmaz, Partially supported by TUCS
COM3-project and ERC Advanced Grant MQC33

###### Abstract

We investigate the computational power of affine automata (AfAs) introduced in
[4]. In particular, we present a simpler proof for how to change the cutpoint
for any affine language and a method how to reduce error in bounded error
case. Moreover, we address to the question of [4] by showing that any affine
language can be recognized by an AfA with certain limitation on the entries of
affine states and transition matrices. Lastly, we present the first languages
shown to be not recognized by AfAs with bounded-error.

###### Keywords:

non-classical models of automata, affine automata, cutpoint languages, bounded
error, compact sets, error reduction

##  1 Introduction

Finite automata are interesting computational models because of their
simplicity, compared to more complex models like pushdown automata or Turing
machines. They also represent a very concrete restriction on computation: they
only have a finite memory. A lot of different automata models have been
studied during the years, such as deterministic [10], probabilistic [8] and
quantum [2] ones. All these models share two common features: the state vector
set is compact and the acceptance function can be interpreted as linear. The
linearity is desirable because of mathematical simplicity, but on the other
hand, it may represent a limitation on the computational power.

Jeandel [5] demonstrated that in the bounded-error acceptance model the finite
automata with compact state set accept only regular languages. Hence the
compactness property of the state set may be one very important limiting the
computational power, but since most known models have a compact state set, it
remains open how much the compactness of the state set actually contributes.

Recently, A. Díaz-Caro and A. Yakaryılmaz introduced a new model, called
affine automata [4], also investigated in [13] and [3]. It is a purely
theoretical model, which means that it cannot be implemented by a physical
device like quantum automata. But it allows us to investigate on the power of
interference caused by negatives amplitudes in the computation, like in the
quantum case. Moreover, this model allows us to study the effect of state set
compactness, since unlike quantum automata, affine ones have an unbounded
state set. In addition, the final operation corresponding to quantum
measurement cannot be interpreted as linear, but it is analogous to
renormalization in Kondacs-Watrous [6] and Latvian [1] quantum automata
models.

In this paper, we present some stability results (Section 3): we show how to
obtain a new AfA from two AfAs by tensoring and direct sum. Then, we present a
simpler proof for how to change the cutpoint for any affine language and an
error reduction method in bounded error case.

Any entry of an affine state or a transition matrix can be arbitrarily away
from zero. Here, by addressing to the question of [4], we show that (Section
4) any affine language can be recognized by an AfA with the restriction that
all the entries of transition matrices are in the interval $[-1,1]$. We also
show that by an additional state we can guarantee that any AfA can start its
computation from the first deterministic state.

Finally, we present (Section 5) the first languages shown not to be recognized
by any bounded-error AfA.

##  2 Preliminaries

We denote the input alphabet $\Sigma$ and the empty string $\varepsilon$.

Probabilistic automata are a generalization of deterministic finite automata
that can make random choices [9]. Formally, a probabilistic finite automaton
(PFA) $P$ is a 5-tuple $P=(E,\Sigma,\\{M_{x}\leavevmode\nobreak\
|\leavevmode\nobreak\ x\in\Sigma\\},e_{s},E_{a}),$ where
$E=\\{e_{1},\dots,e_{k}\\}$ is the finite set of states of $P$,
$\\{M_{x}\leavevmode\nobreak\ |\leavevmode\nobreak\ x\in\Sigma\\}$ is the set
of stochastic transition matrices (all their coefficients are real numbers in
$[0,1]$ and their columns sums up to $1$), $v_{0}$ is the initial
probabilistic state (the probability distribution on the states), and
$E_{a}\subseteq E$ is the set of accepting states. The computation starts in
$v_{0}$, and then the given input, say $w=w_{1}\cdots w_{n}\in\Sigma^{*}$ for
some $n>0$, is read once from left to right symbol by symbol. For each symbol
the corresponding transition matrix is applied:
$\vec{v}_{f}=M_{w}\vec{v}_{0}=M_{w_{n}}\cdots M_{w_{1}}\vec{v}_{0}.$ Remark
that if $w=\varepsilon$, $v_{f}=v_{0}$. The accepting probability of $P$ on
$w$ is given by

$f_{P}(w)=\vec{p}M_{w}\vec{v}_{0},$ (1)

where $\vec{p}=\begin{pmatrix}\delta_{1}&\cdots&\delta_{k}\end{pmatrix}$ and
$\delta_{i}=1$ if $e_{i}\in E_{a}$ and $0$ if $e_{i}\notin E_{a}$.

Affine automata are a generalization of PFAs allowing negative transition
values. Only allowing negative values in the transition matrices does not add
any power (generalized probabilistic automata are equivalent to usual ones
[11]), but affine automata introduces also a non-linear behaviour. The
automaton acts like usual generalized probabilistic automaton until the last
operation, a non-linear operation called weighting.

A vector $\vec{v}\in\mathbb{R}^{n}$ is an affine vector if and only if its
coordinates sums up to $1$. A matrix $M$ is an affine matrix if and only if
all its columns are affine vectors. Remark that if $M$ and $N$ are affine
matrices , then $MN$ is also an affine matrix. In particular, if $\vec{v}$ is
an affine vector, then $M\vec{v}$ is also an affine vector.

Formally, an affine finite automaton (AfA) $A$ is a 5-tuple

$A=(E,\Sigma,\\{M_{x}\leavevmode\nobreak\ |\leavevmode\nobreak\
x\in\Sigma\\},\vec{v}_{0},E_{a})$

where all components exactly the same as for probabilistic automata by
replacing stochastic property with affine one in the initial state and
transition matrices.

As in PFAs, after reading a word $w=w_{1}\cdots w_{n}$, the final state of $A$
is $\vec{v}_{f}=M_{w}\vec{v}_{0}$ like in the probabilistic case, but the
function $f_{A}:\Sigma^{*}\to[0,1]$ computed by $A$ is defined as

$f_{A}(w)=\frac{\sum_{e_{i}\in E_{a}}|(v_{f})_{i}|}{\sum_{e_{i}\in
E}|(v_{f})_{i}|},$ (2)

and referred as the accepting value of $A$ on $w$. Similar to projective
measurements, we can rewrite Eq. (2) as given below. First, we define a
projection matrix based on $E_{a}$: $P_{A}=P=\begin{pmatrix}\delta_{1}&&&\\\
&\delta_{2}&&\\\ &&\ddots&\\\ &&&\delta_{n}\\\ \end{pmatrix},$ $\text{ where
}\delta_{i}=\left\\{\begin{matrix}1\text{ if }e_{i}\in E_{a}\\\ 0\text{
otherwise}\end{matrix}\right..$
Then, we can denote $f_{A}(\cdot)$ as

$f_{A}(w)=\frac{|PM_{w}\vec{v}_{0}|}{|M_{w}\vec{v}_{0}|}.$ (3)

Notice that the final value for PFA $P$ (1) is defined as matrix product
$\vec{v}_{f}\mapsto\vec{p}.\vec{v}_{f}$, which is a linear operation on
$\vec{v}_{f}$. On the other hand, computing final value from $\vec{v}_{f}$ as
in (3) involves nonlinear operations
$\displaystyle\vec{v}_{f}\mapsto\frac{|P\vec{v}_{f}|}{|\vec{v}_{f}|}$ due to
absolute value and normalization of affine states having length greater than
1.

Given a function $f:\Sigma^{*}\to[0,1]$ computed by an automaton (stochastic
or affine), there are different ways of defining the language of an automaton.
The natural one is as follows: A language $L\subseteq\Sigma^{*}$ is recognized
by an automaton $A$ with cutpoint $\lambda$ if and only if

$L=\\{w\in\Sigma^{*}\leavevmode\nobreak\ |\leavevmode\nobreak\
f_{A}(w)>\lambda\\}.$

These languages are called cutpoint languages. In the case of probabilistic
(resp. affine automata), the set of cut-point languages are called stochastic
languages (resp. affine languages) and denoted by $\mathsf{SL}$ (resp.
$\mathsf{AfL}$).

A stronger condition is to impose that accepted and rejected words are
separated by a gap: the cutpoint is said to be isolated: A language $L$ is
recognized by an automaton $A$ with isolated cutpoint $\lambda$ if and only if
there exist $\delta>0$ such that $\forall w\in L,f_{A}(w)\geq\lambda+\delta$,
and $\forall w\notin L,f_{A}(w)\leq\lambda-\delta$.

As we shall see, for affine automata it is always possible to shift the
cutpoint $\lambda\in(0,1)$ to $\lambda=\frac{1}{2}$, and hence this notion of
isolated cutpoint becomes equivalent to the bounded error recognition:
Language $L\subseteq\Sigma^{*}$ is said to be recognized by an automaton $A$
with bounded error if and only if there exists $\varepsilon>0$ such that
$\forall w\in L,f_{A}(w)\geq 1-\varepsilon$, and $\forall w\notin
L,f_{A}(w)\leq\varepsilon$.

The set of languages recognized with bounded error (or isolated cutpoint)
affine automata is denoted by $\mathsf{BAfL}$.

A classical result by Rabin [9] shows that isolated cutpoint stochastic
languages are regular (denoted $\mathsf{REG}$). Rabin’s proof essentially
relies on two facts: 1) the function mapping the final vector into $[0,1]$ is
a contraction, and 2) the state vector set is bounded.

By modifying Rabin’s proof, it is possible to show that also many quantum
variants of stochastic automata obey the same principle [7] bounded-error
property implies the regularity of the accepted languages. In fact, E. Jeandel
generalized Rabin’s proof by demonstrating that the compactness of the state
vector set together with the continuity of the final function are sufficient
to guarantee the regularity of the accepted language if the cutpoint is
isolated [5].

In the affine case however, the vector states do not lie in a compact set, we
cannot prove that $\mathsf{BAfL}=\mathsf{REG}$ like in the probabilistic (or
even quantum) case [5]. In fact, it is even the contrary:
$\mathsf{REG}\subsetneq\mathsf{BAfL}$ [4].

We close this section by three basic facts. The following three operations on
the state sets will be useful, when constructing new automata from the
existing ones:

  * •

$\overline{E}=\\{e_{i}\leavevmode\nobreak\ |\leavevmode\nobreak\ e_{i}\notin
E\\}$ the complement of $E$,

  * •

$E_{a}\times E_{b}=\\{(e_{i},e_{j})\leavevmode\nobreak\ |\leavevmode\nobreak\
e_{i}\in E_{a},e_{j}\in E_{b}\\}$ the Cartesian product of $E_{a}$ and
$E_{b}$,

  * •

$E_{a}\cup E_{b}=\\{e_{i}\leavevmode\nobreak\ |\leavevmode\nobreak\ e_{i}\in
E_{a}\text{ or }e_{i}\in E_{b}\\}$ the union of $E_{a}$ and $E_{b}$.

The following lemma shows how to formulate the above operations by using the
formalism of projection matrices.

###### Lemma 1

Let $E$ be the set of all states, $E_{a},E_{b}\subseteq E$ and $P_{a}$,
$P_{b}$ be the projections associated to them. Then

  * •

$P$ is the projection associated to the complement $\overline{E_{a}}$ if and
only if $P=I-P_{a}$, and,

  * •

$P$ is the projection associated to $E_{a}\times E_{b}$ if and only if
$P=P_{a}\otimes P_{b}$.

###### Lemma 2

Let $E$ be the set of all states, $E_{a},E_{b}\subseteq E$, such that
$E_{a}\cap E_{b}=\emptyset$. Let $P_{a}$ and $P_{b}$ be the projections
associated to them. Then,

  * •

$P$ is the projection associated to $E_{a}\cup E_{b}$ if and only if
$P=P_{a}+P_{b}$, and,

  * •

for any matrix $M$ and vector $\vec{v}$,
$|PM\vec{v}|=|P_{a}M\vec{v}|+|P_{b}M\vec{v}|$.

###### Lemma 3

If $A$ and $B$ are affine matrices, then $A\otimes B$ is also affine.
Moreover, $|A\otimes B|=|A||B|$.

##  3 Stability Results

The main results of this section are stability results. The first are about
the functions of affine automata. They provide a way to prove an error
reduction theorem. We then use this theorem to show the stability of bounded-
error affine languages under intersection and union.

###### Proposition 1

Let $f$, $g$ be functions computed by affine automata, then there exists an
affine automaton $\mathcal{C}$ such that $f_{\mathcal{C}}=f\times g$.

###### Proof

The proof is the same as the stochastic case and essentially relies on the
property of tensor product of Lemma 3. ∎

It is easy to design a 2-state PFA $P$ such that
$f_{P}:\Sigma^{*}\rightarrow\alpha$ for $\alpha\in[0,1]$. Thus:

###### Corollary 3.1

Let $f$ be a function computed by an AfA and $\alpha\in[0,1]$, then there
exists an AfA $\mathcal{C}$ such that $f_{\mathcal{C}}=\alpha f$.

###### Proposition 2

Let $f$, $g$ be functions computed by some AfAs and $\alpha,\beta\geq 0$ such
that $\alpha+\beta=1$, then there exists an AfA $\mathcal{C}$ such that
$f_{\mathcal{C}}=\alpha f+\beta g$.

###### Proof

Let $\mathcal{A}=(E^{A},\Sigma,\\{A_{x}\\},\vec{v}_{0}^{A},E_{a}^{A})$ and
$\mathcal{B}=(E^{B},\Sigma,\\{B_{x}\\},\vec{v}_{0}^{B},E_{a}^{B})$ two
automata such that $f=f_{\mathcal{A}}$ and $g=f_{\mathcal{B}}$. The idea here
is to make two copies of $\mathcal{A}\otimes\mathcal{B}$ working in parallel,
one having the final states of $\mathcal{A}$, the other the final states of
$\mathcal{B}$. We define
$\mathcal{C}=(E^{C},\Sigma,\\{C_{x}\\},\vec{v}_{0}^{C},E_{a}^{C})$ by:

$C_{x}=\left(\begin{array}[]{c|c}\\\ \mathcal{A}_{x}\otimes B_{x}&0\\\ \\\
\hline\cr\\\ 0&A_{x}\otimes B_{x}\\\ \\\
\end{array}\right),\vec{v}_{0}^{C}=\left(\begin{array}[]{c}\\\
\alpha(\vec{v}_{0}^{A}\otimes\vec{v}_{0}^{B})\\\ \\\ \hline\cr\\\
\beta(\vec{v}_{0}^{A}\otimes\vec{v}_{0}^{B})\\\ \\\
\end{array}\right),P^{C}=\left(\begin{array}[]{c|c}\\\ P^{A}\otimes I_{n}&0\\\
\\\ \hline\cr\\\ 0&I_{k}\otimes P^{B}\\\ \\\ \end{array}\right),$

with $P^{A}$, $P^{B}$ and $P^{C}$ be the projections on $E^{A}_{a}$,
$E^{B}_{a}$ and $E^{C}_{a}$. Thus,

$\displaystyle f_{\mathcal{C}}(w)$ $\displaystyle=\frac{\alpha|(P^{A}\otimes
I_{n})(A_{x}\otimes
B_{x})(\vec{v}_{0}^{A}\otimes\vec{v}_{0}^{B})|+\beta|(I_{k}\otimes
P^{B})(A_{x}\otimes
B_{x})(\vec{v}_{0}^{A}\otimes\vec{v}_{0}^{B})|}{(\alpha+\beta)|(A_{x}\otimes
B_{x})(\vec{v}_{0}^{A}\otimes\vec{v}_{0}^{B})|}$
$\displaystyle=\alpha\frac{|P^{A}A_{w}\vec{v}_{0}^{A}|}{|A_{w}\vec{v}_{0}^{A}|}+\beta\frac{|P^{B}B_{w}\vec{v}_{0}^{B}|}{|B_{w}\vec{v}_{0}^{B}|}=\alpha
f(w)+\beta g(w).$

∎

The first consequence of these stability results is a really short proof for
shifting the cutpoint of an affine automaton. Although the construction in [4]
gives a much more compact automata in term of number of states, our
construction is simpler, and does not require as many specific cases.

###### Proposition 3

Let $\mathcal{A}$ be an affine automaton and
$\lambda_{1},\lambda_{2}\in[0,1]$. There exists an affine automaton
$\mathcal{B}$ such that

  * •

$f_{\mathcal{A}}(w)>\lambda_{1}\Leftrightarrow f_{\mathcal{B}}(w)>\lambda_{2}$
and

  * •

$f_{\mathcal{A}}(w)=\lambda_{1}\Leftrightarrow
f_{\mathcal{B}}(w)=\lambda_{2}$.

###### Proof

First we suppose $\lambda_{1}\neq 1$. Let $\mathcal{B}$ the automaton such
that $f_{\mathcal{B}}=\alpha f_{\mathcal{A}}+(1-\alpha)1$, with
$\alpha=\frac{1-\lambda_{2}}{1-\lambda_{1}}$. Then
$f_{\mathcal{A}}>\lambda_{1}\Rightarrow
f_{\mathcal{B}}>\frac{(1-\lambda_{2})\lambda_{1}+\lambda_{2}-\lambda_{1}}{1-\lambda_{1}}=\lambda_{2}$.
And one has the same with $=$ or $<$.

For $\lambda_{1}=1$ it is even simpler, one has just to “resize” the function
by taking $\mathcal{B}$ such that
$f_{\mathcal{B}}=\lambda_{2}f_{\mathcal{A}}$. And then,
$f_{\mathcal{A}}=1\Rightarrow f_{\mathcal{B}}=\lambda_{2}$, and same for $<$.
∎

Using the same kind of construction we can prove that bounded-error mode, it
is always possible to reduce the error. Reducing the error means increasing
the gap between accepted and rejected words. The error probability could even
be made as close to zero as one wants.

###### Lemma 4

Let $f$ be a function computed by affine automaton, then there exists an
affine automaton $\mathcal{B}$ such that $f_{\mathcal{B}}=f^{2}(3-2f)$.

###### Proof

Let $\mathcal{A}=(E,\Sigma,\\{A_{x}\\},\vec{v}_{0},E_{a})$ such that
$f=f_{\mathcal{A}}$. The automaton $\mathcal{B}$ will run 3 copies of
$\mathcal{A}$ in parallel, and its final states are made to accept if 2 or 3
copies of $\mathcal{A}$ accept and reject otherwise (i.e. taking the majority
answer). Formally, $\mathcal{B}=(E\otimes E\otimes
E,\Sigma,\\{B_{x}\\},\vec{v}^{\prime}_{0},E_{a}^{\prime})$ with

$B_{x}=A_{x}\otimes A_{x}\otimes A_{x},$
$\vec{v}^{\prime}_{0}=\vec{v}_{0}\otimes\vec{v}_{0}\otimes\vec{v}_{0},$
$E^{\prime}_{a}=\left(E_{a}\times E_{a}\times
E_{a}\right)\cup\left(\overline{E_{a}}\times E_{a}\times
E_{a}\right)\cup\left(E_{a}\times\overline{E_{a}}\times
E_{a}\right)\cup\left(E_{a}\times E_{a}\times\overline{E_{a}}\right).$

Note that the four sets in parenthesis are all pairwise disjoints. Let $P$ and
$P^{\prime}$ be the projections associated to $E_{a}$ and $E^{\prime}_{a}$.
Then,

$P^{\prime}=P\otimes P\otimes P+(I-P)\otimes P\otimes P+P\otimes(I-P)\otimes
P+P\otimes P\otimes(I-P).$

And by Lemma 1,

$\displaystyle
f_{\mathcal{B}}(w)=\frac{|P^{\prime}B_{w}\vec{v}^{\prime}_{0}|}{|B_{w}\vec{v}^{\prime}_{0}|}$
$\displaystyle=\frac{|PA_{w}\vec{v}_{0}|^{3}+3|PA_{w}\vec{v}_{0}|\left(|A_{w}\vec{v}_{0}|-|PA_{w}\vec{v}_{0}|\right)}{|A_{w}\vec{v}_{0}|^{3}}$
$\displaystyle=f(w)^{3}+3f(w)^{2}(1-f(w))$ $\displaystyle=f(w)^{2}(3-2f(w)).$

∎

######  Proposition 4 (error reduction)

Let $L\in\mathsf{BAfL}$. There exists an affine automaton $\mathcal{A}$ such
that:

  * •

$\forall w\in L,f_{\mathcal{A}}(w)\geq\frac{3}{4}$

  * •

$\forall w\notin L,f_{\mathcal{A}}(w)\leq\frac{1}{4}$

###### Proof

We do not detail the proof, but the idea is simple: mapping $x\to x^{2}(3-2x)$
has attracting points at $x=0$ and $x=1$. Iterating the mapping, any point
$\neq\frac{1}{2}$ will tend to $0$ or $1$. ∎

This technique could be applied to get any constant instead of $\frac{1}{4}$,
to have an error bound as small as one wants.

This error reduction theorem also applies to probabilistic automata, but is
not very interesting because in the probabilistic case it is known that
bounded-error languages are exactly regular languages [5], and hence the error
probability could always be 0. In our case, bounded-error languages are more
complex than regular languages. But thanks to this error reduction, they are
stable under union, intersection, and complement, just like regular languages.

###### Proposition 5

Let $L_{A},L_{B}\in\mathsf{BAfL}$. Then

  * •

$L_{A}\cup L_{B}\in\mathsf{BAfL}$,

  * •

$L_{A}\cap L_{B}\in\mathsf{BAfL}$,

  * •

$\overline{L_{A}}\in\mathsf{BAfL}$.

###### Proof

Let $\mathcal{A}$ and $\mathcal{B}$ be automata recognizing $L_{A}$ and
$L_{B}$ with error bound $\varepsilon$ at most $\frac{1}{4}$ (thanks to
Theorem 4). We define $\mathcal{C}$ and $\mathcal{D}$ such that
$f_{\mathcal{C}}=\frac{1}{2}(f_{\mathcal{A}}+f_{\mathcal{B}})$ and
$f_{\mathcal{D}}=f_{\mathcal{A}}f_{\mathcal{B}}$. Let $w\in\Sigma^{*}$. We
study the 4 possible options depending on the membership of $w$ to $L_{A}$ and
$L_{B}$.

  * •

$w\in L_{A},w\in L_{B}$ (i.e. $w\in L_{A}\cup L_{B},w\in L_{A}\cap L_{B}$)
$\Rightarrow f_{C}\geq\frac{3}{4}$ and $f_{\mathcal{D}}\geq\frac{9}{16}$,

  * •

$w\in L_{A},w\notin L_{B}$ (i.e. $w\in L_{A}\cup L_{B},w\notin L_{A}\cap
L_{B}$) $\Rightarrow f_{C}\geq\frac{3}{8}$ and
$f_{\mathcal{D}}\leq\frac{1}{4}$,

  * •

$w\notin L_{A},w\in L_{B}$ (i.e. $w\in L_{A}\cup L_{B},w\notin L_{A}\cap
L_{B}$) $\Rightarrow f_{C}\geq\frac{3}{8}$ and
$f_{\mathcal{D}}\leq\frac{1}{4}$,

  * •

$w\notin L_{A},w\notin L_{B}$ (i.e. $w\notin L_{A}\cup L_{B},w\notin L_{A}\cap
L_{B}$) $\Rightarrow f_{C}\leq\frac{1}{4}$ and
$f_{\mathcal{D}}\leq\frac{1}{16}$.

Because $\frac{3}{8}>\frac{1}{4}$ and $\frac{9}{16}>\frac{1}{4}$,
$\mathcal{C}$ and $\mathcal{D}$ are deciding $L_{A}\cup L_{B}$ and $L_{A}\cap
L_{B}$ with bounded error.

For the complement one has just to make a copy of $\mathcal{A}$ with accepting
states $\overline{E_{a}}$. The resulting function will be $1-f_{\mathcal{A}}$,
leading to accept the rejected words of $\mathcal{A}$ and vice-versa. ∎

##  4 Equivalent Forms of Affine Automata

General affine automata are hard to study because of the lack of structure of
their transition matrices and state vectors. We provide here some equivalent
forms which have more restrictive properties. These equivalent forms are
useful not only because it provides simpler equivalent models but also because
they provide a way understand the power of affine computation.

The first result is that assuming the initial affine (probabilistic) state as
the first deterministic state does not change the power of AfAs (PFAs).

###### Proposition 6

Let $\mathcal{A}$ be an affine automaton with $n$ states, there exist
$\mathcal{B}$ with $n+1$ states with the initial state $(1,0,\ldots,0)$ and
such that $f_{\mathcal{A}}=f_{\mathcal{B}}$.

###### Proof

Let $\mathcal{A}=(E,\Sigma,\\{A_{x}\\},\vec{v}_{0},E_{a})$. Then,
$\mathcal{B}=(E\cup\\{e^{\prime}\\},\Sigma,\\{B_{x}\\},\vec{v}^{\prime}_{0},E_{a})$,
with $\vec{v}^{\prime}_{0}=(1,0,\ldots,0)^{T}$ and
$B_{x}=\left(\begin{array}[]{c|ccc}0&0&\cdots&0\\\ \hline\cr\\\
A_{x}\vec{v}_{0}&&A_{x}\\\ \\\ \end{array}\right).$ Thus we can deduce
$f_{\mathcal{B}}=f_{\mathcal{A}}$ from
$B_{w}\vec{v}^{\prime}_{0}=B_{w_{n}}\ldots
B_{w_{2}}B_{w_{1}}\vec{v}^{\prime}_{0}=\left(\begin{array}[]{c|ccc}0&0&\cdots&0\\\
\hline\cr\\\ A_{w}\vec{v}_{0}&&A_{w}\\\ \\\
\end{array}\right)\begin{pmatrix}1\\\ 0\\\ \vdots\\\
0\end{pmatrix}=\begin{pmatrix}0\\\ \\\ A_{w}\vec{v}_{0}\\\ \\\ \end{pmatrix}.$
∎

Then we prove that one could also assume that all state vectors and transition
matrices have coefficients only in $[-1,1]$.

###### Proposition 7

Any language in $\mathsf{AfL}$ can be recognized by a AfA $\mathcal{B}$ with
cutpoint $\frac{1}{2}$ such that each entry of affine states during the
computation is always in $[-1,1]$.

###### Proof

Let
$\mathcal{A}=(E=\\{e_{1},\ldots,e_{k}\\},\Sigma,\\{A_{x}\\},v_{0}=(1,0,\ldots,0)^{T},E_{a})$
be an AfA such that $w\in L\Leftrightarrow f_{\mathcal{A}}(w)>\frac{1}{2}$,
and $C=\max_{x,i,j}|(A_{x})_{i,j}|$. Then, $\mathcal{B}$ is as follows:

$\mathcal{B}=(E\cup\\{e_{n+1},e_{n+2}\\},\Sigma,\\{B_{x}\\},\vec{v}^{\prime}_{0},E_{a}\cup\\{e_{n+1}\\})\leavevmode\nobreak\
\leavevmode\nobreak\ \mbox{ with}$
$B_{x}=\frac{1}{2kC}\left(\begin{array}[]{ccc|cc}&&&0&0\\\
&2A_{x}&&\vdots&\vdots\\\ &&&0&0\\\ \hline\cr kC-1&\ldots&kC-1&2kC&0\\\
kC-1&\ldots&kC-1&0&2kC\end{array}\right)\mbox{ and
}\vec{v}^{\prime}_{0}=(1,0,\ldots,0)^{T}.$

Then, with $w=w_{1}\cdots w_{n}$, we can deduce that

$B_{w}=B_{w_{n}}\cdots
B_{w_{2}}B_{w_{1}}=\frac{1}{2(kC)^{n}}\left(\begin{array}[]{ccc|cc}&&&0&0\\\
&2A_{w}&&\vdots&\vdots\\\ &&&0&0\\\
\hline\cr(kC)^{n}-1&\ldots&(kC)^{n}-1&2(kC)^{n}&0\\\
(kC)^{n}-1&\ldots&(kC)^{n}-1&0&2(kC)^{n}\end{array}\right),$

which gives the final values of the states:

$\vec{v}^{\prime}_{f}=B_{w}\vec{v}^{\prime}_{0}=\frac{1}{(kC)^{n}}\begin{pmatrix}\vdots\\\
\vec{v}_{f}\\\ \vdots\vspace{0.6em}\\\ \frac{(kC)^{n}-1}{2}\vspace{0.7em}\\\
\frac{(kC)^{n}-1}{2}\end{pmatrix}.$

Since $|(\vec{v}_{f})_{i}|\leq k^{n-1}C^{n}$, it is clear that
$|(\vec{v}^{\prime}_{f})_{i}|\leq[-1,1]$: the values of the states are
bounded. Now, one has

$f_{\mathcal{B}}=\frac{|PA_{w}\vec{v}_{0}|+\frac{(kC)^{n}-1}{2}}{|A_{w}\vec{v}_{0}|+(kC)^{n}-1},$

and so,

$\displaystyle w\in L\Leftrightarrow f_{\mathcal{A}}>\frac{1}{2}$
$\displaystyle\Leftrightarrow|PA_{w}\vec{v}_{0}|>\frac{1}{2}|A_{w}\vec{v}_{0}|$
$\displaystyle\Leftrightarrow|PA_{w}\vec{v}_{0}|+\frac{(kC)^{n}-1}{2}>\frac{1}{2}\left(|A_{w}\vec{v}_{0}|+(kC)^{n}-1\right)$
$\displaystyle\Leftrightarrow f_{\mathcal{B}}>\frac{1}{2}.$

∎

##  5 The first languages shown to be not in $\mathsf{BAfL}$

This part is dedicated to prove that some languages are not recognizable by
affine automata. This is an adaptation of the proof of Turakainen [12] for
non-stochastic languages. All the difficulty of exhibiting a non-affine
language relies in the fact that a large majority of non-stochasticity proof
are based on the linearity of the automaton, which is not the case in the
affine case. This proof however, is more based on some “regularity” induced by
the matrix-based operations, and number theoretic properties of languages like
$Prime$. Hence it was possible to adapt it for the affine case, where the only
non-linear operation is the final projection.

Let $L\subseteq a^{*}$ be a unary language. We call lower density of $L$ the
limit

$\underline{dens}(L)=\liminf_{n\rightarrow\infty}\frac{\left|\\{a^{k}\in
L\leavevmode\nobreak\ |\leavevmode\nobreak\ k\leq n\\}\right|}{n+1}.$

Let $(\vec{x}_{n})$ be a sequence of vectors in $\mathbb{R}^{k}$ and
$I=[a_{1},b_{1})\times\cdots\times[a_{k},b_{k})$ be an “interval”. We define
$C(I,n)$ as $C(I,n)=\left|\\{\vec{x}_{i}\mod 1\in I\leavevmode\nobreak\
|\leavevmode\nobreak\ 1\leq i\leq n\\}\right|$.

We say that $(\vec{x}_{n})$ is uniformly distributed mod 1 if and only if for
any $I$ of such type,

$\lim_{n\rightarrow\infty}\frac{C(I,n)}{n}=(b_{1}-a_{1})\cdots(b_{k}-a_{k}).$

###### Proposition 8

If $L\subseteq a^{*}$ satisfies the following conditions:

  1. 1.

dens(L) = 0.

  2. 2.

For all $Q\in\mathbb{N}^{*}$, there exist $h\in\mathbb{N}$ and an infinite
sequence $(n_{i})\in\mathbb{N}^{\mathbb{N}}$ such that $a^{h+n_{i}Q}\subseteq
L$ and for any irrational number $\alpha$, the sequence
$\left((h+n_{i}Q)\alpha\right)_{i\in\mathbb{N}}$ is uniformly distributed mod
1.

Then $L$ is non-affine ($L\notin\mathsf{BAfL}$).

###### Proof

Let’s assume for contradiction that $L\in\mathsf{BAfL}$. Then there exists an
affine automaton $A$ with $s$ states such that

$f_{A}(a^{n})=\frac{|PM^{n}\vec{v}|}{|M^{n}\vec{v}|}$

and there exists $\varepsilon>0$ such that

  * •

$\forall w\in L,f_{A}(w)\geq 1-\varepsilon$,

  * •

$\forall w\notin L,f_{A}(w)\leq\varepsilon$.

Note that

$|M^{n}\vec{v}|=\sum_{i=1}^{s}\left|(M^{n}\vec{v})_{i}\right|\geq\left|\sum_{i=1}^{s}(M^{n}\vec{v})_{i}\right|=1\text{
(triangle inequality).}$

Hence the denominator of $f_{A}$ is never $0$, and so $f_{A}$ is continuous.

Using the Jordan decomposition $M=PJP^{-1}$, one has $M^{n}=PJ^{n}P^{-1}$. So
the coordinates $\vec{v}_{i}$ of $M^{n}\vec{v}$ have the form

$\vec{v}_{i}=\sum_{k=1}^{s}p_{ik}(n)\lambda_{k}^{n}$ (4)

where $\lambda_{i}$ are the eigenvalues of $M$ and $p_{ik}$ are polynomials of
degree less than the degree of the corresponding eigenvalue. Let
$\lambda_{i}=|\lambda_{i}|e^{2i\pi\theta_{i}}$, we assume
$|\lambda_{1}|=\cdots=|\lambda_{s^{\prime}}|>|\lambda_{s^{\prime}+1}|\cdots$.
Let $\lambda=|\lambda_{1}|$ be the largest module of all eigenvalues and $r$
be the maximum degree of all polynomials $p_{ik}$, where $k\leq s^{\prime}$.
Then, one can use (4) to write

$|M^{n}\vec{v}|=\sum_{i\in E}|\vec{v}_{i}|=\lambda^{n}n^{r}\left(\sum_{i\in
E}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi
n\theta_{k}}\right|+g_{E}(n)\right)$

where $a_{ik}$ is the coefficient of degree $r$ of $p_{ik}$ (note that one can
have $a_{ik}=0$ for some $a,k$), and $g_{E}$ a function such that
$\lim_{n\rightarrow\infty}g_{E}(n)=0$. Similarly,

$|PM^{n}\vec{v}|=\sum_{i\in
E_{a}}|\vec{v}_{i}|=\lambda^{n}n^{r}\left(\sum_{i\in
E_{a}}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi
n\theta_{k}}\right|+g_{E_{a}}(n)\right).$

Now let $F(n)=f(a^{n})$. Using the previous equations, one has

$\displaystyle F(n)$ $\displaystyle=\frac{|PM^{n}\vec{v}|}{|M^{n}\vec{v}|}$
$\displaystyle=\frac{\lambda^{n}n^{r}\left(\sum_{i\in
E_{a}}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi
n\theta_{k}}\right|+g_{E_{a}}(n)\right)}{\lambda^{n}n^{r}\left(\sum_{i\in
E}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi
n\theta_{k}}\right|+g_{E}(n)\right)}$ $\displaystyle=\frac{\sum_{i\in
E_{a}}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi
n\theta_{k}}\right|+g_{E_{a}}(n)}{\sum_{i\in
E}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi n\theta_{k}}\right|+g_{E}(n)}.$

We define

$G(n)=\frac{\sum_{i\in E_{a}}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi
n\theta_{k}}\right|}{\sum_{i\in E}\left|\sum_{k=1}^{s^{\prime}}a_{ik}e^{2i\pi
n\theta_{k}}\right|}.$

As $\lim_{n\rightarrow\infty}g_{E_{a}}(n)=0$ and
$\lim_{n\rightarrow\infty}g_{E}(n)=0$, one has $G(n)\sim F(n)$, and so,

$\lim_{n\rightarrow\infty}|F(n)-G(n)|=0.$ (5)

We define $A=\\{k\leavevmode\nobreak\ |\leavevmode\nobreak\ 1\leq k\leq
s^{\prime},\theta_{k}\notin\mathbb{Q}\\}$ the indices of the “first”
eigenvalue angles that are not rational. Let $Q$, $h$ and the sequence
$(n_{i})$ be as in the statement. Using the periodic behaviour induced by
rational angle of eigenvalues, and by taking a subsequence of the initial one,
one can also assume that $(n_{i})$ is such that

$G(h+n_{i}Q)=\frac{\sum_{i\in E_{a}}\left|\sum_{k\in
A}a_{ik}e^{2i\pi(h+n_{i}Q)\theta_{k}}+c\right|}{\sum_{i\in E}\left|\sum_{k\in
A}a_{ik}e^{2i\pi(h+n_{i}Q)\theta_{k}}+d\right|}$

with $c$, $d$ some constants.

By assumption, for all $k\in A$, the sequence
$\left((h+n_{i}Q)\theta_{k}\right)_{i}$ is uniformly distributed modulo 1. The
consequence is that the values $e^{2i\pi(h+n_{i}Q)\theta_{k}}$ are dense in
the unit circle. If for some $n$, $G(h+nQ)<\frac{1}{2}$, there exists
$\varepsilon>0$ such that $G(h+nQ)\leq\frac{1}{2}-\varepsilon$. Then, thanks
to the density argument, there are arbitrarily large values of $i$ for which
$G(h+n_{i}Q)\leq\frac{1}{2}-\frac{\varepsilon}{2}$. Since for $i$ sufficiently
large, $|F(h+n_{i}Q)-G(h+n_{i}Q)|\leq\frac{\varepsilon}{2}$ (using (5)), one
has $F(h+n_{i}Q)\leq\frac{1}{2}$, and so $a^{h+n_{i}Q}\notin L$, contradicting
condition 2 of the statement.

Therefore, $G(h+nQ)\geq\frac{1}{2}$ for large enough $n$. Because $G$ is not
identically equal to $\frac{1}{2}$ (if it is the case, $F$ would be as close
to $\frac{1}{2}$ as one wants, which is impossible since $L\in\mathsf{BAfL}$),
again using density, there must be some $\varepsilon>0$ and $k_{0}$ such that
$G(h+k_{0}Q)\geq\frac{1}{2}+\epsilon$.

First if $A=\emptyset$, it means that all the angles of the eigenvalues
$\theta_{1},\dots,\theta_{s^{\prime}}$ are rational. We can then write them as
$\theta_{k}=\frac{l_{k}}{m_{k}}$. Then $G(n)$ takes a finite number of values,
and these values only depend on $(n\mod m_{1}),\dots,(n\mod m_{s^{\prime}})$.
Let’s call $k_{1}=h+k_{0}Q$ the number where $G$ is larger than $\frac{1}{2}$:
$G(n_{1})>\frac{1}{2}$. $G$ has the same value for all $n\in
Z=\\{k_{1}+km_{1}\cdots m_{s^{\prime}}|k\in\mathbb{N}\\}$ (because for $n$ in
this set, the values of all $(n\mod m_{1}),\dots,(n\mod m_{s^{\prime}})$ are
the same). Then, thanks to (5), one has, for $n\in Z$ sufficiently large,
$F(n)>\frac{1}{2}$, so $\\{a^{n}\leavevmode\nobreak\ |\leavevmode\nobreak\
n\in Z,n\geq n_{1}\\}\subseteq L$. And because
$\left|\\{a^{n}\leavevmode\nobreak\ |\leavevmode\nobreak\ n\in Z,n\geq
n_{1}\\}\right|\sim\frac{n}{m_{1}\cdots m_{s^{\prime}}}$, one has
$\underline{dens}(L)>0$, which contradicts condition 1 of the statement.

Next, if $A\neq\emptyset$. Let

$R((x_{k})_{k\in A})=\frac{\sum_{i\in E_{a}}\left|\sum_{k\in
A}a_{ik}x_{k}+c\right|}{\sum_{i\in E}\left|\sum_{k\in
A}a_{ik}x_{k}+d\right|}.$

Note that $G(h+n_{i}Q)=R((e^{2i\pi(h+n_{i}Q)\theta_{k}})_{k\in A})$. Then,
because the sequences $((h+n_{i}Q)\theta_{k})_{i}$ are uniformly distributed
modulo 1, it follows that any value obtained by the function $R((e^{2i\pi
y_{k}})_{k\in A})$ can be approximated by some $G(h+n_{i}Q)$ with arbitrary
precision. The function $R$ is continuous, therefore there exists an interval
$I=(x_{1},y_{1},...)=((x_{k},y_{k}))_{k\in A}$ on which
$R((x_{k}))>\frac{1}{2}+\frac{\varepsilon}{2}$. So, if $n_{i}$ is large enough
and satisfies

$\left((h+n_{i}Q)\theta_{1}\mod 1,\dots\right)=\left((h+n_{i}Q)\theta_{k}\mod
1\right)_{k\in A}\in I,$

then $G(h+n_{i}Q)>\frac{1}{2}+\frac{\varepsilon}{2}$, which implies
$F(h+n_{i}Q)>\frac{1}{2}$ and hence $a^{h+n_{i}Q}\in L$. Now we just have to
prove that the sequence $(h+n_{i}Q)$ is “dense enough” to have
$\underline{dens}(L)>0$, contradicting again condition 1.
Because of uniform distribution imposed by condition 2, one has

$d=\lim_{i\rightarrow\infty}\frac{C(I,h+n_{i}Q)}{h+n_{i}Q}=\prod_{k\in
A}(y_{k}-x_{k})$

And so for $i$ large enough, $\frac{C(I,h+n_{i}Q)}{h+n_{i}Q}\geq\frac{d}{2}$,
with $a^{h+n_{i}Q}\in L$, implying $\underline{dens}(L)>0$. We have proved
that $L$ cannot be affine. ∎

Turakainen [12] proved that $Prime=\\{a^{p}\mid p\mbox{ is prime}\\}$ and
$Poly(q)=\\{a^{q(n)}\mid n\in\mathbb{N},q(n)\geq 0\\}$ (where $q$ is any
polynomial of degree $>2$ with non-negative coefficients) both satisfy the two
conditions of Theorem 8. Hence they are not in $\mathsf{BAfL}$ .

###### Corollary 5.1

$Prime\notin\mathsf{BAfL}$ and $Poly(q)\notin\mathsf{BAfL}$.

##  6 Conclusion

In this paper we demonstrated that even if they are strictly more powerful,
bounded-error languages of affine automata share stability properties with
regular languages (which are bounded-error languages of stochastic automata).

We also showed that the computational power of affine automata does not come
alone from the unboundedness state vector set: the general model of unbounded
state vector set can always be simulated with a bounded state vector set.
Hence some of the computational power of affine automata comes from the
nonlinear nature of the final projection, at least in the case of unbounded-
error computation.

## References

  * [1] Ambainis, A., Beaudry, M., Golovkins, M., Ķikusts, A., Mercer, M., Thérien, D.: Algebraic results on quantum automata. Theory of Computing Systems 39(1), 165–188 (2006)
  * [2] Ambainis, A., Yakaryılmaz, A.: Automata: From Mathematics to Applications, chap. Automata and Quantum Computing (To appear), (arXiv:1507.01988)
  * [3] Belovs, A., Montoya, J.A., Yakaryılmaz, A.: Can one quantum bit separate any pair of words with zero-error? Tech. Rep. 1602.07967, arXiv (2016)
  * [4] Díaz-Caro, A., Yakaryılmaz, A.: Affine computation and affine automaton. In: Computer Science — Theory and Applications. LNCS, vol. 9691, pp. 1–15. Springer (2016), arXiv:1602.04732
  * [5] Jeandel, E.: Topological automata. Theory of Computing Systems 40(4), 397–407 (2007)
  * [6] Kondacs, A., Watrous, J.: On the power of quantum finite state automata. In: FOCS’97. pp. 66–75 (1997)
  * [7] Li, L., Qiu, D., Zou, X., Li, L., Wu, L., Mateus, P.: Characterizations of one-way general quantum finite automata. Theoretical Computer Science 419, 73–91 (2012)
  * [8] Paz, A.: Introduction to Probabilistic Automata. Academic Press, New York (1971)
  * [9] Rabin, M.O.: Probabilistic automata. Information and Control 6, 230–243 (1963)
  * [10] Sipser, M.: Introduction to the Theory of Computation, 3rd edition. Cengage Learning, United States of America (2013)
  * [11] Turakainen, P.: Generalized automata and stochastic languages. Proceedings of the American Mathematical Society 21, 303–309 (1969)
  * [12] Turakainen, P.: On nonstochastic languages and homomorphic images of stochastic languages. Information Sciences 24(3), 229–253 (1981)
  * [13] Villagra, M., Yakaryılmaz, A.: Language recognition power and succintness of affine automata. In: Unconventional Computation and Natural Computation. LNCS, vol. 9726, pp. 116–129. Springer (2016)