extfieldsofstudy:
Computer Science
Mathematics
Medicine

text:
Theory of Deep Convolutional Neural Networks II: Spherical Analysis

Deep learning based on deep neural networks of various structures and architectures has been powerful in many practical applications, but it lacks enough theoretical verifications. In this paper, we consider a family of deep convolutional neural networks applied to approximate functions on the unit sphere $\mathbb{S}^{d-1}$ of $\mathbb{R}^d$. Our analysis presents rates of uniform approximation when the approximated function lies in the Sobolev space $W^r_\infty (\mathbb{S}^{d-1})$ with $r>0$ or takes an additive ridge form. Our work verifies theoretically the modelling and approximation ability of deep convolutional neural networks followed by downsampling and one fully connected layer or two. The key idea of our spherical analysis is to use the inner product form of the reproducing kernels of the spaces of spherical harmonics and then to apply convolutional factorizations of filters to realize the generated linear features.

Introduction
Deep learning has attracted tremendous attention from various fields of science and technology recently. Wide applications including those in image processing [9] and speech recognition [12] have received great successes. Based on deep neural network structures, it has a strong capability of obtaining data features and distinguishes itself from classical machine learning methods. Though it is successful in practical applications, theoretical assurances are still lacking and need to be investigated. Many attempts have been made trying to understand the practical power of deep neural networks [3,18].
The classical (shallow) neural networks to approximate functions or process data on R d take the form 1) where N is the number of neurons called width, are parameters corresponding to connection weights, biases, and coefficients, w k , x = w k · x is the inner product in R d , and σ : R → R is an activation function. Approximation of functions on subsets of R d by shallow neural networks (1.1) was studied well around the late 1980s. See [19,5] and references therein. As a natural extension of shallow nets, fully connected deep neural networks (DNNs) have been developed since the 1990s. A fully connected DNN of input x = (x 1 , x 2 , . . . , x d ) ∈ R d with J hidden layers of neurons DNNs designed by convolutions are called deep convolutional neural networks (CNNs) and have been very successful in image classification and related applications [16]. Such a network is associated with a sequence of convolutional filters w = {w (j) : Z → R} J j=1 , where w (j) is supported in 0, · · · , S (j) for some S (j) ∈ N called filter length. The convolution of such a filter w supported in {0, · · · , S} with another sequence v = (v 1 , . . . , v D ) is a sequence w * v given by which is supported in {1, · · · , D + S}. Then by restricting the convoluted sequence onto its support, for input x = (x 1 , x 2 , . . . , x d ) ∈ R d , a deep CNNs with J hidden layers of neurons {h (j) : R d → R d j } and widths d j = d j−1 + S (j) is defined iteratively by h (0) (x) = x and By inducing a Toeplitz type convolutional matrix ..,D+S,k=1,...,D associated with a filter w of filter length S and D ∈ N given explicitly by a deep CNN can be regarded as a special sparse form of a fully connected DNNs with the full connection matrices F (j) in (1.2) replaced by the sparse Toeplitz convolutional matrices T (j) := T w (j) with D = d j−1 and S = S (j) for j = 0, 1, . . . , J. That is, (1.3) becomes Compared with fully-connected DNNs, deep CNNs reduce the computational complexity by using at each layer a Toeplitz matrix due to the sparsity and convolutional nature. Throughout the paper we take an identical filter length S (j) ≡ S ∈ N implying {d j = d+jS} and take the rectified linear unit (ReLU) activation function It was shown in [33] that the output layer of any fully-connected DNN can be realized by a downsampled deep CNN with free parameters of the same order, and that deep CNNs can approximate ridge functions of the form g(ξ ·x) with univariate functions g and unknown ξ ∈ R d to the same accuracy with much smaller number of free parameters. Universality of approximation by deep CNNs was also established in [32] where rates of approximation were provided for restrictions of functions from the Sobolev space H r (R d ) with a regularity index r > 2 + d/2. The regularity index r is large when d increases and is essentially needed in the analysis there due to the function regularity on the whole Euclidean space R d . Note that the issue of approximating non-smooth functions by fully connected DNNs has been studied in [28,13]. Moreover, the Sobolev space H r (R d ) requires derivatives of various orders to belong to the L 2 space, while the approximation considered in [32] is measured in the L ∞ norm.
In this paper, we overcome the difficulty in the large regularity index and the inconsistency of L 2 and L ∞ norms for the setting with data from the unit sphere S d−1 in R d . With our novel analysis conducted with spherical harmonic expansions, we can present rates of approximating functions from the Sobolev space W r ∞ (S d−1 ) (to be defined below) on S d−1 , with any positive index r, by downsampled deep CNNs defined in [33] followed by two fully connected layers.
where ⌊u⌋ denotes the integer part of u > 0.
After the last CNN layer, we add two fully connected layers h (J+1) , h (J+2) with widths D 1 , D 2 > 0, respectively, connection matrices F (J+1) , F (J+2) and bias vectors b (J+1) , b (J+2) , to be determined. Precisely, and Such a network with many convolutional layers followed by downsampling operations and very few fully connected layers is quite common in practical applications [16,9]. The hypothesis space of functions on S d−1 induced by our network is given by (1.9)

Main Results on Rates of Approximation
Our target is to establish rates of approximating functions in W r ∞ (S d−1 ) by those from the hypothesis space H J,D 1 ,D 2 ,S defined by (1.9). Since the sums of the rows in the middle of the Toeplitz type matrix (1.4) are equal, we impose for the bias vectors {b (j) } J j=1 of the convolutional layers a restriction b (j) For the two fully connected layers we take the widths for some positive integer N ∈ N and connection matrices Here the matrix Ξ D 2 , u takes a block form as Our first main result, to be proved in Section 5, can be stated as follows. A positive parameter τ (which can be arbitrarily small) is needed due to the continuous embedding of the Sobolev space W s 2 (S d−1 ) with s > d−1 2 into C(S d−1 ) found in Proposition 1 below.
Then for any f ∈ W r ∞ (S d−1 ), there exists a deep neural network consisting of J layers of CNNs with filters of length S and bias vectors satisfying (2.1) followed by downsampling and two fully connected layers with widths (2.2) and connection matrices (2.3) such that the hypothesis space H J,D 1 ,D 2 ,S contains a functionf satisfying where C r,d,τ,S is a constant independent of J or f . Moreover, the total number of free parameters N in the network can be bounded as Remark 1. To achieve the approximation accuracy f −f ∞ ≤ ǫ, we only need to take where we denote the smallest integer greater than or equal to u > 0 as ⌈u⌉.
When 0 < r < d − 1, we may set τ = 1 and see that the network achieving the approximation accuracy ǫ has depth J = O ǫ − 2d r and N = O ǫ − 2d r free parameters. When r > d − 1, we may set 0 < τ < r − (d − 1) and find the depth and the number of free parameters of this network are both of orders O (ǫ −2 ), slightly better than that in [32].
Our second main result aims at further demonstrating the superiority of deep CNNs over fully connected networks observed empirically in many practical applications. Motivated by our earlier work [33] on approximating ridge functions of type g(y · x) with y ∈ R d , g : R → R, and additive models (see, e.g., [4,24]) in statistics of the form f (x) = d j=1 g j (x j ) with univariate functions {g j } d j=1 , we consider a family of additive ridge functions of the form with y j ∈ S d−1 , g j : R → R for j ∈ {1, . . . , m}. The following theorem to be proved in Section 5 is about approximating additive ridge functions by deep CNNs followed by downsampling and one fully-connected layer. For 0 < α ≤ 1, denote the space of Lipschitz-α functions on for some 0 < α ≤ 1, then there exists a deep neural network consisting of J layers of CNNs with filters of length S and bias vectors satisfying (2.1) followed by downsampling and one fully connected layer h (J+1) with width (2N + 3)⌊(d + JS)/d⌋ and connection ma- The total number of free parameters N in the network can be bounded as

Remark 2.
To achieve an accuracy ǫ > 0 for approximating an additive ridge function (2.4) on S d−1 , we only need to take The total number of free parameters of the achieving network is of orders O ǫ − 1 α . This is the complexity required by the classical literature on fullyconnected networks to achieve an accuracy ǫ > 0 for approximating univariate functions. This extends our earlier work [33] from the ridge case with m = 1 to an additive ridge case with m ∈ N and hints the superiority of deep CNNs in approximating multivariate functions of special structures.

Spherical Analysis for Deep CNNs
In this section, we introduce ideas of our analysis before proving our main results. We first give a brief review on relevant concepts from spherical harmonic analysis and introduce some classes of functions. More details can be found in [6,29].

Spherical harmonics and Sobolev spaces on spheres
A spherical harmonic of degree n ∈ Z + on S d−1 is the restriction to S d−1 of a homogeneous and harmonic polynomial of total degree n defined on R d . Let H d n denote the set of all spherical harmonics of degree n on S d−1 . It can be found in [6] the dimension of the linear space H d n is The spaces H d n , n ∈ Z + , of spherical harmonics are mutually orthogonal with respect to the inner product of L 2 (S d−1 ). Since the space of spherical polynomials is dense in L 2 (S d−1 ), each f ∈ L 2 (S d−1 ) has a spherical harmonic expansion: is an orthonormal basis of H d n , f n,l is the Fourier coefficients of f given by and Proj n is the orthogonal projection of L 2 (S d−1 ) onto the subspace H d k of spherical harmonics, which has an integral representation: It can be readily shown that Z n (x, y) is the reproducing kernel of H d n independent of the choice of {Y n,l } N (n,d) where C λ n (t) is the Gegenbauer polynomial of degree n with parameter λ > −1/2, see, for instance, [6].
The spaces H d n of spherical harmonics can also be characterized as eigenfunction spaces of the Laplace-Beltrami operator ∆ 0 on S d−1 . Indeed, where λ n = n(n + d − 2) and C 2 (S d−1 ) denotes the space of all twice continuously differentiable functions on S d−1 . As a matter of fact, we may define the fractional power which is a self-adjoint operator on L 2 (S d−1 ) in the sense that is a Hilbert space with the inner product: f n,l g n,l .
In addition, we have the following continuous embedding lemma, see [11] and also [14, Eq. 14, p. 420]. By these lemmas, we know that the infinity norm can be bounded by the Sobolev norm when r > d−1 p .
where c r,d is a constant independent of f .

Near-best approximation and discretization
The best approximation of a function by those from polynomial spaces of various degrees might be nonlinear. A useful tool in spherical harmonic analysis is a linear scheme L n .

6)
where with λ = d−2 2 , l n is a kernel given by It can be found in [6,Chapter 4] that L n is near best, achieving the order of best approximation for f ∈ W r p (S d−1 ).
where c is a constant depending only on the function η in defining L n .
Note that since (−∆ 0 here and in what follows, we denote In fact, ζ n,r is a polynomial of degree at most 2n with expression The novelty here using a fractional power of (−∆ 0 + I) caused by the regularity of f ∈ W r ∞ (S d−1 ) enables us to get an r-dependent error bound for discretizing L n (f ): the larger the regularity index r, the smaller the bound.
To approximate the function L n (f ) by a neural network, we need a stepping stone, discretizing the integral form (3.6) to an empirical version given in terms of a sample y = {y 1 , . . . , y m } ⊂ S d−1 . The following estimate for the error L y n,m (f ) − L n (f ) will be proved by a probability inequality in Section 5. Such a probabilistic argument has been applied in [15].
, then for any n, m ∈ N, there exist y = {y 1 , y 2 , . . . , y m } ⊂ S d−1 such that where for θ ∈ R and n ∈ N we denote and C r,d,τ is a positive constant depending on r, d, τ but not on n, m or f .

Approximating ridge functions by deep CNNs
The last step in our spherical analysis of deep CNNs is to approximate the ridge function L y n,m (f ) by functions from the network with a bound to be proved in Section 5.

11)
where c ′ r,d is a constant depending only on r and d. The total number of free parameters N in the network can be bounded as

Comparison with Related Work
In this section, we give a brief review of related work on rates of function approximation by neural networks.
The fully connected shallow nets (1.1) or multi-layer nets (1.2) have nice approximation properties due to the fully connected nature, which was well studied in a large literature around 30 years ago. When the activation function is a C ∞ sigmoidal type function, approximation rates were obtained by many authors. In particular, in [1] rates were given for functions in f ∈ L 2 (R d ) whose Fourier transformsf satisfy a decay condition R d |w||f(w)|dw < ∞. Another typical result based on localized Taylor expansions asserts [19] that even for shallow nets (1.1) , if for some b ∈ R and some integer ℓ ∈ N \ {1}, the C ∞ activation function σ satisfies σ (k) (b) = 0 for all k ∈ Z + and lim u→−∞ σ(u)/|u| ℓ = 0 and lim u→∞ σ(u)/u ℓ = 1. These conditions required by the localized Taylor expansion approach are not satisfied by ReLU, so the approximation theory developed 30 years ago does not apply to ReLU. The difficulty was overcome in the recent deep learning literature and approximation properties of ReLU nets were established in [15] for ReLU shallow nets and functions satisfying R d |w||f(w)|dw < ∞, in [30,2,22,20] for deep nets and functions from W r ∞ ([−1, 1] d ) with 0 < r ≤ 2, and in [25] for approximation on manifolds. These results are obtained for fully connected nets.
Deep CNNs are different from fully connected nets. They have special sparse convolutional connection matrices (1.4), which leads to sparsity and reduces the computational complexity for structured data. Recently in [32], for functions f on Ω ⊂ [−1, 1] d satisfying f = F | Ω with F ∈ W r 2 R d and an integer index r > 2 + d/2, it was shown that the approximation accuracy f −f ∞ ≤ ǫ can be achieved by a deep CNN of depth 4⌈ 1 ǫ 2 log 1 ǫ 2 ⌉ and at most ⌈ 75 ǫ 2 log 1 ǫ 2 ⌉d free parameters. The linear increment of the free parameter number with respect to d improves the bound in Theorem 1 of [30] which requires at least 2 d ǫ −d/r free parameters and C 0 d 4 (log(1/ǫ)+d) fully connected layers with C 0 > 0 to achieve the same approximation accuracy ǫ. Periodized deep CNNs with different architectures and connection matrices different from the Toeplitz convolutional ones (1.4) were shown in [21,23] to be able to realize the output layer of any fully-connected DNN with free parameters of the same order. The same result was shown for deep CNNs (1.5) in [33].
The index r > 2 + d/2 required in [32] can be very large for processing high dimensional data. Hence approximated functions are required to possess high regularity which is not the usual case in applications. The problem happens because the approximation considered in [32] is measured in the L ∞ norm, while the Sobolev space H r (R d ) requires derivatives of various orders to belong to the L 2 space. This essentially causes the technical difficulty by To overcome the difficulty, we consider the case when the data is from the unit sphere S d−1 . The restriction can be relaxed in this situation through applying spherical harmonic expansions to construct a near best approximation for functions f ∈ W r ∞ (S d−1 ) in L ∞ norm, while the Sobolev embedding theorem is only used in discretizing integrals. In the literature, there have been some other harmonic analysis approaches in dealing with approximation by fully connected neural networks, using ridgelet transforms in [27], local Taylor expansions in [30], and B-spline functions in [28]. Our spherical harmonic analysis approach makes full use of the inner product nature (3.3) of the reproducing kernel of H d n which, after discretizing the polynomial approximation L n (f ) to L y n,m (f ), enables us to represent the linear transformations { y i , x } by deep CNNs with linearly increasing widths, an idea borrowed from our earlier work [32]. A key consequence of our approach is to allow the index r here to be an arbitrarily small positive number, which relaxes the restriction in [32] for the regularity of the approximated function. While the approximation of non-smooth functions is unified for r > 0 and the same order O (ǫ −2 ) for the number of network free parameters to achieve an approximation accuracy ǫ > 0 is kept when r > d − 1, a parameter number of order O ǫ − 2d r is required when 0 < r < d − 1. This is due to our approach of using a Hilbert space W s 2 (S d−1 ) in our probabilistic estimate for the discretization, which makes our rate suboptimal compared with [30,20,24] for r < d − 1. It would be interesting to derive optimal rates of approximating by deep CNNs functions from W r ∞ (S d−1 ) with small r. On the other hand, as stated in Remark 2, for an additive ridge function on S d−1 in the family (2.4), deep CNNs followed by a fully connected layer can extract linear features {y j } m j=1 and then approximate the function efficiently, with the same order of network free parameters as that for approximating a univariate function by fully connected DNNs. This demonstrates the superiority of deep CNNs in approximating functions with structures. It would be of great interest to explore other structures of multivariate functions for which deep CNNs together with network architectures like pooling and parallel channels may have super performance in function approximations and representations. Applying deep CNNs to some practical or empirical problems involving additive ridge functions (2.4) would also help understand advantages of convolutional structures of deep learning in some practical domains.

Proof of the Main Results
This section is devoted to the proof of our main analysis. Our analysis for the error f −f is carried out by means of the bounds for f − L n (f ) ∞ in Lemma 1, L n (f )− L y n,m (f ) ∞ in Lemma 2, and L y n,m (f )−f ∞ in Lemma 3.

Proving the lemma on discretization
To complete our analysis, we first prove Lemma 2 and Lemma 3. The proof of Lemma 2 is based on the following probability inequality for random variables with values in a Hilbert space which can be found in [26].
be independent random drawers of ρ. Then for any 0 < δ < 1, we have with confidence 1 − δ, Proof of Lemma 2. Recall that L n (f ) is defined by (3.6) with f ∈ W r ∞ (S d−1 ) and τ > 0. In applying Lemma 4 we take the Sobolev space W s 2 (S d−1 ) with the smoothness index s = τ + d−1 2 to be the Hilbert space H and the random variable ξ on S d−1 , σ d with values in H given by Then E(ξ) = L n (f ) and 1 m m i=1 ξ(y i ) = L y n,m (f ). To bound the norm ξ = ξ W s 2 , we recall the norm of W s 2 (S d−1 ) given by (3.4) with p = 2 and find for y ∈ S d−1 , . Then by the orthogonality and reproducing properties, where we have used the identity Z k (y, y) = N(k, d) found in [6] as Corollary 1. 2.7 and N(k, d) is the dimension of spherical harmonics H d k . Notice that for k ∈ N, k 2 < 1 + λ k ≤ dk 2 . We find (1 + λ k ) s−r ≤ d max{s−r,0} k 2(s−r) for either s − r ≥ 0 or s − r < 0. Since 0 ≤ η(t) ≤ 1 for t ∈ [0, 2], we can apply (3.1) to estimate the summation as Combining this with the definitions of the norm f W r ∞ (S d−1 ) and the function Λ τ (n), we know that ξ(y) 2 W s 2 can be bounded as This verifies (3.10) by the embedding Proposition 1 with p = 2 and s = τ

Proving the lemma on ridge approximation
The proof of Lemma 3 about approximating the function L y n,m (f ) is conducted by approximating the ridge functions l n ( y i , x ) in (3.9) with y i ∈ S d−1 ⊂ R d by deep CNNs. A key idea in our analysis is to use the inner product form (3.3) of the reproducing kernel of H d n and then to apply convolutional factorizations to realize the generated linear features, which enables us to conduct analysis after removing the restriction on large regularity index r. This idea might be applied to some other learning theory problems [8,10,17,34].
We first apply the following two lemmas proved in [32] implying that the linear function y i , x = y i · x can be realized by deep CNNs by factorizing y i regarded as a sequence into convolutions of filters supported in {0, 1, . . . , S}.
We then construct a fully connected layer to approximate the univariate function l n by continuous piecewise linear functions (splines) spanned by {σ(· − t i )} N i=1 with t i = −1 + i−2 N , based on the following well known result in approximation by splines which can be found in [7] and [31, Lemma 6].
Now we construct bias vectors in the neural networks. We denote for j = 2, · · · , J. The bias vectors satisfy b Hence the components of h (J) (x) satisfy where B (J) = Π J p=1 w (p) 1 . Applying the downsampling operator (1.6) leads to Hence d ≥ m.
We turn to expressing the last two fully connected layers. Of them, h (J+1) is given by with the connection matrix F ( where t := {t 1 < · · · < t 2N +3 } is given in Lemma 7. Note that F (J+1) is a determined matrix without free parameters. Then the first fully-connected layer h (J+1) (x) ∈ Rd (2N +3) of the deep network is for j = 1, . . . , m, while the other blocks are zero vectors.
Take the vector Θ N ∈ R 2N +3 in the connection matrix F (J+2) = Ξ T D 2 ,Θ N of the second fully-connected layer stated in (2.3) in terms of the linear operator L N as , then by the identity (5.4), for j = 1, . . . , m, the jth entry of the product The other entries of the product F (J+2) h (J+1) (x) vanish. Thus, by taking we see from the homogenous property σ(u/N) = σ(u)/N that the last layer h (J+2) is given by which is bounded by c ′′ r,d Λ d−1−r (n) with a positive constant c ′′ r,d depending only on r and d. It follows that ω(ζ n,r , 1 N ) ≤ c d c ′′ r,d n 2 Λ d−1−r (n)/N. Combining this with (5.9), Lemma 7 and the embedding Proposition 1 yields

General error bounds
With the proved bounds for f − L n (f ) ∞ in Lemma 1, L n (f ) − L y n,m (f ) ∞ in Lemma 2, and L y n,m (f ) −f ∞ in Lemma 3, the following bounds for the error f −f follows immediately.

Proving the main results
We are in a position to derive our main results from the general error bounds in Theorem 3.
The proof of Theorem 1 is complete.
Remark 3. When r = d − 1, from the above proof, we can see by taking N = ⌊n d+1 log(n + 1)⌋ that the statement of Theorem 1 still holds except that the bound for the number of free parameters should be replaced by N ≤ (3S + 3) J + 2J L t (g j ) ( y j , x ) .
Combining this with the additive ridge form (2.4) of f and Lemma 7, we know that for x ∈ S d−1 , Then the desired error bound is verified. The total number of free parameters N in the network is the sum of J(S + 1) contributed by w, J(2S + 1) by the bias vectors in the first J layers, 1 by the parameter B (J) in the fully-connected layer, and 2N +1 by the vector {g j (t i )} 2N +2 This proves Theorem 2.

Conclusion and Discussion
In this paper spherical harmonic analysis is conducted rigorously for the approximation theory of deep CNNs followed by downsampling and one fully connected layer or two on spheres. Our analysis provides rates of uniformly approximating functions f ∈ W r ∞ (S d−1 ) with r > 0 by deep CNNs followed by two fully connected layers. To approximate a Lipschitz function in a special additive ridge form, a network with one fully connected layer can be as fast as one for approximating a univariate Lipschitz function, which demonstrates the super power of deep CNNs in approximating or representing functions with special structures. Our spherical analysis relies on a special property of the reproducing kernel of H d n on the sphere. It would be interesting to extend our technique to approximation of non-smooth functions on [−1, 1] d and to L p approximation with 1 ≤ p < ∞.


extfieldsofstudy:
Computer Science
Engineering

text:
BBAND Index: A No-Reference Banding Artifact Predictor

Banding artifact, or false contouring, is a common video compression impairment that tends to appear on large flat regions in encoded videos. These staircase-shaped color bands can be very noticeable in high-definition videos. Here we study this artifact, and propose a new distortion-specific no-reference video quality model for predicting banding artifacts, called the Blind BANding Detector (BBAND index). BBAND is inspired by human visual models. The proposed detector can generate a pixel-wise banding visibility map and output a banding severity score at both the frame and video levels. Experimental results show that our proposed method outperforms state-of-the-art banding detection algorithms and delivers better consistency with subjective evaluations.

INTRODUCTION
Banding/false contour remains one of the dominant artifacts that plague the quality of high-definition (HD) videos, especially when viewed on high-resolution or Retina displays. Yet, while significant research effort has been devoted to analyzing various specific compression related artifacts [1], such as noise [2], blockiness [3], ringing [4], and blur [5], less attention has been paid to analyzing banding/false contours. Given the rapidly growing demand for HD/Ultra-HD videos, the need to assess and mitigate banding artifacts is receiving increased attention in both academia and industry.
Banding appears in large, smooth regions with small gradients and presents as discrete, often staircased bands of brightness or color as a result of quantization in video compression. All popular video encoders, including H.264/AVC [6], VP9 [7], and H.265/HEVC [8] can introduce these artifacts at lower or medium bitrate when coding contents containing smooth areas. Fig. 1 shows an example of banding artifacts exacerbated by transcoding. Traditional quality prediction algorithms such as PSNR, SSIM [9], and VMAF [10], however, do not align well with human perception of banding [11]. The development of a highly reliable banding detector for both original user-generated content (UGC) and Work done during an internship at Google.  transcoded/re-encoded videos, would, therefore, greatly assist streaming platforms in developing measures to avoid banding artifacts in streaming videos. Related Work. There exists some prior study relating to banding/false contour detection. Some methods [12][13][14] take advantage of local features such as the gradient, contrast or entropy to measure potential banding edge statistics. However, methods like these generally do not perform very well when applied to assess the severity of banding edges in video content. Another approach to banding detection is based on pixel segmentation [11,15,16], where a bottom-up connected component analysis is used to first detect uniform segments, usually followed by a process of banding edge separation. These methods are often sensitive to edge noise, though. We do not include block-based processing, as in [17,18], since it is hard to classify blocks where banding and textures coexist. If post-filtering is applied to these blocks, textures near the banding may become over-smoothed.
Our objective is to design an adaptive blind processor which can detect or enhance both "noisy" banding artifacts that arise in original UGC videos, as well as "clean" banding edges in transcoded videos. In this regard, it could be utilized as a basis for the development of pre-processing and Pre-processing  post-processing debanding algorithms. More recent banding detectors like the False Contour Detection and Removal (FCDR) [14] and Wang's method [11] are not designed for this practical purpose, and hence it is essential to devote more research to developing other adaptive banding predictors applicable to pre-or post-debanding implementations.
In this paper we propose a new, "completely blind" [19] banding model, dubbed the Blind BANding Artifact Detector (BBAND index), by leveraging edge detection and a human visual model. The proposed method operates on individual frames to obtain a pixel-wise banding visibility map. It can also produce no-reference perceptual quality predictions of videos with banding artifacts. Details of our proposed banding detector are given in Section 2, while evaluation results are given in Section 3. Finally, Section 4 concludes the paper.

PROPOSED BANDING DETECTOR
A block diagram of the first portion of the proposed model, which generates a pixel-wise banding visibility map (BVM), is illustrated in Fig. 2. Based on our observation that banding artifact appears as weak edges with small gradient (whether "clean" or "noisy"), we build our banding detector (BBAND), by exploiting existing edge detection techniques as well as certain visual properties. A spatio-temporal visual importance pooling is then applied to the BVM, as shown in Fig.  3, yielding "completely blind" banding scores for both individual frames and the entire video.

Pre-processing
We have observed that re-encoding videos at bitrates optimized for streaming often exacerbates banding in videos that already exhibit slight banding artifacts that may be barely visible, as shown in Fig. 1 § . We thereby deployed self-guided filtering [20], which is an effective edge-preserving smoothing process, to enhance banding edges. We deemed the guided § The exemplary frames used in this paper are from Music2Brain (YouTube channel). Website: https://bigweb24.de. Used with permission. to be a better choice than the bilateral filter [21], since it better preserves gradient profile, which is a vital local feature used in our proposed framework. Image gradients are then calculated by applying a Sobel operator after pre-smoothing, yielding a gradient feature map.

Banding Edge Extraction
Inspired by the efficacy of using the Canny edge detector [22] to improve ringing region detection [23], we performed a similar procedure to extract banding edges. After pre-filtering, the pixels are classified into three classes depending on their Sobel gradient profiles: pixels having Sobel gradient magnitudes less than T 1 are labeled as flat pixels; pixels with gradient magnitudes exceeding T 2 are marked as textures. The remaining pixels are regarded as candidate banding pixels (CBP), on which the following steps are implemented to create a banding edge map (BEM). (We used (T 1 , T 2 ) = (2, 12)). 1. Uniformity Check: Only the CBPs whose neighbors are either flat pixels or CBPs are retained for further processing. 2. Edge Thinning: Non-maxima suppression [22] is applied to each remaining CBP along its Sobel gradient orientation to better localize the potential bands. 3. Gap Filling: If two candidate pixels are disjoint, but able to be overlapped by a binary circular blob, the gap between the two points is filled by a proper banding edge. 4. Edge Linking: All connected CBPs are linked together in lists of sequential edge points. Each edge is either a curved line or a loop. 5. Noise Removal: Linked edges shorter than a certain threshold are discarded as visually insignificant. 6. Edge Labeling: The resulting connected banding edges are labeled separately, defining the ultimate BEM. The colored edge map in Fig. 2 shows a BEM extracted from an input frame. The banding edges are well localized.

Banding Visibility Estimation
Staircase-like banding artifacts appear similar to Mach Bands (Chevreul illusion), where perceived edge contrast is ex- Visual Importance Pooling  Fig. 3: Flowchart of the second portion (Section 2.4) of the proposed BBAND model, which produces banding scores on both frames and whole videos. aggerated by edge enhancement by the early visual system [24]. Explanations of the illusion usually involve the center-surround excitatory-inhibitory pooling responses of retinal receptive fields [25]. Inspired by the psychovisual findings in [26], we developed a local banding visibility estimator based on edge contrast and perceptual masking effects. The estimator processes the BEM and yields an element-wise banding visibility map (BVM).

Basic Edge Feature
Banding artifact presents as visible edges. As described earlier, we use the Sobel gradient magnitude as an edge visibility feature. Since edge visibility is also affected by content, we also model visual masking as it may affect the subjective perception of banding.

Visual Masking
Visual masking is a phenomenon whereby the visibility of a visual stimulus (target) is reduced by the presence of another stimulus, called a mask. Well-known masking effects include luminance and texture masking [23,27]. Here we deploy a simple but effective quantitative model of the effect of masking on banding edge visibility. Local Statistics: At each detected banding pixel in the BEM, compute local Gaussian-weighted mean and standard deviation ("sigma field") on the original un-preprocessed frame: where (i, j) are spatial indices at detected pixels in the BEM with corresponding original pixel intensity I(i, j), and w = {w k, |k = −K, ..., K, = −L, ..., L} is a 2D isotropic Gaussian weighting function. We use the µ(·) and σ(·) feature maps to estimate the local background luminance and complexity. The window size in our experiments was set as 9×9.
Luminance Masking: We define a luminance visibility transfer function (VTF ) to express luminance masking as a function of the local background intensity. We have observed that banding artifacts remain visible even in very dark areas, so we only model the masking at very bright pixels. A final luminance masking weight is computed at each pixel as where µ(i, j) is calculated using (1). (α, β) is a pair of constants chosen to adjust the shape of the transfer function. We used (α, β, µ 0 ) = (1.6×10 −5 , 2, 81) in our implementations. Texture Masking: We also define a texture visibility transfer function (VTF t ) to capture the effects of texture masking. The VTF t is defined to be inversely proportional to local image activity [23] when an activity measure (mean "sigma field") rises above threshold λ 0 . The overall weighting function is formulated as and where σ(i, j) is given by Eq. (2), and γ is a parameter that is used to tune the nonlinearity of VTF t . The values of (γ, λ 0 ) = (5, 0.32) were adopted after careful inspection. Cardinality Masking: The authors of [11] have shown that edge length is another useful banding visibility feature in a subjective study. We accordingly define the following transfer function which weights banding visibility by edge cardinality: where E(i, j) = {E ∈ BEM|(i, j) ∈ E} is the set of banding edges passing through location (i, j), and c 0 is a threshold on minimal noticeable edge length, above which banding edge visibility is positively correlated to normalized edge length. M and N denote the image height and width, respectively. We used parameters (c 0 , η) = (16, 0.5) in our experiments.

Visibility Integration
The overall visibility of an artifact depends on the visual response to it modulated by a concurrency of masking effects.
Here we use a simple but effective product model of feature integration at each computed banding pixel to obtain the banding visibility map: where w(·)'s are the responsive weighting parameters that scale the measured edge strength (Sobel gradient magnitude) |G(i, j)| at location (i, j).

Making a Banding Metric
Previous authors [27][28][29][30] have studied the benefits of integrating visual importance pooling into objective quality model, generally aligning with the idea that the overall perceived quality of a video is dominated by those regions having the poorest quality. In our model, we apply the worst p% percentile pooling to obtain an average banding score from the extracted BVM, where p = 80 is employed in experiments. Banding usually occurs in non-salient regions (e.g., background) while salient objects catch more of the viewer's attention. We thereby use the well-known spatial information (SI) and temporal information (TI) to indicate possible spatial and temporal distractors against banding visibility. SI is computed as the standard deviation of the pixel-wise gradient magnitude, while TI as the standard deviation of the absolute frame differences on each frame [31]. These are then mapped by an exponential transfer function to obtain weights: Finally, we construct the frame-level BBAND index by applying visual percentile pooling and SI weights to BVM: where K p% is the index set of the largest p th percentile nonzero pixel-wise visibility values contained in the BVM of frame I. We also obtain the video-level BBAND metric by averaging all frame-level banding scores, weighted by perframe TI, respectively: w TI (TI n ) · Q BBANDI (I n ). (10) Fig. 3 shows the entire workflow of the BBAND indices.

SUBJECTIVE EVALUATION
Other implemented parameters in our proposed BBAND model are (a SI , b SI , a TI , b TI ) = (10 −6 , 3, 2.5 × 10 −3 , 2), respectively, after empirical calibration, and we've found these selected parameters generally perform well in most cases. We evaluated the BBAND model against two recent banding metrics, Wang [11] and Baugh [16], on the only existing banding dataset, created by Wang et al. [11]. It consists of six clips of 720p@30fps videos with different levels of quantization using VP9. The Spearman rank-order correlation coefficient (SRCC) and Kendall rank-order correlation coefficient (KRCC) between predicted scores and mean opinion scores (MOS) of subjects are directly reported for the evaluated methods. We also calculated the Pearson linear correlation coefficient (PLCC) and the corresponding root mean squared error (RMSE) after fitting a logistic function between MOS and predicted values [32]. Table 1 summarizes the experimental results, and Fig. 4 plots the fitted logistic curves of MOS versus the evaluated banding models. These results have shown that the proposed BBAND metric yields highly promising performance regarding subjective consistency.

CONCLUSION AND FUTURE WORK
We have presented a new no-reference video quality model called the BBAND for assessing perceived banding artifacts in high-quality or high-definition videos. The algorithm involves robust detection of banding edges, a perceptioninspired estimator of banding visibility, and a model of spatial-temporal visual importance pooling. Subjective evaluation shows that our proposed method correlates favorably with human perception as compared to several existing banding metrics. As a "completely blind" (opinion-unaware) distortion-specific quality indicator, BBAND can be incorporated with other video quality measures as a tool to optimize user-generated video processing pipelines for media streaming platforms. Future work will include further improvements of BBAND by integrating with more temporal cues, and its applications to address such banding artifacts via debanding pre-processing or post-filtering.


Development and reconstruction of district water supply systems

The aim of the study is to develop a methodology for the comprehensive optimization of promising schemes of water supply systems, ensuring the selection of optimal water withdrawals taking into account water quality, structure and parameters of transporting facilities, as well as operating modes. An algorithm of theoretical, practical and methodological recommendations on optimizing promising schemes of water supply systems is presented. A new optimization technique for group water supply systems is proposed, which allows the complex to solve the following problems: to optimize the location and performance of water intake structures from surface and underground water sources; to optimize installation sites, productivity and composition of water treatment facilities; choose the route of group pipelines and their optimal diameters.

Introduction
Our country has huge water resources. However, due to their uneven distribution in many regions, there is an acute shortage of water, and where it is available, its quality does not always meet the requirements of SanPiN 2.1.4.1074-01 and GOST "drinking water". As a result of man-made impacts, the hydrological situation deteriorates every year, water sources are polluted with untreated sewage, and due to significant fluctuations in the level of surface and underground water, many water intakes operate on the verge of failure or are completely exposed. One of the ways to solve this problem is the construction and development of group water supply systems (WSS). A large number of such systems have already been built in the southern and central regions of our country. Many of them are complex, capital-intensive water management systems designed to transport and supply water over long distances to heterogeneous and dispersed consumers over a wide area. When designing and developing such structures, the most important tasks are: justification of the location and performance of water intake and treatment facilities, parameters of pumping stations and control tanks; selecting the optimal route and diameter of the pipeline system, taking into account the increased reliability, manageability, and seismic stability of individual structures and WSS as a whole.
optimal allocation of investments in the reconstruction of existing and construction of new group water supply facilities.
To date, we have accumulated a wealth of experience in the design and construction of WSS, developed various methods for justifying their structure, parameters and modes of operation. However, all of them do not allow us to fully solve these problems in a comprehensive manner and therefore require further development and development of new approaches. The existing methodology for optimal design is based on a comparison of cost, technical, reliability, and environmental characteristics of the analyzed options for the construction and operation of WSS. The main technological parameters of the designed water supply systems are taken into account through their cost in the optimization criteria, which is used to select the best option for the subsequent construction and operation of the WSS. Reliability, manageability and seismic stability are taken into account in the form of restrictions of a certain type, both on technical parameters and on investments in the construction and operation of WSS.

Methods
A common criterion for optimization and comparison of options is the cost of construction and operation of the WSS (by one year), which is determined by the following process [1]: where З прgiven the costs; Кinvestments in the network; Еthe coefficient of efficiency of investment, which in the conditions of market economy is identified with bank interest; ЭЗannual operating costs (thousand rubles/year), defined as follows [2,3,4]: where Самamortization charge (takes 0.05*K through pipelines, 0.09*K for pumping stations and surface water intake, 0.14*K for wells and underground water intakes); Ск.рthe cost of capital repairs (0.046*К), rubles/year; Ст.рcost of current repairs (0.01*К),rubles/year; Сэлкcost of electricity, thousand rubles/year: (3) where элкunit cost of electricity (rub per kWh), х -the amount of conveyed water (m 3 /s), Нthe pressure pumping station (m. the height of a column.). Сфзпpayroll of service personnel per 1 km of water pipe length: where ЗПср -average projected salary for the company, thousand rubles/month.; х, system performance, in m 3 /s. Ссвconsists of insurance contributions to the pension Fund (0.22*Сфзп), social insurance fund (0.029*Сфзп) and the health insurance fund (0.051*Сфзп); Сн.вthe tax on water (thousand rubles/year) is determined: where Сн.стtax rate in rubles for 1 m 3 waters, collected from surface and underground sources (in Irkutsk in 2018, this rate is 0.122 rubles per 1 m 3 ) ; Спрother expenses, thousand rubles/year, accepted in the amount of 20% of the amount of depreciation charges (Сам) and the cost of staff salaries (Сфзп), Taking into account (2), the operating costs will be as follows: To specify one-time (capital) costs, we will use the information provided in the consolidated construction price standards [5,6]. For outdoor water supply networks with a length of 1 km made of polyethylene pipes, when developing dry soil to a depth of 3 m (which is typical for the Irkutsk region), it is not difficult to obtain the following dependence of investment on the diameter of the pipe wire: ; vspeed in m/s; xflow rate in m 3 /s.
(8) Calculation of energy costs requires separate consideration. In the theory of hydraulic circuits [7], the law of energy conservation for an arbitrary hydraulic circuit is formulated as follows: "All the energy brought in, minus the energy for pouring water, is spent on overcoming friction forces": where J1a set of nodes of pumping stations, J2a set of nodes of water consumption, Hjpiezometric heads of pumping stations, Pjpiezometric water heads of consumers, in m. hihead losses on sections of the water supply network, xi -the flow of transported water on the network sections, Qj.water consumption of the consumer j, m 3 /s. Therefore, the annual cost of electricity for the water supply system can be represented as the sum of the cost of electricity for each section of the network: , for х in m 3 /s. For simplicity of presentation and calculation of pressure losses along the length of the pipeline we use the formula of F. A. Shevelev for plastic pipes [8]: where Уihydraulic gradient, Linetwork segment length, in m.
, or taking into account (7): The given estimated costs for the section of the network i in thousand rubles.: As shown by the research, the optimization criterion in the form of the estimated costs shown is mainly aimed at minimizing one-time capital investments. Although it is known that the operating costs for the life cycle of water disposal systems exceed one-time investment by 10-20 times. Obviously, when justifying the project for the reconstruction and development of the WSS, it is necessary to take into account the costs of the entire life cycle of the system (LCC). For a WSS consisting of n sections, these costs over the lifetime of the system can be represented as the following terms [9][10][11]: where Т -the service life of the water supply system in years is assumed to be equal to the service life of the longest-lasting element of the water supply system, Кinvestments equal to one-time costs, plus the cost of restoring individual elements of the system whose life span is less than the estimated (multiple costs).
Taking into account the values of operating costs (4), the life cycle costs (in thousands of rubles) can be represented: or taking into account (8) As a result, the given costs and life cycle costs are obtained as a function of the flow rate and speed of the transported water, the length of the pipeline and the specific cost of electricity. SP 31.13330.2012 and a number of guidelines recommend that the speed of water movement be taken in the range of 1.5-3 m/s. However, in each case, the optimal speed can be determined by calculation. To do this, you need to take a partial derivative of the reduced costs and life cycle costs by speed and equate it to zero. Then the resulting equation is resolved relative to the speed.
To select the optimal speeds, we will conduct numerical experiments on the reduced costs and life cycle costs, varying the unit values of electricity costs and pipeline lengths. The unit cost of electricity will change from the lowest for the Irkutsk region -1.06, to the maximum for Chukotka -8.2 rubles per kWh. the results of numerical experiments are presented in table. 1, which implies that the electricity tariff has a significant impact on the optimal speed values. The length of the networks practically does not affect the speed of water movement. The life cycle cost criterion reduces the optimal speed but 30%. Therefore, optimal speeds can be taken in the range from 1 to 3 m/s, depending on the cost of electricity in the region for which the water supply system is being designed.
Consider the construction region of the Irkutsk region. The length of the water supply network is 1 km made of polyethylene pipes. Let's construct the functions of the reduced costs as a function of water flow, while varying the speed from 0.5 to 3.5 m/s. These functions and their linear approximations are presented in picture 1.

Picture 1. Dependence of the reduced costs on the flow rate and speed water flows through pipelines and their linear approximations
For life cycle costs, these dependencies are also shown in picture 2. Similar calculations and transformations can be made for underground water intakes with treatment facilities. In recent years, with the advent of new technologies for water treatment, pumping equipment, block configuration of treatment facilities and pumping stations, their pricing policy and efficiency have changed significantly, expensive auxiliary rooms and enclosing structures are no longer required, and maintenance personnel and operating costs have been significantly reduced due to automatization. At the same time, the configuration and cost of treatment facilities depend significantly on the quality of the source water and, therefore, must be determined for each specific case.
Energy costs will be as follows: The estimated costs for network section i, in thousands of rubles: Consider a sewage treatment plant for iron removal and softening (removal of calcium and magnesium salts) produced by CC "Rosa" -typical water treatment stations based on pressure filters for de-Ironing and softening. The principle of degreasing-softening is based on catalytic oxidation and filtration of impurities in the loading layer. К = 178085 . х Operating cost: In this case, the energy costs will be as follows:  the problem of optimizing the structure and parameters of group water pipes as a problem of finding the maximum flow of the minimum cost [13,14,15]: when nnumber of network sections, Сi -cost per unit of flow. For example, at a speed of 3 m/s for the reduced cost criterion, the specific flow costs will be as follows: for new pipeline sections of the network 4381,0 thousand rubles per 1m 3 /s on 1 km length; for existing pipeline sections of the network (operating costs) 2463.0 thousand rubles per 1m 3 /s on 1 km length; for the device of wells 21224.0.0 thousand rubles per 1m 3 /s; for underground water intakes with water treatment system 1165500.0 thousand rubles per 1m 3 /s. To solve this problem, the following method is proposed in this paper [16,17,20]. A redundant scheme is being constructed, which is formed on the basis of the overlay of several pre-developed options for the construction of the WSS. The redundant scheme can represent: all possible places for water intakes (from open and underground sources) with local and centralized water treatment facilities, and without them (if the sources have drinking water quality); all possible routes for passing and supplying water to consumers by pressure pipeline sections.
The transport network is built on the basis of the received redundant scheme. To build it, all waterconsuming nodes are closed using fictitious branches to a common node tflow outlet, and watersource nodes are closed to a common fictitious node s-flow inlet. at the same time, restrictions on their throughput (upper and lower) are assigned for each section of the transport network. For fictitious stream entry branches, the upper limits correspond to the maximum possible water intake from water sources, and for fictitious stream exit branches, the upper limits correspond to the water needs of localities or subscribers. Restrictions on the flows of projected pipeline sections are not assigned. For existing sections of the network, the upper limit on flows is determined based on the optimal values of water flow rates: For each section of the network, the cost per unit of flow is determined and assigned. For fictitious branches of input streams, this is the cost of intake and water treatment of one m 3 of water, for real pipeline sections of the network, this is the cost of transporting one m 3 of water per 1 km of pipewater. For fictitious sections of the flow outlet, costs are not assigned, or are accepted as costs for further transportation of one m 3 of water to specific subscribers. For existing sections of the network, the flow unit costs correspond to operating costs.
Taking into account the transport network built in this way the problem of finding the maximum flow of the minimum cost is solved.

Results and discussions
As a result of optimization, the location and performance of water intake and treatment facilities, the optimal route of the pipeline system and its parameters are determined. At the same time, there are options for a fully centralized or decentralized water supply system.
Based on the proposed method, optimization calculations were performed and research was performed on the example of the Cheremkhovsk group water pipeline in the Irkutsk region. As possible sources of water supply, the following are considered: a new surface water intake from the


s2fieldsofstudy:
Computer Science
Mathematics
Engineering

text:
On detectability of cyber-attacks for large-scale interconnected systems ?

The paper deals with the key problem of detecting cyber-attacks in the context of large-scale systems (LSS). As these systems grow in size and complexity, cyber-attacks may target limited parts of the LSS, leading to tackling the problem in a decentralized way. We analyze the properties of distributed detection schemes under both local and interconnection attacks, and show that they are vulnerable to attacks that exploit the structure of the interconnections between subsystems. We also provide conditions and strategies that may be adopted to make the distributed control architecture regulating the LSS robust to this class of cyber-attacks.

INTRODUCTION
Many large-scale critical infrastructure and industrial control systems can be usefully described by the class of largescale dynamic systems (Lunze, 1992), in turn modeled as the interconnection of smaller subsystems. Indeed, this partition into subsystems may be done either to satisfy certain implementation objectives (such as a reduction of computational or communication costs associated with a centralized controller design), or because the plant itself is distributed over a large area. This has lead over the past decades to formulate distributed architectures to operate LSS (Šiljak and Zečević, 2005), where each subsystem is equipped with a local controller.
In order to achieve global objectives, controllers can share information over a communication network, allowing for coordination. In addition, these architectures increasingly rely on the integration of networked sensors and other socalled cyber infrastructure, leading to an increased concern in their cyber-security (Cárdenas et al., 2008). It has been ascertained in the security literature (Teixeira et al. (2015); Dibaji et al. (2019) and references therein) that malicious attackers can detrimentally affect the performance of a system by injecting disruptive signals into the plant's actuator, sensor, or communication channels, as real instances of attacks have proven (see Lee et al. (2016); Sobczak (2019) for recent evidence). Given this risk, together with the critical nature of LSS, it is deemed necessary to ensure their safety and security. Specifically, implementation of distributed monitoring architectures has attracted recent attention in literature (Anguluri et al., 2018;Gallo et al., 2018;Barboni et al., 2019), for which the distributed fault detection and isolation literature has been foundational (see for example Shames et al. (2011); Blanke et al. (2016); Boem et al. (2017a) to name but a few).
In this paper, we address the distributed detection of attacks on the communication network between controllers, named interconnection attacks, and highlight the structural limitations caused by the interconnection between subsystems. Furthermore, we include local attacks in our analysis, as in certain cases they can be modeled as interconnection attacks. The property of detectability of interconnection attacks is shown to be related to that of input observability of neighbors' states through interconnections, which can be seen as input to a subsystem's dynamics. Comprehensive results on input observability can be found in Hou and Patton (1998), or more recently Boukhobza et al. (2009) tackled a similar problem in a distributed context from a graph theoretic point of view.
Other works relevant to detectability of attacks are Pasqualetti et al. (2013) and Weerakkody et al. (2017). In the cited works, attackability is studied from a centralized point of view, and properties are formulated in terms of left-invertibility, which may not hold for subsystems with weakly coupled interconnections in LSS. Hence, we offer the following contributions: a. a generalized model of attacks which are local to one subsystem but target either the local control loop or the communication links;

Notation
Given a matrix A ∈ R p×n , rank A = p, A † is its the Moore-Penrose pseudo-inverse. I is the identity matrix of dimensions appropriate to the context. For a given matrix M , [M ] is used to define its related structural matrix, where specific entries are fixed at zeros, i.e., if [M ] (ij) = 0 then M (ij) = 0 also, while all others are set as free parameters. For some matrix Q, ker Q ⊥ denotes the space orthogonal to ker Q, the null space defined by Q. This decomposition is complete, hence any vector x =x + x ⊥ , wherex ∈ ker Q and x ⊥ ∈ ker Q ⊥ . Throughout the paper the superscripts a and h respectively indicate the attacked and nominal components of some signal x, i.e. in x = x h + x a , x h denotes its healthy component and x a the attacked one. With col[·], row[·], diag[·] we intend column and row concatenation of vectors or matrices, and block-diagonal concatenation of matrices.

Modeling Large-Scale Systems
We consider a large-scale system as a network of N subsystems S i , each of which is interconnected with a set N i ⊆ N {1, . . . , N }, called the set of neighbors of S i . The dynamics of the state of the LSS can be partitioned into N equations, each describing the dynamics of S i as: where x i ∈ R ni , u i ∈ R mi , and y i ∈ R pi are respectively the state, the control input, and the output of S i . The term j∈Ni A ij x j models the interconnection between different subsystems. Matrices A ii , B i , A ij , and C i are of appropriate dimensions. Furthermore, without loss of generality, rank C i = p i . Assumption 1. The pair (A ii , C i ), ∀i ∈ N is observable.
Assumption 2. For some full column rankÃ ij such that Im(Ã ij ) = Im(row[A] ij∈Ni ), there are no zero dynamics fromÃ ij to C i . Remark 3. In this paper we focus on the structural properties of a distributed diagnostic tool. As such no uncertainties are included in (1).
We consider that controllers C i are designed to implement a distributed control algorithm, for which u i , i ∈ N is computed both from locally measured outputs y i , and information received from neighboring controllers C j , j ∈ N i . In order to implement such an algorithm, it is necessary to introduce a peer-to-peer communication network between the controllers, as represented in Figure 1. Assumption 4. The topology of the network of physical interconnections between subsystems is mirrored by that of the communication network, and the networks are undirected.
We consider two scenarios: a) Subsystems share their local measurements y i . b) Subsystems share locally computed state estimatesx i .
The state estimatex i of x i is such that the estimation error i x i −x i → 0 asymptotically as t → ∞. How to compute such an estimate is out of the scope this paper. We introduce the following vector to model the information received by C i : wherex i col[x] j∈Ni is a vector grouping all the states of the neighboring subsystems S j , j ∈ N i , acting as an external input to S i . Notation (2) allows to consider the two information exchange scenarios in a unified way.
In the following, with some abuse of notation, we refer to its j-th block asC i,(j) . Assumption 5. If the following condition does not hold: controllers C i , i ∈ N transmit their local estimatex i . Remark 6. We note here that Assumption 5 is required to allow for the possibility of obtaining a local estimatex i such that i → 0 as t → ∞ without making the estimate independent of the interconnections.
To simplify notation, from (2) we can define an augmented output vector y i , containing both local measurements and transmitted signals:

Modeling cyber-attacks
Definition 7. (Local attack). Attacks between the controller C i and plant S i . These attacks alter the control input signals u i and the output measurements y i . Definition 8. (Interconnection attack). Attacks on information transmitted between controllers C j , or some subset of the neighboring subsystems j ∈ N i ⊆ N i , and C i . The attacker alters the signals contained inỹ i .
These two classes of attacks presuppose different architectures, and different attack resources. Specifically, local attacks presume that the control loop of S i is closed over a communication network, and the attacker is capable of corrupting the signals exchanged between C i and the networked sensors/actuators. On the other hand, in interconnection attacks, it is the communication network over which the controllers exchange information that is compromised. In Figure 1 we show how these different classes of attacks may affect the system. Note that we do not consider the possibility of attacks on the physical interconnections of the LSS.
Following the framework described in Teixeira et al. (2015), we model cyber-attacks in communication networks as additive signals. Specifically, the state dynamics, measurement and received communication outputs in (1) and (2) are defined as: where Moreover φ i models an interconnection attack onỹ i received by C i , µ i a local attack on the actuators, and η i one on the sensor measurements, as represented in Figure 1. We consider both B a i and D a i to be such that Im(B a i ) ⊆ Im(B i ) and Im(D a i ) ⊆ Im(C i ). Note that in (5) we have combined the modeling framework necessary to implement either local or interconnection attacks. Appropriate definition of u a i allows for modeling a number of different attack scenarios (Teixeira et al. (2015)). We denote the starting time of the attack as T a . Assumption 9. For t ≥ T a , an attack input function u a i = 0, for a single i ∈ N . Furthermore, it is defined either as Remark 10. Assumption 9 implies that at any given time either (u i , y i ) orỹ i can be trusted, i.e. are not corrupted by a malicious attacker.

Distributed attack detection architectures
Given the possible vulnerability of the LSS's communication networks, it is necessary to equip it with a monitoring layer. Specifically, as shown in Figure 1, we consider the addition of local model-based diagnostic modules D i capable of detecting local and interconnection attacks.
We choose to consider distributed architectures, as they are capable of integrating properties of scalability, low computational complexity, and privacy amongst subsystems within the diagnosis layer of the LSS. As such, we are interested in the analysis of its properties, and specifically of those limitations that arise following the partition of the LSS into subsystems, while having constraints on the information available to each local diagnoser. Let us define the set of information used in the design of D i : Definition 11. (Information set). The information set I i {M i , γ i , υ i } represents the information available to the diagnoser D i , where M i is a tuple defining the dynamics, and γ i and υ i are output and input signals. Definition 12. (Reduced-information distributed architecture). We define an architecture as "reduced-information" if the information set of D i is When D i is designed according to this definition, it is said to be "reduced-information" itself.
As D i does not have direct access to x i and x j∈Ni , to achieve detection the diagnoser defines an appropriate residual r i given I i . This residual is such that in nominal conditions r i → 0 as t → ∞, while if under attack this is detected once r a i is larger than some suitably chosen threshold.
Differently from other works in literature, we are not here focused on the detectability properties of the diagnoser with respect to local attacks. Rather, we are interested in the analysis of stealthiness of interconnection attacks with respect to the structure of the partitioned LSS. To clarify the explanation throughout the remainder of the paper, let us define the following classes of attacks 1 : Definition 13. (Stealthy attacks). An attack is stealthy if no monitor detects the attack, i.e. for an attack on S k , u a k = 0, r a i = 0, for any k ∈ N and all i ∈ N . Definition 14. (Locally stealthy attack). An attack is locally stealthy if [µ i , η i ] = 0 implies r a i = 0. Definition 15. (Structurally stealthy attacks). An interconnection attack on S i is said to be structurally stealthy to D i if, given I i , φ i = 0 leads to r a i = 0.

DETECTABILITY LIMITATIONS OF REDUCED-INFORMATION D i
In this section, we highlight the limitations of D i with respect to the partitioned structure of the LSS, given a reduced-information distributed architecture. We focus our attention to those attacks that are not locally detectable, namely, interconnection attacks and locally stealthy attacks. The first are analyzed as they influencẽ y i as received by C i , not visible to D j , j ∈ N i , which are transmitting the information. Of the latter, we focus on local covert attacks, inspired by covert attacks in Smith (2015), introduced in a distributed scenario by Barboni 1 The classes of attacks defined here can be seen as subsets of undetectable attacks in Pasqualetti et al. (2013).
Preprints of the 21st IFAC World Congress (Virtual) Berlin, Germany, July 12-17, 2020 et al. (2019), and shown to be locally stealthy to D i , as summarized in the following Lemma: Lemma 16. (Barboni et al. (2019)). Given a local covert attack active for time t ≥ T a , where µ i is freely chosen by the attacker to achieve some objective, and η i is defined as the output of the following linear time invariant system with initial condition x a j (T a ) = 0: a diagnostic module D i is not capable of detecting it.
Although Lemma 16 holds for all local diagnosers D i , in Barboni et al. (2019) it was shown that, through suitable design of the residual generator, it is possible for diagnosers of neighboring subsystems D j to detect local covert attacks. As such, they share similar properties to interconnection attacks. Proposition 17. Given a local covert attack on S j such that Lemma 16 holds, it is equivalent to µ j = η j = 0 and an interconnection attack In (8), with abuse of notation, φ i,(j) are the components of φ i corresponding to the signal received by C i from C j . Remark 18. Although in this paper we limit our analysis to local covert attacks, for both space restrictions and ease of explanation, the method exploited in Proposition 17 can be applied to any locally stealthy attack.
As such in the remainder, we consider that for some unknown j ∈ N i for t ≥ T a , φ i,(j) = 0, modeling both interconnection and local covert attacks on S j . We now introduce the main result of this paper. Theorem 19. Given a subsystem S i with dynamics as in (5), equipped with a detector D i designed with information set I i , suppose that, starting at some time T a > 0, the subsystem is exposed to an interconnection attack φ i,(j) = 0, for some j ∈ N i . Such an attack is structurally stealthy to local diagnoser D i if and only if Proof.
[Sufficiency]: In order to detect an attack, D i must generate a residual such that, given an attack oñ y i,(j) , ifC j x j =ỹ i,(j) , then r i → 0.
As the control input u i is known to D i (given I i in (6)) and as it cannot be altered maliciously (Assumption 9), its effect on y i will be disregarded. An appropriate estimator can be defined such that, nominally, r i → 0 is satisfied. Consider the map from x j to y i : Given that x j and M j are unknown to D i , the neighbor's state must be estimated (statically) in order for r i → 0 to be satisfied in nominal conditions: such a static estimate isx j (ỹ i,(j) ) =C † jỹi,(j) . Since rankC j = p i ≤ n i , this introduces error ε i,(j) (ỹ i,(j) ) x j −C † jỹi,(j) = 0. However, this does not affect the residual generator in nominal behavior, as jC j x j ∈ kerC j ⊆ ker A ij , given Assumption 5, and therefore r i → 0, even though ε i,(j) = 0. Hence,x j (ỹ i,(j) ) is an appropriate estimate of x j . Following the results on input observability in Hou and Patton (1998) any componentx j = 0,x j ∈ ker(A ij ) does not affect y i , as x j is in fact not input observable from y i . To prove sufficiency, suppose that φ i,(j) is designed by the attacker such that (9) holds. Due to linearity, the estimated output affected by the attack,ŷ a i , is: Hence, for t ≥ T a , r a i = 0, asŷ a i = 0, ∀t ≥ T a , and as such the attack is stealthy.
[Necessity]: We prove necessity by contradiction. Let us suppose thatC † j φ i,(j) ∈ ker A ij , but that the attack is indeed structurally stealthy.
The residual generator comparing y i andŷ i (ỹ i,(j) ), as t → ∞: Given thatC † j is full column rank by design, and that condition (9) is assumed not to hold,C † j φ ⊥ i,(j) = 0. Furthermore, following Assumption 2, there are no zero dynamics betweenÃ ij and C i . Hence, r a i = 0, which contradicts the hypothesis.
The result of Theorem 19 shows some similarities with respect to the class of zero-dynamics attacks in literature. Indeed, asx i can be considered as an exogenous input to S i , φ i can in turn be seen as an attack on the measurements of the inputs, and the class of attacks defined in Theorem 19 can be interpreted as a subset of zero-dynamics attacks. However, rather significantly, the attacks satisfying the stealthiness condition in Theorem 19 exploit the fact that A ij is not full column rank, rather than requiring knowledge of the dynamics of S j to exploit its transmission zeros. The above result is also relevant because attacks in this scenario are acting upon S i through the computation of u i by C i , which proves fundamental when considering the design of control architectures robust to such attacks.
If we consider an attacker with knowledge of the structure of the interconnection between subsystems, rather than A ij itself, a special case of Theorem 19 is derived. Such a scenario is relevant as not only said structure might be easier to obtain, but also because, in fact, subsystems in LSS are often described as being weakly coupled (Lunze, 1992). As a consequence, the interconnection matrices are often not structurally full rank, and the following holds. Corollary 20. If rank[A ij ] < n j , it is sufficient that be satisfied for attack to be structurally stealthy.
Note furthermore that, with the appropriate model knowledge and disruption resources, it would be possible for an interconnection attack to satisfy condition (9) even if rank A ij = n j . Indeed, given an attack on a number of interconnections such that j∈ Ni n j > n i is satisfied, and given knowledge of matrices A ij , for j ∈ N i , an attacker would be able to satisfy Theorem 19. A similar result was presented in Boem et al. (2017b), where the properties of a specific distributed detection architecture were analyzed.

CONTROLLER ROBUSTNESS TO STRUCTURALLY STEALTHY ATTACKS
In this section, we are interested in analyzing the interaction between structurally stealthy attacks and distributed control architectures. By affecting only the received measurementsỹ i , interconnection attacks φ i do not have a direct impact on the dynamics of S i . However, through appropriate definition of φ i , the control signals u i may be altered, thus maneuvering the state of S i to an operating point undesired by the operator. It is therefore important to define control laws such that u i is robust to structurally stealthy interconnection attacks. Definition 21. Consider a LSS with an interconnection attack φ i that satisfies Theorem 19, for some φ i,(j) , j ∈ N i . Let C i compute u i κ i (y i ), for some operator κ i 2 . We say that the control architecture is robust if ∀i ∈ N , for all attacks φ i,(j) = 0 stealthy to D i .
In order to ensure control architecture C i is robust, there are two strategies that can be pursued: • ensure κ i ([0 ,ỹ i ] ) = 0 for allỹ i such thatC † jỹi,(j) ∈ ker A ij ; • augment information set I i so that the class of stealthy interconnection attacks is reduced.
This section is dedicated to the former, while Section 5 will address the latter. Before proceeding, we analyze the special case in which φ i,(j) models a local covert attack on S j . Indeed, although in Proposition 17 we have shown that they can be interpreted as interconnection attacks, they affect the state of S j directly, and could therefore have an impact on the behavior of the LSS in two ways: i. by propagating the "attacked state" to other subsystems; ii. by having neighbors' controllers react to behavior that is not nominal. We address these cases separately. Proposition 22. If S j is subject to a local covert attack, controllers C i , i ∈ N , are robust to it.
Proof. In Proposition 17, we assumed the state of a subsystem under a covert attack to be nominal. This can be seen as the definition of a new nominal state and a shift of the attack action entirely onto the information transmission between subsystems. With respect to the actual nominal state from the perspective of C j , however, we have thatỹ i,(j) = y h j , which in fact satisfies (12), as φ i,(j) = 0. Proposition 23. If S j is subject to a local covert attack satisfying (9), the large-scale system composed of subsystems S k , k ∈ N \{j} is not influenced by the attack.
Hence, the dynamics of the large-scale system as a whole is, from (5): where A is a block-matrix where each ik-th entry is j∈N \{j} , and A kj groups the interconnection terms A ij for all S i such that j ∈ N i . We decompose x = x h +x a , where x a is the state affected by all attacks on the LSS. From Proposition 22, u a = 0, as for all i ∈ N j ,ỹ i,(j) = y h j . Hence, to prove x a = 0, i.e. that the state of the LSS once S j is subtracted is not affected by the attack, we must show that A kj x a j = 0. Recall condition (9) and that φ i,(j) = −C j x a j from Proposition 17. Thereforȇ C † jC j x a j ∈ ker A ij , implying that x a j ∈ ker A ij , ∀i ∈ N j . Hence, x a j ∈ ker A kj . As such, x a j will not influence the dynamics of any of its neighbors.
Having established that local covert attacks that are structurally stealthy do not propagate to the rest of the LSS, we evaluate the case in which there is an interconnection attack affectingỹ i . Consequently, the control inputs computed by C i should depend only on measurementsỹ i,(j) associated to those components of x j which directly influence S i through physical interconnection. This translates into ensuring that the arguments of the control function lie in ker(A ijC † ) ⊥ , as shown in the next result. Specifically, we introduce matrix Γ i diag[Γ] ij∈Ni , with Γ ij suitably designed for each j ∈ N i so that the following holds: Proposition 24. If each controller is such that with each Γ ij designed such that C i is robust to all structurally stealthy attacks.
Proof. By setting u i as in (15), it is clear that it is no longer required to detect φ i,(j) , but rather we are restricting the sensitivity of our residual to Γ ij φ i,(j) . As such, rather than (9), we are interested in verifying whether an attack function φ i,(j) such that C † j Γ ij φ i,(j) ∈ ker A ij holds. By setting Γ ij such that (15) holds, no attack function such that Γ ij φ i,(j) = 0 exists such that As such, any attack function such that κ i ([0 , Γ i φ i ] ) = 0 will not be stealthy to D i , and controllers are therefore robust to structurally stealthy attacks.
A possible definition of Γ ij such that Proposition 24 holds is to set it as a basis of ker(A ijC † j ) ⊥ . This is indeed quite conservative, requiring C i to only exploit measurements of neighbors' state components physically coupled with S i .

AUGMENTING I i FOR DETECTABILITY OF D i
A robust control law as presented in the previous section might be infeasible or in contrast with design constraints. As a result, in this section, we show how the addition of information on dynamics of neighboring subsystems S j , j ∈ N i to I i may assist with the design of robust controllers. Indeed, under certain conditions, this permits the estimation of x j directly, and to not rely on the input observability ofx i .
To illustrate the idea, we start by assuming that all of the model information M j , j ∈ N i is available to D i , showing how a residual generator may be capable of detecting structurally stealthy interconnection attacks. Given Assumption 1, it is possible to augment residual generator D i to include additional knowledge of the dynamics of S j .
By extending D i to detect additive attacks onỹ i,(j) , attackers must satisfy further constraints in order to be stealthy. Specifically, the following holds: Proposition 25. Attack φ i,(j) is stealthy to D i if and only if it is both structurally stealthy given A ij , and can be seen as a zero-dynamics attack on S j throughỹ i,(j) .
Proof. If D i is composed of r i and r i,(j) , j ∈ N i , it holds that for an attack φ i,(j) = 0 must satisfy both r a i = 0, r a i,(j) = 0, as if either of these do not hold D i will detect an attack. Following Theorem 19, φ i,(j) must be a structurally stealthy attack for r i → 0. As attack signal φ i,(j) is an additive signal toỹ i,(j) = y j , being a zero-dynamic attack is known to be a necessary and sufficient condition for it to be stealthy to a residual generator exploiting M j (Pasqualetti et al., 2013).
In general, there are conditions on the dynamics of S j that must be satisfied for D i to compute convergent estimates. While out of the scope of this paper to discuss such conditions in a general way, for completeness we refer the reader to Gallo et al. (2018). In that work, distributed diagnosers exploiting model knowledge of neighbors are designed to detect interconnection attacks in islanded DC microgrids.

CONCLUSIONS
Analysis of stealthy attacks in partitioned LSS have been presented, within the framework of interconnection attacks. We have derived conditions under which such attacks are not detectable, and have shown how controllers can be made robust to them, thus preventing the impact to be propagated to the entire LSS. Detectability, furthermore, can be recovered by properly extending the information set locally available to each diagnoser. Future studies will be focused on defining such a set so that it is minimal while ensuring estimation convergence.