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.