s2fieldsofstudy:
Computer Science
Mathematics

text:
Sparsest Input Selection for Controllability of Singular Systems via a Two-Step Greedy Algorithm

In this paper, the problem of determining the sparsest input matrices to ensure controllability of linear singular systems is investigated. Firstly, it is shown that, determining the sparsest input matrices to ensure reachable controllability or complete controllability is NP-hard, even when the system ‘singularity’ is arbitrarily large. Secondly, submodular functions for singular systems are built, upon which greedy algorithms are developed to approximate the sparsest input matrices with guaranteed performance bounds for the case where there is no restriction on the number of independent inputs. Thirdly, a two-step greedy algorithm is proposed for determining the sparsest input matrices with a given number of inputs to ensure controllability. Compared with the existing algorithms for sparsest input selections, the proposed algorithm achieves better trade-off between the approximation performances and computation efficiency, which are demonstrated by two simulation examples.

I. INTRODUCTION
Input selections under the objectives to meet/optimize certain system performances have long been active but challenging issues in control community [1]. In a recent paper [2], it is first shown that, given an autonomous system as (1), it is NP-hard to determine the minimal number of state variables that need to be actuated by an input to ensure system controllability, where x(t) is the state vector, and A is state transition matrix. A more general proposition is that, given a collection of possible input matrix columns, it is NP-hard to choose the minimal number of columns to form an input matrix so that the resulting system can be controllable [2], [3]. An input selection problem for a more general class of linear systems is considered in this paper, which is known as the singular system or the descriptor system. Analysis and synthesize of descriptor systems have received significant attention due to their widespread application in modeling The associate editor coordinating the review of this manuscript and approving it for publication was Donatella Darsena . and control of variously actual systems, such as electrical systems, biological systems, economical systems, hybrid systems and intelligent transportation systems [4]- [8], [10], [11]. Additionally, the descriptor system involves the usual state-space formulation as a special form, and can establish some dynamical systems without state-space formulation. This paper focuses on determining the sparsest input matrices with a given number of inputs that ensure controllability of singular systems. Two points make this problem worthy of extra attention. One is that, criteria for complete controllability of singular systems are a little more complicated than that of the nonsingular systems described by the standard state-space model. Some related concepts, like modes, controllable subspaces, etc., are often defined upon the standard canonical forms of singular systems, which make extending input selection methods from nonsingular systems to singular ones not so straightforward. The other one is that, as argued in [9], the problem of determining the sparsest input matrices with guaranteed performance bounds to ensure controllability is far from completely settled for nonsingular systems with multiple eigenvalues. It is left for us to find algorithms with guaranteed bounds, acceptable VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ approximation performances and relatively low computational complexity. Related Work: Controllability and observability of largescale networked dynamical systems have drawn the attention of control scientists [12]- [14]. As mentioned earlier, determining the minimal number of actuated states to ensure controllability for nonsingular systems has been proven to be NP-hard [2]. However, there is a simple greedy algorithm which can maximize the rank increase of controllability matrix in each iteration. It can achieve a logarithmic approximation factor and can obtain the best performance in polynomial time. Moreover, it has a closed form solution with the minimum input number to ensure controllability, which is equal to the maximum geometric multiplicity of the state transition matrix [15]. Note that, this is similar to actuator and sensor deployment problems investigated under structural framework in [12], [16]. And more specifically, system matrices are fixed with zeros or free parameters. For example, the matching theory [12] is used to provide the minimum input number required for the controllability of structure, which enable the problem of determining the minimum actuated states with structural controllability guarantee in polynomial time [16].
Apart from focusing on the binary concept of controllability, researchers also proposed some energy related metrics to quantify controllability and found appropriate actuator selection strategies to satisfy certain performance criteria [17]. Approaches for input selection based on this concept, including submodular optimization [3], a method of the relaxed control energy metric [18], have been developed in the literature. These investigations extend the binary concept of controllability to quantitative one, which deepens our insight into controllability of singular system. However, these methods either have higher computing complexities or have lower approximation performance.
Concepts of reachable controllability and complete controllability for descriptor systems seem to be first developed in [19]. Afterwards, various criteria have been proposed, including some rank based criteria [19], graph characterization for structural controllability [20], and matroid theory based criteria for mixed descriptor systems which entries of system matrices are either fixed constants, or free parameters [21]. Input selections for controllability of mixed descriptor systems have been considered in [22], under a dimensionless assumption defined therein. The main difference between problems considered in that paper and this one is that, the objective in [22] is to determine the minimal number of dedicated inputs, that is, an input is dedicated if it actuates only one state, whereas we will show that for a given singular system with fixed system matrices, selecting the minimal number of dedicated inputs is less complicated than selecting the sparsest input matrices with limited number of inputs.
Main Contributions: In this paper, we consider the problem of determining the sparsest input matrices with given number of inputs to ensure controllability of singular systems.
The contributions of this paper are in three aspects. Firstly, it is proven that determining the sparsest input matrices to ensure reachable controllability or complete controllability is NP-hard, even when the system 'singularity' is arbitrarily large. Secondly, submodular functions for singular systems are built, upon which greedy algorithms are developed to approximate the sparsest input matrices with guaranteed performance bounds for the case where there is no restriction on the number of independent inputs. Thirdly, a two-step greedy algorithm is proposed for determining the sparsest input matrices with a given number of inputs to ensure (reachable or complete) controllability. This algorithm has much less search space than the simple greedy one, and remarkably, has a submodularity-ratio like provable approximation bound. Compared to the existing methods, numerous numerical simulations show that this algorithm achieves better trade-off between the approximation performances and computation efficiency.
The organization of this paper is as follows. Section II gives the problem formulation, and Section III presents some preliminaries. Section IV gives computational complexities of the problems considered in this paper. Afterwards, Sections V and VI deal with algorithms for problems considered in this paper. Some numerical examples are provided in Section VII. Section VIII finally concludes this paper.
Notations: R, C, Z and N denote the set of real, complex, integral and nonnegative integral numbers, respectively. For a matrix M , M ij denotes its (i, j)-th entry if no confusion is made. ||M || 0 denotes the number of nonzero entries in a matrix M . M † denotes the Moore-Penrose pseudo-inverse of M . Let J , J 1 , J 2 be a set of integers, then M J represents the submatrix of M comprising the columns with indices given by J , and M J 1 ,J 2 the submatrix of M comprising the rows indexed by J 1 and columns indexed by J 2 . |S| denotes the cardinality of a set S.
we denote the block diagonal matrix (compact matrix) with the i-th diagonal block (row block) being M i .

Consider the following linear time invariant (LTI) system
where x(t) ∈ R n is the state vector and u(t) ∈ R l is the input vector at time t, A ∈ R n×n , B ∈ R n×l are respectively the state transition matrix and input matrix. The constant matrix E ∈ R n×n may be singular. System (2) or (E, A, B) is called the linear singular system or descriptor system. Obviously, when E is nonsingular, this system collapses to the standard state-space system [23]. If det(sE − A) = 0 for at least one s ∈ C, where E may be singular, the pair (E, A) is called regular. System (2) is said to be solvable, if for any admissible inputs u(t) ∈ R l , there is a unique state solution x(t) satisfying Equation (2). It is known that system (2) is solvable, if and only if det(sE − A) = 0 for at least one s ∈ C. In this paper, we always assume that system (2) is solvable.
In this work, we are interested in determining the sparsest input matrices to ensure controllability of the linear singular systems. Before proceeding that, the notion of reachable controllability (R-controllability) and complete controllability (C-controllability) for singular systems is introduced.
The admissible initial state is the initial state vector x(0) ∈ R n such that Equation (2) admits a differentiable solution of x(t) and u(t). The set of all admissible states is denoted by R 0 , which is called the reachable set.
Definition 1 (R-controllability [19]): System (2) is said to be R-controllable, if for any two states x 0 and x 1 in the set of admissible initial states R 0 , there exists a finite time T and a differentiable control input u(t)| 0≤t≤T , such that Definition 2 (C-controllability [19]): System (2) is said to be C-controllable, if for any two states x 0 , x 1 ∈ R n , there exists a finite time T and a differentiable control input For notation simplicity, we just say (E, A, B) is R-controllable (C-controllable) if system (2) is so.

Lemma 2 (Conditions for C-Controllability
The problems considered in this paper can be formulated as following. Problem 1: Given the pair (E, A) in (2), determine Problem 2: Given the pair (E, A) in (2) and an integer l ≤ n such that there exists a B ∈ R n×l making (E, A, B) C-controllable (resp. R-controllable). Determine Problem 2 differs from Problem 1 in that, there is a restriction on the available number of independent inputs. Such scenario occurs, for example, in a leader-follower multi-agent system where the number of leaders is limited, or in a circuit system where the number of independent voltage sources is limited. It is shown in [9] that, these two problems are not necessarily equivalent for a nonsingular system where A has multiple eigenvalues. It is straightforward to see that, Problem 1 is always feasible, as the input matrix B = I n always makes (E, A, B) C-controllable and R-controllable. Besides, Problem 1 can be seen as a special case of Problem 2 by setting l = n.

III. PRELIMINARIES
This section presents some preliminaries for our further derivations.

A. GENERALIZED EIGENVALUES
The generalized eigenvalue problem of two matrices M , N ∈ R n×n is to find a scalar λ and the corresponding vector φ such that φM = λφN . We call λ the generalized eigenvalue and φ the corresponding left eigenvector of the pair (M , N ). There are many numerical algorithms to determine the generalized eigenvalues. See [24] and the references therein.
Using the generalized eigenvalue arguments, the following lemma states an equivalent condition for (3) to hold.
Lemma 3: Suppose that the solvability condition holds for system (2), i.e., det(sE − A) = 0 for at least one s ∈ C.
Assume that there are m distinct generalized eigenvalues For each s i , let X i be the matrix consisting of the maximum number of linearly independent row vectors of φ, where φ satisfies φE = s i φA (i.e., X i is a basis matrix of the left null space of s i E − A). Then, (3) holds, if and only if X i B is of full row rank for each i = 1, . . . , m.
Proof: See the appendix.

B. SUBMODULARITY
Let be a finite set and 2 be its the power set. A set function f : 2 → R assigns a real scalar to each subset of .
A nonincreasing function f is a set function such that for all

Definition 3 (Submodularity [25]):
A set function f : 2 → R is submodular if for any sets F 1 , F 2 with F 1 ⊆ F 2 and any element ω ∈ \ F 2 , we have that and supermodular if the reversed inequalities in (6) hold. Lemma 4 [25]: Suppose that the set function f : 2 → R is submodular and nonincreasing. For the optimization problem min F ⊆ |F| under the constraint f (F) ≥ α, where α is a given threshold, the simple greedy algorithm (i.e., in each iteration choosing the element with the maximum increase in f (F)) achieves a solution F gre that satisfies the following inequality where F * is the optimal solution, F T −1 is the set returned at the second-last iteration of the greedy algorithm.

C. MATROID
Matroid is a structure that captures and generalizes the notion of linear independence in vector spaces. Given a finite set and elements together with a family F = {F 1 , F 2 , · · · } of subsets of , the pair ( , F) is a matroid if: From the above definition, is also named as the ground set, and a member of F is referred to as an independent set. For a matrix E, a linear matroid can be defined as N (E) = ( , F), where is the set of indices of columns of E, F is the collection of indices of columns of E which are linearly independent, i.e., F = {J : rank(E J ) = |J |}. Given two matroids N 1 = ( , F 1 ) and N 2 = ( , F 2 ), the matroid intersection N 1 ∩N 2 is defined as the collection of all common independent sets of N 1 and N 2 . The cardinality of N 1 ∩ N 2 , represented by ρ(N 1 ∩ N 2 ), is defined as A structured matrix is a matrix whose entry is either a fixed zero or can take an arbitrary real value independently. The set of n × m structured matrices is denoted by {0, * } n×m , where * denotes the entry which can take values independently. For a matrix M whose entries are functions of some free parameters, its generic rank, indicated by grank(M ), is the maximum rank it can achieve as the function of its free parameters.
Lemma 5 [26]: The generic rank of XB can be computed via matroid intersection algorithm, which has polynomial time complexity. From Lemma 5, it is easy to see that grank(XB) = ρ(N (X ) ∩ N (B )).

IV. COMPLEXITY ANALYSIS
In this section we prove that both Problem 1 and Problem 2 are NP-hard for the singular systems. We will first give the computational complexities and then present their proofs.
Theorem 1: When the purpose is to make the system C-controllability, Problem 1 and Problem 2 are NP-hard, even if the rank deficiency of E is arbitrarily large.
Theorem 2: When the purpose is to make the system R-controllability, Problem 1 and Problem 2 are NP-hard, even if the rank deficiency of E is arbitrarily large.
Since by setting l = n, Problems 1 and 2 are equivalent, which are both equivalent to finding the sparsest diagonal input matrix to ensure system controllability, 1 for the proofs of these two theorems it suffices to prove the NP-hardness of the latter problem. In the following, we will use e p×q to denote a p × q matrix whose entries are all one, 0 n×n the n × n zero matrix, and e n i the ith standard basis vector in the n-dimensional vector space. We will say (E, A) is k C-controllable (R-controllable), if there exists a diagonal matrix B with no more than k nonzero entries such that (E, A, B) is C-controllable (R-controllable).
Our proofs adopt some similar ideas from [2]. The reductions start from the set cover problem.
where k is called the cardinality or size of this set cover. It is known that, given an integer k, determining whether there exists a set cover with size k is NP-complete. The complete proofs are given in the following.
Proof of Theorem 1: From this construction, matrix X I,I is block triangular with nonzero diagonal entries, and hence invertible, , such a i always exists. Let X ⊥ ∈ R (h+p+q)×q be a basis matrix that spans the right null space of X . Since the first h + p columns of X are linearly independent, it can be validated that Based on these constructions, it can be validated that XE = XA. Therefore, X i , which is the ith row of X , is the left generalized eigenvector of (E, A) associated with the generalized eigenvalue i. Moreover, there is a monomial s p+h det X ⊥ J 1 ,J 2 in det(sE − A) that cannot be zeroed out by other monomials. Hence, (E, A) is regular. As rank(E) = h + p, E is singular and its rank deficiency can be arbitrarily large. Obviously, the above mentioned A can be constructed in polynomial time by solving linear matrix equations. In the following, we will prove that the system associated with (E, A) is k + q C-controllable, if and only if the collection C has a set cover with size k.
For the one direction, assume that C has a set cover Let the input matrix B ∈ R (h+p+q)×(h+p+q) be such that B jj = 1 for j = i 1 , . . . , i k , and B jj = 1 for j = h + p + 1, . . . , h + p + q, and other entries are zero. Then, from Lemma 2, it can be seen that (E, A, B) is C-controllable.
For the other direction, assume that (E, A) is k + q C-controllable. This means that, there exists a diagonal matrix B ∈ R (h+p+q)×(h+p+q) with no more than k + q nonzero entries, so that (E, A, B) is C-controllable. According to Lemma 2, to make rank([E, B]) = h+p+q, the (h+p+1)th, (h + p + 2)th, . . ., (h + p + q)th diagonal entries of B must be nonzero. On the basis of that, to make (E, A, B) C-controllable for a diagonal B, it is necessary and sufficient to make the support 2 of the (h + i)th row of X intersects with J , for each i ∈ {1, . . . , p}, where J = {j : B jj = 0, j ∈ {1, . . . , p+h}} and |J | ≤ k. As every element of {1, . . . , p} is contained in C 1 ∪ · · · ∪C h , for each j ∈ J ∩ {h+1, . . . , h+p}, we can find a replacement j * ∈ {1, . . . , h}, such that C jj * = 1 and thus the support of the jth row of X intersects with J \{j} ∪ {j * }. As |J | ≤ k, we finally have that the collection C has a cover set with size no more than k. The NP-hardness of Problem 1 for regular singular systems then follows the NP-completeness of the set cover problem.
Proof of Theorem 2: Again we use the reduction from the set cover problem. The reduction is almost the same as that of the C-controllability case. The only difference is that, there is no need to consider the condition (5). A little bit modifications are sufficient to fit for this purpose. To this end, let E = diag{I h+p+1 , 0 q×q }, where q > 1 and the 'singularity' q can be arbitrarily large. Parameter matrix It can be verified that the first (h + p + 1) columns of X are linearly independent. Hence, X is of full row rank. Let E = diag{I h+p+1 , 0 q×q }. Afterwards, construct matrix A ∈ R (h+p+q+1)×(h+p+q+1) in the similar way to the C-controllability case. It can be verified that, (E, A) is regular. Then, following similar arguments to the proof of Theorem 1, it can be claimed that (E, A) is k R-controllable, if and only if the collection C has a set cover with size no more than k + 1. Remark 1: As mentioned above, our main ideas follow those of [2], but with some technical differences. More specific, to generalize their proofs to the singular case, as the rank deficiency of E (which can be seen as the 'singularity' of the system) can be arbitrarily large, it takes extra attention to guarantee that (E, A) is regular in the constructions. Besides, the conditions for C-controllability are more complicated than those of the Popov-Belevitch-Hautus (PBH) test for the nonsingular state-space systems.

V. SUBMODULARITIES FOR SINGULAR SYSTEMS AND SOLUTIONS TO PROBLEM 1
In this section, we identify some submodular functions for linear singular systems that facilitate approximation algorithms for Problem 1, as well as the two-step greedy algorithm in the next section for Problem 2.
First, we consider the dimension of reachable subspace. Unlike the linear nonsingular case, the reachable state set of linear singular systems is often hard to be represented in the original forms of their system parameters [19]. A transformation is utilized to obtain the standard canonical form of system (2). To be specific, there exist two nonsingular 2 The support of a vector is the set of indices of nonzero entries of this vector. matrices P and Q, such that where E 2 ∈ R n 2 ×n 2 is nilpotent, i.e., E k 2 = 0 for some integer k, and integers n 1 + n 2 = n. The standard canonical form of system (2) is given as follows: According to [17], the reachable subspace R rs of system (2) is given by where ⊕ denotes direct sum.
Proposition 1: Let dimR rs (B) be the dimension of reachable subspaces of system (2) with input matrix B. For any given B ∈ R n×l , dimR rs (B S ) is submodular on S ⊆ {1, . . . , l}.
Proof: See the appendix. However, utilizing the dimension of reachable subspaces based on the standard canonical form to approximate Problem 1, may be not natural and cause extra computational burdens. In the following, we give another submodular function based on the generalized eigenvectors.
Let m, s i | m i=1 and X i | m i=1 be defined as in Lemma 3. Moreover, let r i be the number of rows of X i . It is valid that The greedy algorithm based onf (S) (similar to Algorithm 1, the termination condition isf (S) = m i=1 r i ) achieves an O(log rank(E))-approximation for Problem 1 when the objective is to ensure R-controllability.
Proof: See the appendix. This section has identified some submodular functions which facilitate approximation algorithms for Problem 1 on singular systems. In fact, these results are consistent with nonsingular systems. Algorithms in this section are the first step towards the two-step greedy algorithm developed in the next section.

VI. TWO-STEP GREEDY ALGORITHM FOR PROBLEM 2
We now proceed with Problem 2 when there is a fixed number of available inputs. A two-step greedy algorithm with provable performance bounds will be proposed.

A. FEASIBILITY CONDITION AND NONSUBMODULARITY
For the ease of derivations, let r m+1 be the dimension of the left null space of E (i.e., r m+1 = n − rank(E)), and X m+1 be the matrix consisting of a collection of linearly independent row vectors which span the left null space of E. Then, from Lemma 3, system (2) is C-controllable, if and only if X i B is of full row rank for i = 1, . . . , m + 1.
The following proposition gives the feasibility condition for Problem 2. It is easy to see that P(i) has zero Lebesgue measure in R n×l max . Define a set P by P=R n×l max \∪ m+1 i=1 P(i). As m+1 is countable, P is open and dense in R n×l max . As a result, a real valued matrix always exists in P, which makes the associated system C-controllable.
Proposition 2 generalizes the known result that, the minimal number of inputs to guarantee controllability of a standard state-space described system equals the maximum geometric multiplicity of the system state transition matrix [13]. As a result of the above mentioned proposition, Problem 2 is feasible if and only if l ≥ max{r max , n − rank(E)}. We denote l max max{r max , n − rank(E)}. In the following, we always assume that Problem 2 is feasible, i.e., l ≥ l max .
It is natural to consider the following function: whereB ∈ {0, * } n×l , as a performance function 3 for It can be verified that, Problem 2 has the same optimal value as the following problem: In fact, although solutions to the above mentioned optimization problem are structured matrices, which characterize the zero-nonzero patterns of the input matrices, it has been proven in [11] that once we have the feasible zero-nonzero patternsB of the input matrices, a numerical realization ofB making each X i B full of row rank, thus making the system C-controllable, i = 1, . . . , m + 1, can be determined in polynomial time. Hence, in what follows, we mainly focus on the zero-nonzero patterns of the input matrix B, which is denoted byB.
With a little abuse of notation, givenB 1 ,B 2 ∈ {0, * } n×l , we sayB 1 ⊆B 2 , ifB 1ij = 0 impliesB 2ij = 0. ByB 2 \B 1 we denote a structured matrix whose (i, j)th entry is nonzero only ifB 2ij = * andB 1ij = 0. An elementē i,j ∈B 1 is such that e ∈ {0, * } n×l , and there is only one nonzero entry inē, which is the (i, j)th entry with (i, j) satisfyingB 1ij = 0. LetB full be the n × l structured matrix whose entries are all nonzero. It is easy to see that g(B) is a set function onB ⊆B full . 4 However, g(B) is usually neither submodular nor supermodular. This can be seen from the following example.

B. SUBMODULARITY-RATIO FAILS TO GIVE NON-TRIVIAL PERFORMANCE BOUNDS FOR G(B)
Recently, submodularity-ratio, which is used to measure how far a set function is from submodularity, has been adopted to give some performance bounds for greedy algorithms for nonsubmodular functions [27]. Definition 4 (Submodularity-Ratio, [27]): Leth be a nonnegative set function on V . For any S, T ⊆ V , defineh S (T ) = h(T ∪ S) −h(S). The submodularity-ratio ofh is defined as The submodularity-ratio γh is such that 0 ≤ γh ≤ 1, and is 1 if and only ifh is submodular. The γh usually gives some nontrivial performance bounds for the greedy algorithm of a nonsubmodular function, but fails to so do when γh = 0. 5 However, the following example shows that, g(B) indeed has a zero submodularity-ratio.
Example 2 Example 1 Continued, Zero Submodularity-Ratio of g(B): Consider the parameters in Example 1. It holds that g(

C. A TWO-STEP GREEDY ALGORITHM
We propose an alternative for the simple greedy algorithm, which is a two-step greedy algorithm. This algorithm is an overall greedy algorithm, but consists of two steps. The first step is to approximate a minimal set S such that (E, A, I S ) is C-controllable. Here S is the set of directly actuated states. The second step is to select elements from (B full ) S,C , where C = {1, . . . , l}, that is, to select elements from the rows of B full indexed by S. The detailed description is given in Algorithm 2. Remarkably, this algorithm has a submodularityratio like provable approximation bound, as shown in the following theorem.
Theorem 4: The two-step greedy algorithm achieves an O(l max log(n))-approximation for Problem 2.
To prove Theorem 4, we first prove the following two intermediate results.
The above lemma means that, before the termination of Algorithm 2, in each iteration when an element is added toB, the increase of the associated performance function (i. e., g(B)) is at least one.
Let us have a brief comparison between the two-step greedy algorithm and the simple greedy one. The worst-case approximation ratio O(l max log(n)) for the former algorithm has a form similar to the performance bound of using nontrivial submodularity-ratios for nonsubmodular functions. By contrast, currently we have not obtained a non-trivial performance bound for the simple greedy algorithm. Suppose that the second step of the two-step algorithm and the simple greedy algorithm both terminate after N ite iterations. Then, the ratio of size of total search spaces of these two algorithms is n|S|+|S|lN ite nlN ite , where S denotes the set returned in the first step of Algorithm 2. Usually |S| is much less than n [2]. Hence, this ratio is generally much smaller than 1. It is expected that, the two-step greedy algorithm is more computationally efficient than the simple greedy one. Our numerical experiments in the next section will show that, the two-step greedy algorithm almost has the same approximation performance as that of the simple greedy one, while is more computationally efficient.
Remark 2: By changing g(B) to beḡ(B) = m i=1 r i , Algorithm 2 can be directly modified to fit for selecting the sparsest input matrices with a given number of inputs to ensure R-controllability. Details are omitted due to its obviousness.

VII. NUMERICAL RESULTS
We present two numerical examples to show the main results of this paper.
First, we calculate a descriptor system defined by an electric circuit in [20]. In this example, the system dynamics are defined by Eẋ(t) = Ax(t) + Bu(t), where Bu(t) are the alternative input options, and the system state is defined by the voltage and current at each element of the circuit network of Fig.1. We have For each element of the series circuit, R 1 (R 2 ), L 1 (L 2 ), C are the resistance, inductance and the capacitance respectively. Let the parameters be R 1 = R 2 = 100, L 1 = L 2 = 10, C = 200, where their physical units are standard  for C-controllability. Substitute the above B 1 and B 2 into Lemma 1 and Lemma 2 respectively, by Lemma 3, it can be found that rank(s j E − A, B i ) = 8 for j = 1, 2, 3, i = 1, 2, and rank[E, B 2 ] = 8. This means that, the associated system is R-controllable under the input matrix B 1 (corresponding to putting a voltage source on the resistor R 1 ), and C-controllable under the input matrix B 2 (corresponding to putting a voltage source on the inductance L 2 and four current sources to control the currents i 1 , i 2 , i 3 and i 4 ).
Next,we present some numerical results to show the efficiency and effectiveness of the proposed algorithms in this paper. To compare our two-step greedy algorithm with algorithms in the existing literature, consider the problem of constructing the sparsest input matrices for controllability of nonsingular systems, i.e., E = I . We generate system state transition matrices using the same method as in [9], which can generate square matrices with multiple eigenvalues.
More precisely, the maximum geometric multiplicity of eigenvalues is set to be 3, which means l max = 3. For each dimension of states n, ranging from 10 to 400, 10 independent plants are generated. For each plant, set the number of inputs l to be 3, 4 and 5 respectively, and use the graph coloring based algorithm proposed in [9], Algorithm 2 and the simple greedy algorithm (based on g(B)), respectively, to determine the zero-nonzero structure of input matrices that ensure controllability. The averaged computation time, and the averaged sparsity of input matrices are given in Figs. 2, 3 and 4 respectively for l = 3, 4 and 5. All computations are implemented with a personal notebook computer with Intel(R) Core(TM) i7-5500U CUP and 12.0 GB RAM. The corresponding codes for numerical experiments could be found in the Github [28].
From Figs. 2, 3 and 4, the two-step greedy algorithm achieves comparable approximation performances to the simple greedy algorithms, while costs much less time than the latter algorithm. The graph coloring based algorithm in [9] has almost the same time consumption as that of the twostep greedy algorithm. However, when the system dimension increases or the number of independent inputs l gets closer to l max , its approximation performances seem not as good as those of the two-step greedy algorithm. It can be declared that, the two-step greedy algorithm may be preferable, as it achieves better trade-off between the approximation performances and time complexity.

VIII. CONCLUSION
The problem of determining the sparsest input matrices with given number of inputs to ensure controllability of singular systems is considered. It is proven that determining the sparsest input matrices to ensure R-controllability and C-controllability is NP-hard, even when the system 'singularity' is arbitrarily large. Because of lack of submodularity and a nonzero submodularity-ratio, it is hard to give a non-trivial approximation bound of the simple greedy algorithm for this problem. Therefore, a two-step greedy algorithm is proposed. This algorithm has much less search space than the simple greedy one, and remarkably, has a submodularity-ratio like provable approximation bound. Through numerous numerical simulations, it is shown that this algorithm achieves better trade-off between the approximation performances and computation efficiency compared to the existing methods. As for future work, it is interesting to derive some non-trivial approximation bounds for the simple greedy algorithm in solving the problems in this paper. Proof of Proposition 1: Partition P as P = col{P 1 , P 2 } with P 1 ∈ R n 1 ×n , P 2 ∈ R n 2 ×n . Then, According to Theorem 7 of [3], dim < A 1 |P 1 B S > and dim < E 2 |P 2 B S > are both submodular on S ⊆ {1, . . . , l}. Since the summation of two submodular functions is still submodular, dimR s (B S ) is submodular on S ⊆ {1, . . . , l}.
Proof of Theorem 3 and Corollary 1: From Proposition 2.1.9 of [26], given a constant matrix M with l columns, rankM S is submodular on S ⊆ {1, . . . , l}. Notice that the summation of two submodular functions is still submodular. Based on these facts, f (S) andf (S) are both submodular functions. Noticing that f (S) andf (S) are both non-decreasing, the performance bounds in Theorem 3 and Corollary 1 then follow from Lemma 4.


An Efficient Medical Image Retrieval and Classification using Deep Neural Network

Background/Objectives: The main objective of this work is to obtain an efficient brain tumor image retrieval and classification using Deep Neural Network (DNN).Methods/Statistical analysis: The features from the medical images are extracted by using tamura feature extraction, Local Ternary Pattern (LTP) and Histogram of Oriented Gradients (HOG). Subsequently, an Infinite Feature Selection (Inf-FS) technique is incorporated to select optimum features from feature vector, which leads to improve the classification process using sparse auto encoder based DNN. Furthermore, the retrieval performance of the proposed method is improved by Euclidean Distance technique. Findings: An Open Access Series of Imaging Studies (OASIS) and Contrast Enhanced Magnetic Resonance Image (CE-MRI) datasets are utilized to analyze the proposed method. The sparse auto encoder based DNN classification scheme yields an overall accuracy of 95.34% in OASIS dataset and 99.87% in CEMRI dataset with improved sensitivity, specificity, error rate. The retrieval performance of proposed technique is assessed in terms of Average Retrieval Precision (ARP) and compared with two existing methods such as Local Mesh Vector Co-occurrence Pattern (LMVCoP) and Content Based Image RetrievalConvolutional Neural Network (CBIR-CNN). The ARP of the proposed method for CE-MRI and OASIS dataset is 98.33% and 88.25% that is high when compared to the CBIR-CNN, LMVCoPmethod.Novelty/Applications: An appropriate feature selection using Inf-FS and DNN based nonlinear feature data classification are used in the applications of medical image retrieval.

• There are three texture feature extraction techniques such as Tamura, LTP and HOG that are used and then hybridized to enhance classification performance of MRI images. The main advantage of using Tamura is feature derivation based on the psychophysical context. Moreover, the HOG is chosen because of its geometric and photometric transformations and LTP is highly robust to noise. • Additionally, the optimum features are chosen from the feature vector using Inf-FS technique as it helps to assess important features along with its subsets of features from the images. These attributes are learned through DNN based classification scheme. • Moreover, the reduction of irrelevant features from the feature vector helps to improve the classification accuracy. The Euclidean based image retrieval technique is utilized to retrieve brain MRI images from a huge database.

Literature survey
This section provides a literature survey of recent trends in the medical image retrieval system. Generally, the process and analysis over the brain tumor of MRI images are challenging task for the computer vision scientists. Moreover, the selection of the brain tumor is also considered as a crucial step during brain tumor classification. On the other hand, there are various healthcare applications are developed such as classification of stomach infections (26) , recognition of blood cells (27) and Wireless Capsule Endoscopy for abnormality detection (e.g., Ulcer, gastric infection) (28) (29) (30) . In (31) developed a brain parcellation tool wherein all the anatomical images were converted into a standardized index through image transformation and atlas-based parcellation technique. The primary progressive aphasia was selected because of the apparent degeneration of cells at various degrees and locations in patients. Furthermore, the Partial Least Squares-Discriminant Analysis (PLS-DA) and Principal Component Analysis (PCA) were used to perform individual classifications over various medical images. Several critical information such as distance/ clinical severity was not included in the analysis.
In (32) presented the DNN and Extreme Learning Machine (ELM) based brain image classification using prominent features. The contourlet transform and Zernike moments were used to extract the texture and shape features respectively. Then the Genetic Algorithm (GA) and Particle Swarm Optimization (PSO) based feature selection was utilized to hybridize the features. The redundant features were removed by using the GA and PSO based feature selection. The DNN and ELM based methods classify only the brain tumor images from the database, it fails to classify the type of brain tumor.
In (33) presented the combination of Local Pattern Descriptor (LPD) and Gray-level Co-occurrence Matrix (GLCM) as a feature extraction approach to extract the features from medical imageries. This work implements a Local Mesh Vector Cohttps://www.indjst.org/ occurrence Pattern (LMVCoP) by integrating the Local Vector Co-occurrence Pattern (LVCoP) and Local Mesh Co-occurrence Pattern (LMCoP). An OASIS database consisting of MRI brain images was utilized to evaluate LMVCoP methods. In clinical applications, fast decision making was achieved based on LMVCoP based image retrieval. This method fails to retrieve in the extensive high dimensional datasets.
In (34) developed the Co-Active Adaptive Neuro Fuzzy Inference System (CANFIS) classifier and Support Vector Machine (SVM) for diagnosing the brain tumor. The features from LTP and GLCM were learned using SVM and CANFIS classifiers to categorize the brain images as normal and abnormal. From the classified abnormal image, the tumor region was segmented using normalized graph cut approach. The value present outside of the learning database was not effectively determined by CANFIS and a huge volume of database was required for the classification model.
In (35) presented the CBIR system for brain tumor images. The T1-weighted CE-MRI dataset images were taken for validating the image retrieval performance. Here, the feature extraction was performed using the deep Convolutional Neural Network (CNN) model (VGG19). Additionally, the retrieval performances were enhanced by using the transfer learning and block-wise fine-tuning strategy. The results show that this method works well only for small datasets, particularly in the medical imaging domain, but fails to achieve in large datasets.
In (36) presented the deep learning-based method for microscopic brain tumor detection and tumor type classification. The brain tumor was extracted by using the 3D Convolutional Neural Network (CNN) architecture of the two different datasets such as BraTS2015 and BraTS2018. Next, the extracted tumors were transmitted to the pretrained CNN model for extracting the features. An optimum feature was obtained by giving the extracted features to the correlation-based selection method. Hence the final classification was achieved based on the validation using feed-forward neural network. However, the classification time was increased for the 3D CNN based tumor segmented image.
In (37) developed the hybrid technique to extract the tumors the MRI image. This hybrid technique was categorized into five steps such as de-noising of an image, the extraction of the tumor, feature selection, feature fusion, and classification. Initially, the image denoising was carried out in two different datasets such as BRATS2013 and private dataset using the curvelet transformation. The ant colony optimization was used along with the Thresholding method to extract the tumor based MRI scans. A multi type features such as shape and texture were extracted from the tumor images. Next, the PCA reduced skewness method was used to reduce the irrelevant features. Further, the selected features were given to the SVM to classify the tumor portions. But, the SVM classifier was required 70% of features to classify the brain tumors.
In (38) presented the 3D reconstruction algorithm for CT image features by considering the multi-threaded deep learning calculations. Next, the segmentation of image feature texture was obtained using the DeepLabv2 algorithm. From the CT image sub-block content, an optimal feature volume data was acquired that used to solve the issue of packet loss while processing the feature extraction. Next, the Compute Unified Device Architecture (CUDA) was used to obtain the multi-thread calculation of CT image features that used to resolve the serial processing of the contents using multiple sub-blocks. This CUDA was created longer reconfiguration time when it processed with multiple slice CTs.
The machine learning method helps in extracting the features and it performs the classification. Accordingly, the complexity and time consumption was increased, when the extracted features were highly dimensional (39) (40) . The feature selection and reduction algorithms were beneficial for both the accuracy and execution time in medical imaging applications. Additionally, the selection of significant feature minimizes the overall complexity of the system (41) . Hence, the proposed uses the Inf-FS method to select the optimal features from the feature vector. This relevant feature selection from the feature vector is used to improve the classification accuracy of sparse auto encoder based DNN.

Proposed method
In the proposed method, an effective image retrieval over the medical images is performed using optimal feature extraction, DNN and Euclidean distance techniques. The steps processed in this medical image retrieval system are given as follows: 1) CLAHE based image enhancement, 2) Feature extraction, 3) classification using sparse auto encoder based DNN and 4) Euclidean distance based image retrieval. Figure 1 depicts the overall architecture of the proposed method.

Data acquisition
In the proposed method, two different databases such as OASIS (42) and CE-MRI (43) are used to analyze the image retrieval performance. Generally, the OASIS repository contains a collection of cross-sectional images under the age group of 18 to 96. Additionally, the CE-MRI dataset has three types of brain tumor images such as Glioma, Meningioma, and Pituitary tumor. The collected images from the database are subject to the CLAHE method for enhancing the image quality. The sample images from the OASIS and CE-MRI are shown in Figure

Image enhancement using CLAHE
The CLAHE method preprocesses all the collected images from the OASIS and CE-MRI dataset. The nonlinear property of the CLAHE method is used to modify the intensity values to enhance the contrast of images. The clipping limit of CLAHE is utilized to overcome the problem of noise amplification wherein histogram is clipped at predefined value to limit the amplification process. In this approach, an original image is separated as non-overlapping contextual regions and two main parameters such as the size of the block and clip limit are used to control the image quality. The process of CLAHE is described as follows: • Initially, the original image (I) is divided into non-overlapping contextual regions. In that case, the total amount of image tiles equals to M × N. • For every contextual region, histogram is calculated for the array image. • The clipped histogram and gray level mapping are generated in CLAHE. The pixels are equally separated at each gray level and the average number of pixels is given by equation (1) where, N avg represents the average amount of pixel and N gray represents the amount of gray level in the contextual region. The amount of pixels present in X and Y dimension are N r X and N r Y respectively. The clip limit is expressed in the equation (2).
where, the N CL represents the actual clip limit and N clip is the normalized clip limit in the range of [0,1]. For example, the preprocessed image using CLAHE is shown in Figure 3. The preprocessed images (I ′ ) are given to the feature extraction for subsequent process. https://www.indjst.org/

Feature extraction
In this proposed method, three different feature extraction methods such as tamura, LTP and HOG descriptor are utilized to extract the features from the preprocessed image. The tamura feature extraction is generally used for extracting the texture features which is used to describe the image textures based on the coarseness, contrast and directionality. The main advantage of the LTP is resilient to the noise and the image's orientation is based on the global features from the HOG descriptor. The combination of these features is essential to enhance the classification accuracy of the different classes of MRI brain images.

Tamura feature extraction
Tamura is basically a texture-based feature extraction utilized to represent the image content and it is extensively used in image searching and retrieval. In tamura, three features are extracted such as coarseness, contrast and directionality.

• Coarseness
The way of computing the size of texture element is known as coarseness and it can be expressed as: where, the difference of average intensity in horizontal and vertical directions are represented as k = argmax{E 1 , E 2 , .., E k , .., E L }. E and S best (i, j) = 2 k . Moreover, this average intensity E is calculated by using the grey level position of image I ′ .

• Contrast
The intensity difference at an image's local region is defined as contrast that is expressed in equation (4).
where, the standard deviation is represented as σ , the squared variance is represented as σ 4 and represented as k = argmax{E 1 , E 2 , .., E k , .., E L }. E and S best (i, j) = 2 k . Moreover, this average intensity E is calculated by using the grey level position of image I ′ . fourth moment of the mean is µ 4 . Here, the value of n 0 for contrast is fixed as 0.25. https://www.indjst.org/

• Directionality
The directionality provides information about global feature that is extracted from the small region and these features are represented based on the local edge's histogram as shown in equation (5) where, the quantized direction code is represented as ϕ , the histogram peak value is P, the Pth peak range among valleys is denoted as W P and the histogram's Pth peak position is represented as ϕ P .

Local ternary pattern
The LTP is generally an extension of LBP wherein the neighbor pixel value is encoded into 3 valued codes based on the user threshold. The obtained binary pattern from the LTP is utilized as the element of local image descriptor. Therefore, the LTP allocates the label to each pixel p i of an image (I ′ ) based on thresholding value. The central pixel p c and adjacent values are used to assign the label. Equation (6) and (7) expresses the LTP operator.
where, the user-defined threshold is represented asw.

HOG descriptor
HOG descriptor generally depends on the gradient directions of the pixel wherein the pixel's HOG descriptors are computed by considering a small spatial region that is specified as a cell. Here, the feature vector is taken from the concatenation of 1D histogram and this feature vector is used for further process. The preprocessed image I ′ is separated into cells as N × N pixels and these pixels are subjected to further processing by hog descriptor wherein the pixel's orientation of the gradient can be expressed as: Here, the orientations of θ j i , i = 1, . . ., N 2 related to the same cell j are quantized and accumulated in feature vectors. Moreover, the obtained histograms are arranged and concatenated in a unique HOG histogram which is considered as the final resultant of HOG feature vector.
The feature vector (F) is generated by using the features from the tamura, LTP and HOG descriptor. The following equation (9) represents the generated feature vector

Infinite feature selection
The prime advantage of Inf-FS approach is to highlight the importance of particular features including all possible feature subsets. The appropriate feature selection is used for constructing the low dimensional feature vectors to reduce the time consumption and to increase classification accurateness. The set of features given to the feature selection is F = { f (1) , f (2) , . . . , f (n) }. Ultimately, the graph G is fully connected, such that each vertex is related to the features and edges which specifies the model pairwise relations between the features. The G is represented by Z and the nature of edges is specified by each element z i j of Z that is as represented by equation (10).
Where, the loading coefficient is α ε[0, 1], standard deviation through the samples is σ (i) and the correlation coefficient is represented as a spearman coefficient. The pairwise energy shows that feature between f (i) and f ( j) is discriminative. https://www.indjst.org/ The finite path among the vertices i and j is represented as γ = {u 0 = i, u 1 , . . . , u l−1 , u l = j} that is the subset of feature pairs over the path and energy of γ is given by equation (11) The i th feature energy is computed by maximizing the path length to infinity. The i th feature energy is expressed in the following equation (12).
Where, the vertices set is represented as U; path length among i and j is Q l i, j . The divergence is happening due to the sum of infinite Z l . The real-valued regularization factor of equation (13) is used to assure the infinite sum convergence.
The v R (i) is calculated effectively by utilizing the convergence property of geometric power series of the matrix that is given in the equation (14).
The quantity marginalization is used to obtain the final energy score for each feature which is expressed in the equation (15).
The selected feature ranks are obtained by the v R(i) and energy values are arranged in the descending order. The selected features are represented as Y . This Inf-FS is utilized to select only noteworthy features from the feature set i.e., 2585 features are selected from the 8615 features and these selected features are given to the DNN for improving the classification.

Sparse autoencoder based DNN for classification
The autoencoder is applied in each layer of the DNN (44) . The given features are used to compute each label's probability related to the residue. The input feature given to this encoder based DNN is Y = (y 1 , y 2 , . . . , y i , ..y M } , y i εS L , where the amount of features is L and M is the length. The input to this sparse autoencoder is selected features from the infinite feature selection. In this neural network, the sigmoid function is used as an activation function. The feature representation for the input matrix Y is shown in the following equation (16).
The conventional autoencoder detects the approximation to the identity function. In sparse auto encoder, the amount of active neurons is constrained by adding the sparse penalty term in the objective function. The equation (17) shows the cost function of the sparse autoencoder.
Where, the average sum of squares is represented in the first term and the number of examples in the training set is L. The relative weight of regularization term is controlled by using λ . The number of hidden neurons is s 2 , the sparsity penalty weight is β and the Kullback−Leibler divergence is KL (.) . The Kullback−Leibler divergence is expressed in the equation (18).
Therefore, this sparse auto encoder based DNN is helpful for classifying different classes of MRI images. https://www.indjst.org/

Euclidean distance based medical image retrieval
In image retrieval, the Euclidean distance measurement technique computes the similarity measure between the query image and retrieved image. This Euclidean distance based similarity measure is used due to its efficiency and effectiveness. The square root of the sum of the squared absolute differences are calculated to measure the distance among the two vectors of images. The matching between the query image and images from the database is expressed as: In this proposed system, an effective image retrieval and classification of image types are achieved by using an appropriate feature extraction and classification. Additionally, the Inf-FS is used to select optimum features among generated feature vectors. The DNN based classification is used in this medical image retrieval due to its robustness against huge databases.

Experimental setup and performance evaluation
The experimental setup and performance evaluation of the medical image retrieval system is described in this section. The proposed classification and retrieval method is applied over two different databases such as OASIS and CE-MRI. The simulation tool namely MATLAB R2019a used to implement the medical image retrieval with 4GB RAM as hardware requirements. In both the databases, 80% of the images are taken for training and 20% the images are used for testing the retrieval performances. The performance metrics used to analyze the proposed method is accuracy, sensitivity, specificity and error rate and it is explained as follows: Accuracy: Accuracy is defined as the percentage of correctly identified brain images from the OASIS and CE-MRI dataset. Sensitivity: Sensitivity is the amount of actual positive images that are correctly identified during classification. Specificity: Specificity is the amount of actual negative images that are correctly identified during classification. Error rate: The error rate is defined as the amount of error occurred during the image retrieval.

Average retrieval precision (ARP):
It is defined as the ratio between the average of the total number of relevant images retrieved and total retrieved images.

Performance evaluation of proposed method
The proposed method's performance over the OASIS and CE-MRI datasets are analyzed in this section. The performances are evaluated as classification and retrieval results over medical images. In classification, the MRI images are classified as seven classes such as meningioma, glioma, pituitary, group1, group2, group3 and group4. The images exist in each group are randomly selected based on its ventricular shape. The classification performance is improved by extracting the optimal features from the images using tamura, LTP and HOG descriptor. The classification results are assessed in terms of accuracy, sensitivity, specificity and E-rate. Classification results of the proposed method are validated with two different classifiers namely Random Forest Classifier (RFC) and Neural Network (NN). Moreover, the classification results are analyzed for individual features and hybrid features. For the hybrid features, the classification is performed by with and without feature selection. Also, based on the output of the classifier, TP, TN, FP, FN values are estimated and from these obtained values, various performance metrics like accuracy, sensitivity, specificity and error rate are calculated for both OASIS and CE-MRI dataset. Tables 1 and 2 , the accuracy of neural network and RFC for each of the single features is comparatively less. However, for hybrid features employing Infinite feature selection has better performance than without infinite feature selection method. Forexample, the accuracy of RFC for OASIS dataset is 61.39 % that is high when compared to the accuracy of NN i.e. 60.83%. https://www.indjst.org/ Similarly, the performance analysis for OASIS dataset is depicted in Table 3 . From Table 3, the proposed method has better classification accuracy of 95.34% with feature selection that is high when compared to the accuracy of NN and RFC that is 60.83% and 61.39% respectively. The combination of high level feature (LTP, HOG) and low level feature (Tamura) helps to improve the classification accuracy. The performance of the classifier with Inf-FS provides better performance than the classifier without Inf-FS, because the Inf-FS is used to avoid the irrelevant features obtained from the LTP, HOG and Tamura. Here, the Inf-FS selects the relevant features from the feature vector based on the ranking order. Finally, the performance metrics for each of the classifier w.r.t feature selection mechanism for hybrid features are graphically presented as depicted in Figure 4 . The performance analysis for the CE-MRI dataset is shown in Tables 4, 5 and 6 respectively. From Table 6 , the proposed technique has better classification accuracy than the RFC and NN. For example, the accuracy of the proposed method for CE-MRI dataset is 99.87% with feature selection that is high when compared to the accuracy of NN and RFC that are 85.79% and 96.67% respectively and the same has been graphically presented in Figure 5 . The relevant feature selection of Inf-FS and https://www.indjst.org/ Here the Euclidean distance based image retrieval is used to retrieve the medical images. The Euclidean distance is computed between the query image vector with trained image feature vector. The retrieval results of proposed method is calculated in terms of Average Retrieval Precision (ARP). Similar to classification, the retrieval results are carried out for two different datasets with seven different classes, namely meningioma, glioma, pituitary, group1, group2, group3 and group4.   The ARP for OASIS and CE-MRI datasets are shown in Tables 7 and 8 respectively. The ARP is mainly calculated to analyze the retrieval performances for OASIS and CE-MRI datasets. Also, Figures 5 and 6 shows that the proposed method retrieves the appropriate images from the OASIS and CE-MRI dataset respectively. The ARP for OASIS dataset varied for each class due to the similar sizes of images present in the group1, group2, group3 and group4. Similarly, the ARP for CE-MRI dataset also varied for each class due to similar tumor portions of meningioma, glioma and pituitary tumor.

Comparative analysis
The average retrieval precision of the proposed method is compared to the existing techniques to know the effectiveness of the proposed method over medical image retrieval. The existing methods used for the comparison are LMVCoP (33) and CBIR-CNN (35) . In LMVCoP (33) , two different feature extraction approaches are used such as LDP and GLCM. The LMVCoP is a combination of the LVCoP and LMCoP. The deep CNN based method is used for feature extraction and transfer learning, block-wise fine-tuning strategies are utilized to perform the medical image retrieval. The image retrieval of LMVCoP (33) and CBIR-CNN (35) are analyzed over the OASIS and CE-MRI datasets respectively. Table 9. ARP comparison of proposed method with existing method OASIS dataset CE-MRI dataset LMVCoP (33) Proposed method CBIR-CNN (35) Proposed method 71.58 88.25 96.13 98.33 The comparative performance of the retrieval results is given in Table 9 . Here, the comparison is taken in terms of ARP for analyzing the retrieval performance in OASIS and CE-MRI datasets. For example, the graphical illustration of the ARP performance for OASIS dataset is given in Figure 8 . From the results, the proposed method has improved performance than the LMVCoP (33) and CBIR-CNN (35) during the medical image retrieval. The ARP of the proposed method is obtained as 98.33 for CE-MRI dataset, due to the combination of low and high level features and selection of relevant features from the feature vector. The higher ARP in CE-MRI dataset results in high classification accuracy and precision. Both LMVCoP (33) and CBIR-CNN (35) failed to achieve higher ARP in large datasets. The reason behind the higher retrieval performances are DNN and Euclidean distance based image retrieval over the variety of images. Since, the classification is improved by using an appropriate feature extraction and Inf-FS based feature selection.

Conclusion
In this work, an effective medical image retrieval and classification is obtained by using the Euclidean distance based retrieval and sparse auto encoder based DNN respectively. The classification performance is enhanced by using Inf-FS over the feature vectors of tamura features, LTP and HOG. This Inf-FS is used to remove the irrelevant features from the feature vector that helps to improve the classification accuracy. In that, tamura features define the texture features of the image and HOG descriptor provides the orientation of the image. Moreover, the Euclidean based image retrieval gives better performance over medical images due to its efficiency. The classification performance of the proposed method is analyzed with two classifiers namely NN and RFC. The accuracy of the proposed method through the CE-MRI dataset is 99.87% that is high when compared to the NN and RFC. The proposed method with Inf-FS for CE-MRI dataset has higher accuracy of 99.87% when compared to the without Inf-FS i.e., 99.17%. Moreover, the ARP of the proposed method for OASIS dataset is 88.25% that is high when compared to the LMVCoP method. In the future, an ensemble feature selection method can be used to improve the classification performances of the brain tumor.


Facial Expression Recognition Under Partial Occlusion from Virtual Reality Headsets based on Transfer Learning

Facial expressions of emotion are a major channel in our daily communications, and it has been subject of intense research in recent years. To automatically infer facial expressions, convolutional neural network based approaches has become widely adopted due to their proven applicability to Facial Expression Recognition (FER) task.On the other hand Virtual Reality (VR) has gained popularity as an immersive multimedia platform, where FER can provide enriched media experiences. However, recognizing facial expression while wearing a head-mounted VR headset is a challenging task due to the upper half of the face being completely occluded. In this paper we attempt to overcome these issues and focus on facial expression recognition in presence of a severe occlusion where the user is wearing a head-mounted display in a VR setting. We propose a geometric model to simulate occlusion resulting from a Samsung Gear VR headset that can be applied to existing FER datasets. Then, we adopt a transfer learning approach, starting from two pretrained networks, namely VGG and ResNet. We further fine-tune the networks on FER+ and RAF-DB datasets. Experimental results show that our approach achieves comparable results to existing methods while training on three modified benchmark datasets that adhere to realistic occlusion resulting from wearing a commodity VR headset. Code for this paper is available at: https://github.com/bita-github/MRP-FER

I. INTRODUCTION
Facial expressions of emotion are a major channel in daily communications to transmit and enhance information not provided by speech [1]. In recent years, machine-based analysis of expressions from human faces has attracted increasing attention as it has broad applications in various areas such as healthcare, emotionally sensitive robots, driver fatigue monitoring, interactive game design and social marketing [2], [3], [4]. However, automatic Facial Expression Recognition (FER) in an unconstrained situation is still difficult. One of the major obstacles for accurate FER outside laboratory environments is partial occlusion in the face [5]. It is very likely that some parts of the face become obstructed by sunglasses, a hat, hands, hair, and the like. In addition, recent advances in technology, specifically in Virtual Reality (VR) has shifted the way we communicate and interact with each other and the environment. However, in a virtual reality setting a large part of the face is occluded by a head-mounted display which blocks facial expressions, there by restricting engagement [6]. Thus, for VR systems to provide rich social interaction, it is vital to be able to recognize and represent these expressions.
Occlusion can substantially change the visual appearance of the face and severely decrease the performance of FER systems. The presence of occlusion increases the difficulty of extracting discriminative features from occluded facial parts due to inaccurate feature location, imprecise face alignment, or face registration error [7]. Types of partial occlusion can be divided into three main categories, temporary, systematic and hybrid [8], [9]. Temporary occlusion occurs when part of the face being obscured temporarily by self-occlusion (e.g., hand or head movement), other objects and change in environmental condition (e.g., illumination and lighting). Systematic occlusion can be caused by the existence of one or more facial components (e.g., hair, mustache, or a scar), or by people using different accessories (e.g. VR headset, glasses, hat, or mask). These types of occlusions are potentially more damaging since they result in whole features of relevance to judging facial expression being obscured. Hybrid occlusion arise in the presence of both systematic and temporary occlusion at the same time.
While lots of studies have been conducted on automatically inferring human expressions from images and video sequences, most of them focused on non-occluded situation and to our knowledge, few of them focused on user's expression in virtual reality environments. Moreover, although the state-of-the-art techniques in FER systems are highly effective for controlled laboratory environments, the existing approaches do not achieve the same accuracy for applications like VR systems in which severe occlusions exist. Furthermore, recently deep learning methods such as Convolutional Neural Networks (CNNs) have outperformed statistical methods. FER performance improvements have been provided by these methods. However training a deep architecture from scratch requires a lot of training data to ensure proper feature learning, and has the difficulty Figure 1. Example of placing a VR patch on a detected face. of adjusting many system parameters. Moreover, it needs expensive computational capacity. Existing CNN models are less accurate when handling severe systematic occlusion like VR setting where features of the upper half of the face are completely invisible.
In this paper we mainly focus on facial expression recognition in presence of a severe occlusion where the user is wearing a head-mounted display in a VR setting. A unique aspect of the systematic occlusion from arising from a VR headset is that it can be mathematically modeled. Since commodity VR headsets are of known size and shape, we can simulate the occlusion arising from these headsets, rather than trying to collect a new dataset of people wearing such headsets, thus saving significant man hours. In this way, we can employ transfer learning to utilize pre-trained FER networks while simulating occlusion from VR headsets.
The rest of the paper is organized as follows section II reviews related works, section III presents the proposed approach, including our geomtetric model for simulating occlusion and transfer learning-based networks. Section IV describes the experiments. Finally, conclusions are drawn in section V.

II. RELATED WORKS
Since partial occlusions are a common problem in facial recognition systems and also very frequent in real applications, we review the related approaches considering the similar task to ours which is facial analysis with occluded faces.
Some approaches attempt to minimize occlusion impact by reconstructing occluded parts. Cornejo et al. [10] proposed a methodology which is robust to occlusions through the Weber Local Descriptor (WLD). They used RPCA technique for reconstructing occluded facial expressions over the training set and projected all testing images into the space created by RPCA. The WLD descriptor is applied over the entire facial expression image for extracting textural features. For each feature vector they applied feature dimensionality reduction techniques, such as PCA and LDA sequentially and used SVM and K-NN for classification. Wang et al. [11] proposed a novel framework for FER under occlusion by fusing the global and local features. In global aspect, information entropy is employed to locate the occluded region and principal Component Analysis (PCA) method is adopted to reconstruct the occlusion region of image. Then, the occluded region replaced with the corresponding region of the best matched image in training set and Pyramid Weber Local Descriptor (PWLD) feature is extracted. At last, the outputs of SVM are fitted to the probabilities of the target class by using sigmoid function. For the local aspect, an overlapping block-based method is adopted to extract WLD features, and each block is weighted adaptively by information entropy, Chi-square distance and similar block summation methods are then applied to obtain the probabilities which emotion belongs to. Finally, fusion at the decision level is employed for the data fusion of the global and local features based on Dempster-Shafer theory of evidence.
Deep Convolutional Neural Networks have been pushing the frontier of face recognition over past years. To address the occlusion issue, Li et al. [12] proposed end-to-end trainable Patch-Gated Convolution Neutral Network (PG-CNN), a CNN with attention mechanism which tries to focus on different regions of the facial image and weighs each region according to its unobstructed-ness. It decomposes the feature maps of the whole face to multiple sub-feature maps to obtain diverse local patches and encoded them as a weighed vector by a patch-gated unit using attention net considering its unobstructed-ness. Both local and global representations of facial features are concatenated to serve as a representation of the occluded face, which would alleviate the impact caused by the lack of local expression information. Another approach to improve the ability of classification in occluded condition is transferring features. Xu et al. [13] developed a hybrid model based on deep convolutional networks to increase the robustness of transfer features from deep models. They built a deep CNN composed of four convolutional layers with max-pooling to extract features, followed by a multi-class SVM for emotion classification. To increase robustness to occlusion further, they improve model by merging high-level features of two trained deep ConvNets with the same structure, one is trained on non-occluded images, and the other is trained on same database with additive occluded samples. Non-occluded facial expression classifier can be used as a guidance to facilitate the process of an occluded facial expression classifier. Pan et al. [14] built two deep convolutional neural networks with the same architecture, these two networks are first pretrained with the supervised multi-class cross entropy losses. After pretraining, parameters of non-occluded network are fixed, and the occluded network is further fine-tuned under the guidance of non-occluded network. To incorporate guidance into the network during training, they introduce similarity constraint and loss. Also a mask learning strategy can be adapted to find and discard corrupted feature elements. Song et al. [15] designed a Pairwise Differential Siamese Network (PDSN) to establish a mask dictionary by exploiting the differences between the top convolutional features of occluded and nonoccluded face pairs. Each item of this dictionary captures the correspondence between occluded facial areas and corrupted feature elements, which is named Feature Discarding Mask (FDM). When dealing with a face image with random partial occlusions, they generate its FDM by combining relevant dictionary items and multiply it with the original features to eliminate those corrupted feature elements from recognition.The framework proposed in [22] utilizes an efficient filtering method to reduce the search space of face retrieval to increase scalability while keeping representative occluded faces within the search space.
To our knowledge just two studies focused on FER in VR setting. Georgescu and Ionescu [16] proposed a model for recognizing the facial expression of a person wearing a virtual reality (VR) headset which essentially occludes the upper part of the face. They trained VGG-f and VGG-face models, on modified training images in which the upper half of the face is completely occluded. This forces the neural networks to find discriminative clues in the lower half of the face. They proceed by fine-tuning the networks in two stages. In the first stage they fine-tune the CNN models on images with full faces. In the second they further fine-tune the models on images in which the upper half of the face is occluded. Hickson et al. [17] also considered facial expression recognition in VR setting. They proposed an algorithm to infer expressions by analyzing only the images of a person's eyes captured from an eye-tracking enabled VR head-mounted display through an IR gazetracking camera. In terms of taking user's face images, as not all VR headsets include eye-tracking sensors, we adopt similar approach to [16], considering external camera in VR system which is more generalizable and cost-effective. However, a key limitation in the approach in [16] is that they simply covered the upper half of the face without taking into consideration realistic occlusion arising from VR headsets, which we attempt to achieve in this work.

A. Simulation of partial occlusion
In a virtual reality setting, the upper half of the face is occluded by a VR headset. Since, there are no public stan-dard occlusion face images database containing systematic occlusion or persons wearing VR headset, we established VR-occluded images by masking upper region on the standard facial expression images. To simulate the occlusion, face detection is applied on gray-scale images based on a modification to the standard Histogram of Oriented Gradients and Linear SVM-based method for object detection. In order to estimate 68 coordinates of facial landmarks that map to facial structure, we followed the approach described in [18] trained on the iBUG 300-W face landmark dataset [19]. For applying VR patch, we initialize VR dimensions to 207.1 × 98.6mm based on the Samsung Gear VR headset. Typically, a resized set of training images have faces filling various percentages of the image, and it is unlikely that a fixed-sized covering mask fits all the training images. To uniformly scale the VR patch on the training images, we use the distance between the two temporal bones of the facial landmarks as a reference length. Then, we generate the polygonal occluding patch by setting the midpoint of the line passing through eye centre points as the centre coordinate of the VR headset. Moreover, to account for face rotations, we align the resized patch with the axis running through eye centres to reflect accurate simulation of a VR headset on various situations. We obtain the angle of incline by determining the inverse tangent function of changes in ycoordinates to changes in x-coordinates of eye centre points, and then we utilize the rotation matrix to rotate corner points of the blocking patch about its central pivot point on the coordinate plane accordingly. Figure 1 demonstrates the occlusion simulation process. This geometric model provides a more realistic occlusion resulting from wearing a VR headset, rather than simply covering the upper half of the face, as done in [16].

B. Transfer Learning
In CNNs, images from different datasets share similar low-level features after the convolution process. As it is costly to train a network from scratch, especially on a small-scale dataset, an alternative training strategy for a new dataset is utilizing parameters transferred from pretrained models and fine-tuning them on the basis of the new  [20] is a convolutional neural network, trained on the VGG-face dataset. It comprises of 13 convolutional layers, each containing a linear operator followed by one or more non-linearities such as ReLU and max pooling. Last three blocks are fully connected followed by a softmax layer of 2622 units for classification.

2) ResNet:
The residual neural network [21] implements the identity shortcut connection concept or residual block that skips one or more layers. ResNet won the first place on the ILSVRC 2015 classification task. ResNet follows VGG's convolutional layer design. The residual block has two 3 × 3 convolutional layers with the same number of output channels. Each convolutional layer is followed by a batch normalization layer and a ReLU activation function. Then, it skips these two convolution operations and add the input directly before the final ReLU activation function. Generally, the residual mapping is often easier to optimize in practice. we proceed our work by fine-tuning the models on occluded images generated through applying our geometric simulation model as described before. The custom VGG-face has 13 convolutional and 4 max-pooling layers as VGG-face, interleaved with dropout layers. More specifically, starting from second block, after each convolutional layer, a dropout layer is added with rate starting from 0.1 to 0.55. The dropout layers are effective in avoiding the model to over fit. After all the convolution layers, 2 dense layers are added, each with 512 hidden nodes, followed by a 50% dropout layer. The final dense layer is followed with a soft-max layer consisting of 8 nodes to generate the output. In the ResNet-50 architecture, layers pre-trained on VGGFace dataset are transferred to our deep CNN model, followed by 2 dense layers with 512 hidden nodes and replacing the last 1000 fully connected softmax layer by a 8 fully connected softmax layer. Figure 2 demonstrates the customized architechture of VGG-face.
To evaluate our models, we use three popular facial expression datasets, FER+, RAF-DB, and AffectNet. These datasets cover different scales of face images and contain images with temporary occlusions.
1) FER+: The FER+ dataset [23] is derived from FER2013 which is a large-scale and real-world dataset. It contains 8 emotional classes including anger, disgust, fear, happiness, neutral, sadness, surprise and contempt. Comparing to FER2013, some of the original images are relabeled, while other images, e.g. not containing faces, are completely removed. All face images in the dataset are aligned and resized to 48 × 48. The dataset consists of 28,709 training images, 3,589 validation images and 3,589 test images.
2) RAF-DB: Real-world Affective Faces Database [24] is an in-the-wild dataset which contains 30000 facial images posted in social networks and annotated by 40 annotators. The dataset includes 7 classes of basic emotions including surprised, fearful, disgusted, happy, sad, angry and neutral and 12 classes of compound emotions. In our experiment, only images with basic emotions were used, including 12,271 images as training data and 3,068 images as test data.
3) AffectNet: The AffectNet dataset [25] is one of the largest annotated datasets which contains more than one million images. It has 450000 manually annotated images with 8 basic emotional classes as FER+ and 3 other categories related to the intensity of valence and arousal including none, uncertain and non-face.

B. Training
During training, the input to our CNNs is set to a fixed sized 224 × 224 image. We rescaled all training, validation and test images to 224 × 224 pixels. All images are normalized by min-max normalization and the models are trained using data augmentation, which is based on including horizontally flipped images. The training is carried out by optimising the multinomial logistic regression objective using mini-batch gradient descent based on back-propagation with momentum. The batch size is set to 64 and momentum is set to 0.9. The training was regularised by weight decay, the L2 penalty multiplier set to 5 × 10 −4 and the ratio for dropout regularisation was set to 0.5. The learning rate was initially set to 10 −2 , and then decreased by a factor of 10 when the validation set accuracy stopped improving. In case of the ResNet50 architecture we also regularised training by setting max-norm kernel constraint.

C. Results
For each model, we train our network on each dataset separately with 5-fold cross-validation and report the accuracy numbers in Table I training models from scratch provides less accurate results in comparison with transfer-learning based approaches. The VGG-Face architecture pre-trained on VGG-Face dataset and fine-tuned on FER+ achieved accuracy of 79.98% which is the highest among the tested models in this paper. comparing the results with the results reported in [16], we see that we provide slightly better accuracy on AffectNet and slightly worse on FER+. However, it is important to note that our geometric occlusion simulation model adheres to more realistic occlusion resulting from a commodity VR headset (Samsung Gear VR), while in [16] the occlusion was generated by simply covering the upper half of the face. So these results are not directly comparable, however it shows that a simulated occlusion with transfer learning can provide promising results by utilizing existing benchmark datasets.

V. CONCLUSION
In this paper, we propose a transfer learning-based method for addressing the facial expression recognition in presence of severe occlusion where the user is wearing a headmounted display in a VR setting. We modified three existing benchmark datasets to mimic VR occlusion by utilizing a geometric simulation model based off of the Samsung Gear VR headset. We test two popular architectures, VGG and ResNet pretrained on the VGG-Face database and further fine-tune parameters of the models on the occluded images. We evaluate our method on FER+, AffectNet and RAF-DB datasets. Our method achieves performance comparable to the state-of-the art results. For future work, we will study how to integrate approaches using external cameras and methods utilizing internal cameras in VR system to improve the accuracy further and provide a rich engagement experience for virtual reality users.