Sparse State-Space Realizations of Linear Controllers

Sparse State-Space Realizations of Linear Contr ollers Y aozhi Du, Jing Shuang (Lisa) Li Abstract — This paper provides a no vel appr oach for finding sparse state-space realizations of linear systems (e.g., con- trollers). Sparse controllers are commonly used in distributed control, where a controller is synthesized with some sparsity penalty . Her e, motiv ated by a modeling problem in sensorimo- tor neuroscience, we study a complementary question: given a linear time-inv ariant system (e.g., controller) in transfer function form and a desired sparsity pattern, can we find a suitably sparse state-space realization for the transfer function? This problem is highly nonconvex, but we propose an exact method to solve it. W e show that the problem reduces to finding an appropriate similarity transform from the modal realization, which in turn reduces to solving a system of multiv ariate polynomial equations. Finally , we leverage tools from algebraic geometry (namely , the Gr ¨ obner basis) to solve this problem exactly . W e provide algorithms to find real- and complex-valued sparse realizations and demonstrate their efficacy on several examples. I . I N T R O D U C T I O N State-space realization is the process of constructing a state-space model for a gi ven linear time-in variant (L TI) system in transfer function form. A sparse state-space real- ization consists of sparse system matrices. Sparse realizations are advantageous in engineering. T o regulate dynamics of large systems such as a power grid, applying centralized controllers leads to high computation and communication burden; this can be mitigated by distributed (i.e., sparse) controllers [1], [2]. A less con ventional application of sparse realizations is in sensorimotor modeling. Control theory is frequently applied to explain the behavior of animals; that is, the input-output sensorimotor behavior of animals can be modeled using a controller such as a transfer function [3]. It follows that such a controller is implemented by the nervous system. Since networks of neurons can be thought of as dynamical systems, the problem of translating behavior (i.e., transfer function) to neural implementation (i.e., network dynamics) can also be thought of as a realization problem. Furthermore, since biological neural networks are sparsely connected, [4], this becomes a sparse realization problem. In this paper , motiv ated by the biological angle presented abov e, we aim to solve the following problem: given input- output dynamics of an L TI system in transfer function form, and desired sparsity patterns, find a state-space realization for this system that adheres to the desired sparsity patterns. Our work focuses strictly on the realization of the transfer function without inv estigating the design of this transfer function. This is different from existing works in distributed The authors are with the Department of Electrical Engineering and Computer Science, Univ ersity of Michigan. { yzdu, jslisali } @umich.edu control, which solve for transfer function and sparse real- ization simultaneously using conv ex optimization [1], [5], [6] or iterative approaches [2], [7]. Few works focus on the sparse realization problem in isolation; [8] is one such example but it is inexact and unable to handle arbitrary sparsity constraints due to parameterization limitations. Our approach is as follows. First, we introduce relev ant definitions (Section II). W e then demonstrate that solving the sparse realization problem is equiv alent to solving a system of multiv ariate polynomial equations (Section III), and that tools from algebraic geometry [9]–[11] can be adapted to solve this problem (Section IV). W e conclude with simulations demonstrating the efficacy of our approach (Section V) and directions of future work (Section VI). Notation: δ ij is an n ˆ n matrix with the p i, j q -th entry equal to 1 and zeros ev erywhere else. k r x 1 , ¨ ¨ ¨ , x k s is the set of multiv ariate polynomials consisting of variables x 1 , ¨ ¨ ¨ , x k with coef ficients in k “ R or C . I I . P R O B L E M S E T U P Consider a linear time-in v ariant (L TI) system represented as a rational proper (i.e., causal) transfer function K p s q . Consider the follo wing state-space system: 9 z p t q “ Az p t q ` B y p t q u p t q “ C z p t q ` D y p t q (1) System Σ “ p A, B , C , D q is a realization of K p s q when K p s q “ C p sI ´ A q ´ 1 B ` D (2) Note that we deliberately choose our notation to suggest that K p s q is a controller with input y p t q (e.g., the output of some L TI plant) and output u p t q (e.g., control input), which reflects our motiv ating problem. Ho we ver , all results in this paper naturally carry over to realizations of generic L TI systems, not just controllers. Furthermore, while (1) is written in continuous-time, the realization problem is equiv alent in discrete-time; thus, all results from this paper are directly applicable to discrete-time systems as well. W e are interested in sparse realizations (i.e., A, B , C , D are sparse). W e first define some relev ant concepts. Definition 1. A minimal realization Σ “ p A, B , C , D q , A P C n ˆ n of K p s q is a real ization with the smallest possible n , In the rest of the paper, we only consider minimal real- izations. Definition 2. F or LTI system Σ “ p A, B , C , D q with A P C n ˆ n , B P C n ˆ m , C P C q ˆ n , and D P C q ˆ m , its support is S “ p A , B , C , D q wher e A P t 0 , 1 u n ˆ n , B P t 0 , 1 u n ˆ m , C P t 0 , 1 u q ˆ n , D P t 0 , 1 u q ˆ m and: A ij “ # 1 , A ij ‰ 0 0 , A ij “ 0 B ij “ # 1 , B ij ‰ 0 0 , B ij “ 0 C ij “ # 1 , C ij ‰ 0 0 , C ij “ 0 D ij “ # 1 , D ij ‰ 0 0 , D ij “ 0 (3) Definition 3. A state-space r ealization Σ “ p A, B , C, D q of K p s q satisfies a support S if: @ i, j s.t. A ij “ 0 , A ij “ 0 @ i, j s.t. B ij “ 0 , B ij “ 0 @ i, j s.t. C ij “ 0 , C ij “ 0 @ i, j s.t. D ij “ 0 , D ij “ 0 (4) Definition 4. Rational pr oper transfer function K p s q is realizable on a support S “ p A , B , C , D q if ther e exists a r ealization Σ of K p s q s.t. Σ satisfies S . W e are interested in the following problems: Problem 1. Given a rational pr oper transfer function K p s q and support S “ p A , B , C , D q with dimensions that match a minimal r ealization, find a r ealization Σ “ p A, B , C , D q of K p s q s.t. Σ satisfies S . Problem 2. Given a rational pr oper transfer function K p s q and support S “ p A , B , C , D q with dimensions that match a minimal r ealization, find a r ealization Σ “ p A, B , C , D q of K p s q s.t. A , B , C , D are real-valued matrices and Σ satisfies S . T ri vially , solutions to Problem 2 are solutions to Problem 1; howe ver , finding solutions to Problem 2 are more in volv ed (see Section IV). In the remainder of the paper , we will first propose an approach to solve Problem 1, then introduce additional tools to solve Problem 2. I I I . E Q U I V A L E N C E O F S PA R S E R E A L I Z A T I O N T O M U LTI V A R I A T E P O LY N O M I A L E Q U A T I O N S In this section, we will explain how Problems 1 and 2 can be equiv alently formulated as finding solutions to systems of multiv ariate polynomial equations. W e will deri ve: 1) A necessary and sufficient criterion that determines existence of solutions of Problem 1. 2) A sufficient criterion that guarantees solutions of Prob- lem 2 exists. First, we recall two well-known results on realizations. Theorem 1 ( [12, Prop 2.3]) . F or two minimal r ealizations of K p s q written as Σ “ p A, B , C , D q and ˜ Σ “ p ˜ A, ˜ B , ˜ C , D q , ther e exists a unique in vertible matrix T P C n ˆ n s.t. ˜ A “ T ´ 1 AT , ˜ B “ T ´ 1 B , ˜ C “ C T (5) Note that matrix D is the same for all minimal realizations of K p s q ; its sparsity is fixed and cannot be altered through similarity transforms. This prompts the follo wing corollary: Corollary 1. Given rational pr oper transfer function K p s q and support S “ p A , B , C , D q , we have D “ lim s Ñ8 K p s q . Then, if D i, j s.t. D ij “ 0 and D ij ‰ 0 , K p s q is not r ealizable on support S . As shown above, it is trivial to determine whether sparse realizations exist when K p s q is non-strictly proper (i.e., D ‰ 0 ). For the rest of the paper , we assume that the con- ditions of Corollary 1 are not met, i.e., @ i, j s.t. D ij “ 0 , lim s Ñ8 K ij p s q “ 0 , so that K p s q may be realizable on support S . W e no w focus on A , B , C . W e can alw ays find a minimal realization of K p s q with diagonal (i.e., extremely sparse) A : Lemma 1. F or rational proper transfer function K p s q with input dimension m , output dimension p , ther e exists minimal modal realization Σ p “ p A p , B p , C p , D p q of K p s q of the form: A p “ » — – p 1 . . . p n fi ffi fl , B p “ » — – b 11 ¨ ¨ ¨ b 1 m . . . . . . b n 1 ¨ ¨ ¨ b nm fi ffi fl , C p “ » — – c 11 ¨ ¨ ¨ c 1 n . . . . . . c p 1 ¨ ¨ ¨ c pn fi ffi fl , D p “ » — – d 11 ¨ ¨ ¨ d 1 m . . . . . . d p 1 ¨ ¨ ¨ d pm fi ffi fl (6) A p P C n ˆ n , B p P C n ˆ m , C p P C q ˆ n , D p P C q ˆ m such that K p s q “ C p p sI ´ A p q ´ 1 B p ` D p . Mor eover , when p i P R ( i “ 1 , ¨ ¨ ¨ , n ) and coefficients in K p s q ar e real numbers, we can find modal realizations with r eal-valued A p , B p , C p , D p . For the remainder of the paper we limit ourselves to considering K p s q with simple (i.e., non-repeated) poles, which admit modal realizations. Modal realizations have sparse A and dense B and C . W e can use Theorem 1 and Lemma 1 to equiv alent pose Problems 1 and 2 as problems of determining the existence of similarity transform T . Corollary 2. A solution to Pr oblem 1 e xists if and only if for any modal r ealization Σ p “ p A p , B p , C p , D p q of K p s q , ther e exists an invertible T P C n ˆ n , such that Σ T “ p T ´ 1 A p T , T ´ 1 B p , C p T , D p q satisfies S . A solution to Pr oblem 2 exists if and only if for any modal r ealization Σ p “ p A p , B p , C p , D p q of K p s q , ther e exists an in vertible T P C n ˆ n , such that Σ T “ p T ´ 1 A p T , T ´ 1 B p , C p T , D p q satisfies S , and T ´ 1 A p T , T ´ 1 B p , C p T , D p ar e real-valued matrices. In vertible matrices T are products of inv ertible elementary matrices; this can be seen from row reduction. There are three types of elementary matrices. They correspond to row (or column) addition, multiplication of a row (or column) by a scalar , and permutation of rows (or columns). An elementary addition matrix is an identity matrix with one of f-diagonal entry replaced by a scalar x : E x i 0 j 0 “ I n ` xδ i 0 j 0 , p E x i 0 j 0 q ´ 1 “ I n ´ xδ i 0 j 0 (7) Left multiplying a matrix M by elementary addition matrix E x i 0 j 0 means x times the j 0 -th row is added to the i 0 -th ro w of M . An elementary permutation matrix E p i 0 j 0 is an identity matrix with rows i 0 and j 0 permuted. Left multiplying a matrix M by E p i 0 j 0 permutes the i 0 -th and j 0 -th rows of M . An elementary multiplication matrix E m i 0 is an identity matrix with a ro w i 0 multiplied by a scalar a . Left multiplying a matrix M by E m i 0 multiplies the i 0 -th row of M by scalar a . W e can apply a sequence of row operations (i.e., elemen- tary matrix multiplications) to con vert in vertible matrix T to the identity matrix. W e can first apply row permutation, then apply row addition and finally apply row multiplication. Each of those elementary operation are equiv alent to left multiplying an elementary matrix. W e denote the product of permutation matrices as T p , the product of elementary addi- tion matrices T a and the product of elementary multiplication matrices as T m . Then, we can in vert T as follows T m T a T p T “ I , T ´ 1 “ T m T a T p (8) The order of row additions is not unique, so T a is not unique. W e use the following method to construct a valid, fixed T a to be used for the rest of the paper; subsequent analysis holds for other constructions of T a as well. Giv en T p T , first cancel out the nonzero entries of rows 2 , ¨ ¨ ¨ , n in the first column by adding the appropriate scalar times the first row to each row . Then, cancel out the nonzero entries of rows 1 , 3 , ¨ ¨ ¨ , n in the second column by adding the appropriate scalar times the second row to each row . Then, cancel out the nonzero entries of rows 1 , 2 , 4 , ¨ ¨ ¨ , n in the third column, and so on. Repeat this process; the last step is to cancel out the nonzero entries of rows 1 , ¨ ¨ ¨ , n ´ 1 of the n -th column by adding a proper scalar times the n -th row to each row . For the rest of the paper , T a represents the operation described abov e. When T is in vertible, such a sequence of operations is always possible when T p is properly chosen. T p , like T a , is not unique; we use P p to denote the set of possible v alues of T p , with |P p | “ n ! . W e discuss the choice of T p in subsequent sections. In total, n p n ´ 1 q ro w additions were made by T a . By (7), each entry of an elementary addition matrix and its in verse are a polynomial of the matrix’ s parameter (i.e., the scalar v alue x ) with real coefficients. W e denote those parameters as x 1 , ¨ ¨ ¨ , x n p n ´ 1 q . Then, each entry of T a and T ´ 1 a are polynomials of x 1 , ¨ ¨ ¨ , x n p n ´ 1 q with real coeffi- cients, written as: r T a s ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s , “ T ´ 1 a ‰ ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s . These can be combined into a matrix of polynomials: T a “ T a p x 1 , ¨ ¨ ¨ , x n p n ´ 1 q q . By Corollary 2, a solution to Problems 1 and 2 exists if and only if we can find a T m , T p and T a . In fact, we do not need to solve for T m because T m does not change the sparsity: Theorem 2. Rational proper transfer function K p s q is r ealizable on a support S if and only if D invertible T p , T a s.t. for ˜ A “ p T a T p q A p p T a T p q ´ 1 , ˜ B “ p T a T p q B p , ˜ C “ C p p T a T p q ´ 1 , ˜ Σ “ p ˜ A, ˜ B , ˜ C , D p q satisfies S . Pr oof. By Corollary 2 and (8), K p s q is realizable on S if and only if we can find ˜ Σ a “ p T m ˜ AT ´ 1 m , T m ˜ B , ˜ C T ´ 1 m , D p q . Since T m is a diagonal and in vertible, ˜ Σ satisfies support S if and only if ˜ Σ a satisfies S . Thus, K p s q is realizable on S if and only if we can find ˜ Σ “ p ˜ A, ˜ B , ˜ C , D p q , which prov es the claim. Based on Theorem 2, we propose the following criterion to determine if the Problems 1 and 2 are feasible. Corollary 3. Given modal r ealization Σ p “ p A p , B p , C p , D p q and support S : 1) Pr oblem 1 is solvable if and only if D T p P P p and x i P C for T a “ T a p x 1 , ¨ ¨ ¨ , x n p n ´ 1 q q such that ˜ Σ “ p ˜ A, ˜ B , ˜ C , D p q satisfies S , wher e ˜ A “ p T a T p q A p p T a T p q ´ 1 , ˜ B “ p T a T p q B p , ˜ C “ C p p T a T p q ´ 1 . 2) Assume K p s q is rational with r eal-valued coefficients in its numerator and denominator , and A p , B p , C p , D p ar e r eal matrices. Then, Pr oblem 2 is solvable if D T p P P p and x i P R for T a “ T a p x 1 , ¨ ¨ ¨ , x n p n ´ 1 q q such that ˜ Σ “ p ˜ A, ˜ B , ˜ C , D p q satisfies S , wher e ˜ A “ p T a T p q A p p T a T p q ´ 1 , ˜ B “ p T a T p q B p , ˜ C “ C p p T a T p q ´ 1 . Pr oof. For the first claim, as T a “ T a p x 1 , ¨ ¨ ¨ , x n p n ´ 1 q q , existence of T a is equiv alent to existence of x i . Then by Theorem 2 the claim holds. For the second claim, we only need to show that ˜ A , ˜ B , ˜ C are real-v alued. When x i P R , as r T a s ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s , T a P R n ˆ n . As T p P R n ˆ n and Σ p contains real-valued matrices, ˜ A , ˜ B , ˜ C , are real-valued as well, which prov es the claim. Corollary 3 implies that we can treat Problems 1 and 2 as finding solutions to multiv ariate polynomial equa- tions. Since r T a s ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s and “ T ´ 1 a ‰ ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s , we hav e that each entry of ˜ A , ˜ B and ˜ C is a multiv ariate polynomial, expressed as ˜ A ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s , ˜ B ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s , ˜ C ij P R r x 1 , ¨ ¨ ¨ , x n p n ´ 1 q s . ˜ Σ satisfies S if and only if: ˜ A ij “ 0 @ i, j s.t. A ij “ 0 ˜ B ij “ 0 @ i, j s.t. B ij “ 0 ˜ B ij “ 0 @ i, j s.t. C ij “ 0 (9) W e can observe that (9) is a system of polynomial equations. Thus, we can analyze feasibility of Problems 1 and 2 by determining existence of solutions to (9). While Problems 1 and 2 are concerned with finding real- izations for specific v alues of K p s q , we notice that for certain supports (e.g., the support of the modal realization), we can find a solution to (9) and correspondingly , a realization for any K p s q with an appropriate order . This motiv ates the following definitions: Definition 5. A support S with A P t 0 , 1 u n ˆ n is generally realizing if for all K p s q with n or less poles, K p s q is r ealizable on S . It is obvious that when S is the support of a modal realization, controllable canonical realization or observable canonical realization, S is generally realizing. The following lemma states that if a support is generally realizing, any support that is “strictly denser” than this support is also generally realizing. Lemma 2. If S 0 “ p A 0 , B 0 , C 0 , D s 0 q is gener ally r ealizing and S “ p A , B , C , D s q satisfies: @ i, j r A 0 s ij ‰ 0 ù ñ A ij ‰ 0 @ i, j r B 0 s ij ‰ 0 ù ñ B ij ‰ 0 @ i, j r C 0 s ij ‰ 0 ù ñ C ij ‰ 0 @ i, j r D s 0 s ij ‰ 0 ù ñ D s ij ‰ 0 (10) Then S is generally r ealizing. W e now present a nov el generally realizing support for single-input single-output (SISO) systems: Lemma 3. A generally r ealizing support S “ p A , B , C , D q for SISO systems with n states is: A ij “ # 1 i “ 1 or j “ 1 or i “ j 0 otherwise B i 1 “ # 1 i “ 1 or i “ n 0 otherwise C ij “ # 1 j “ 1 or j “ n 0 otherwise , D “ “ 1 ‰ (11) Pr oof. As D “ 1 , we only need to consider the strictly proper part of K p s q . Any proper K p s q with simple poles can be expressed as: K p s q “ ř n i “ 1 c i s ´ p i ` d . A minimal modal realization is: A p “ » — – p 1 . . . p n fi ffi fl , B p “ » — – 1 . . . 1 fi ffi fl , C p “ “ c 1 ¨ ¨ ¨ c n ‰ , D p “ “ d ‰ (12) Define T 0 – ś 2 p n ´ 1 q i “ 1 E i with E i “ # I n ` x i δ n,i ` 1 i “ 1 , ¨ ¨ ¨ , n ´ 2 I n ` x i δ i ` 1 , 1 i “ n ´ 1 , ¨ ¨ ¨ , 2 p n ´ 2 q (13) Since the modal realization is a minimal realization and the system is SISO, c n ‰ 0 . Then, we pick: x i “ # c i ` 1 c n , i “ 1 , ¨ ¨ ¨ , n ´ 2 ´ 1 , i “ n ´ 1 , ¨ ¨ ¨ , 2 p n ´ 2 q (14) Consider p A c , B c , C c , d q where A c “ T 0 A p T ´ 1 0 , B c “ T 0 B p , C c “ C T ´ 1 0 ; this is a realization of K p s q which satisfies S . Therefore, any K p s q is realizable on S . For Problem 1, if the provided support S is generally realizing, then solutions to the problem tri vially exist. Deter- mining whether a gi ven support S is generally realizing is beyond the scope of this work; we next return to our original problems, which have now been shown to be equi v alent to solving systems of multi variate polynomial equations. I V . S O L V I N G M U LT I V A R I A T E P O LYN O M I A L E Q UA T I O N S W I T H G R ¨ O B N E R BA S I S In the previous section, we demonstrated that solving Problems 1 and 2 are equi v alent to finding solutions to multiv ariate polynomial equations. This has been studied in the field of algebraic geometry [9]–[11]. In this section, we lev erage these results to introduce algorithms that compute complex- and real-valued solutions to multiv ariate polyno- mial equations, which will then be used to solve Problems 1 and 2, respectively . For the duration of the section, the number of v ariables is denoted by k P N . As an example, for the polynomial equations generated by (9), k “ n p n ´ 1 q . A. Computing complex-valued solutions W e introduce theory on the existence of complex-v alued solutions to multi v ariate polynomial equations and propose Algorithm 1 that computes them. This algorithm is guaran- teed to solve Problem 1 if solutions exist. Giv en H – t h 1 , ¨ ¨ ¨ , h η u where h i P C r x x , ¨ ¨ ¨ , x k s . The ideal [9, Chapter 1.4] generated by h i written as I “ x h 1 , ¨ ¨ ¨ , h η y is the set: I “ t g P C r x 1 , ¨ ¨ ¨ , x k s | g “ η ÿ i “ 1 f i h i , f i P C r x 1 , ¨ ¨ ¨ , x k su (15) Ordering of monomials, which determines the relati ve order of two arbitrary monomials x α and x β , is crucial for polynomial operations. The algorithms applied in our paper rely on lexicographic order [9, Chapter 2.2]. A key property of such an order is that relati ve ordering x 1 ą x 2 ą ¨ ¨ ¨ x k holds. In a polynomial f “ ř α P U c α x α P C r x 1 , ¨ ¨ ¨ , x k s with U being the set of possible degrees of x , the leading term is c β x β where @ α P U zt β u , x β ą x α . The leading term is denoted by L T p f q . For a fixed monomial ordering, a Gr ¨ obner basis [9, Chapter 2.5] of an ideal I is a set of polynomials G “ t g 1 , ¨ ¨ ¨ , g ζ u such that the ideal generated by leading terms of g j , written as x L T p g 1 q , ¨ ¨ ¨ , L T p g ζ qy equals the ideal gen- erated by the set of leading terms for nonzero polynomials in I , written as x L T p I qy . When G is a Gr ¨ obner basis, x g 1 , ¨ ¨ ¨ , g ζ y “ I . A reduced Gr ¨ obner basis is a special type of Gr ¨ obner basis with minimal size (i.e., | G | ). It is unique for a giv en monomial ordering. Intuiti vely , we can think of the process of con verting a Gr ¨ obner basis to a reduced Gr ¨ obner basis as being similar to con verting a matrix to its reduced echelon form. Computing the Gr ¨ obner basis is desirable because we can use it to determine existence of solutions to multiv ariate polynomial equations and and compute them when the y exist. For an ideal I “ x h 1 , ¨ ¨ ¨ , h η y generated by polynomials h i , i “ 1 , ¨ ¨ ¨ , η and G “ t g 1 , ¨ ¨ ¨ , g ζ u being a Gr ¨ obner basis of I , solutions (possibly complex-v alued) to h i “ 0 , i “ 1 , ¨ ¨ ¨ , η and g j “ 0 , j “ 1 , ¨ ¨ ¨ , ζ are: V C p I q – t x P C k | f p x q “ 0 , @ f P I u (16) For the rest of the paper, we use V C p I q to denote the set of complex-v alued solutions. The existence of complex solutions is gi ven by the following claim [9, Chapter 4.1]: V C p I q “ H if and only if reduced Gr ¨ obner basis of I is t 1 u (17) There are multiple software tools to compute Gr ¨ obner bases. For example, in MA TLAB, the command gbasis returns the reduced Gr ¨ obner basis of a giv en ordering. If the reduced Gr ¨ obner basis is t 1 u , no solutions exists. Otherwise, at least one solution (possibly complex-v alued) exists. Giv en a Gr ¨ obner basis, we can find solutions using the elimination theorem [9, Chapter 3.1]. When G l “ t g 1 , ¨ ¨ ¨ , g ζ u is a Gr ¨ obner basis of and ideal I with respect to a term order satisfying @ τ ‰ l , x τ ą x k , l “ 1 , ¨ ¨ ¨ , k , the set G l l – G l X k r x l s (18) are polynomials that constraints x l . Any x l that solves g j “ 0 , j “ 1 , ¨ ¨ ¨ , ζ solves the polynomial equations generated by G l l . Leveraging this property , we can apply the following technique to find solutions. For l “ 1 , ¨ ¨ ¨ , k , generate G l l and use polynomials in G l l to find a solution candidate of ˜ x l . When G l l “ H , we can pick any ˜ x l P C . W e then plug ˜ x l into G and re-compute the reduced Gr ¨ obner basis. When the new reduced Gr ¨ obner basis is t 1 u , we need to find a new candidate of ˜ x l . Otherwise, we proceed. Properties of the reduced Gr ¨ obner basis guarantees that we can always find a v alid ˜ x l . This technique is summarized in Algorithm 1. Giv en a set of multiv ariate polynomial equations, it is guaranteed to return a solution when solutions exists. When solution does not exists, it returns “infeasible”. B. Computing r eal-valued solutions Next, we introduce theory on the existence of real-valued solutions. W e want to solve for h i p x q “ 0 , i “ 1 , ¨ ¨ ¨ , η , h i P R r x 1 , ¨ ¨ ¨ , x k s (19) The set of real-v alued solutions to (19) is: V R p I q – V C p I q X R k (20) In R r x 1 , ¨ ¨ ¨ , x k s , monomial x α – x α 1 1 x α 2 2 ¨ ¨ ¨ x α k k is of order α “ p α 1 , α 2 , ¨ ¨ ¨ , α k q P N k . For generic polynomials, we can impose some ordering and rewrite coefficients in a Algorithm 1 Computing complex-v alued solutions to a set of multi variate polynomial equations Input: H – t h 1 , ¨ ¨ ¨ , h η u Output: Solution x “ p x 1 , ¨ ¨ ¨ , x k q P C k or Infeasible . 1: Compute reduced Gr ¨ obner basis G of H 2: for j “ k to 1 do 3: if G “ t 1 u then 4: retur n Infeasible . 5: G j Ð G X C r x j s 6: if G j “ H then 7: Randomly pick ˜ x j P C . 8: else 9: Find ˜ x j that solves equations generated by G j 10: Substitute ˜ x j into G to yield G 1 11: Compute reduced Gr ¨ obner basis ˜ G of G 1 12: while ˜ G “ t 1 u do 13: Repeat Line 6-11 to generate another ˜ x j 14: x j Ð ˜ x j , G Ð ˜ G . 15: retur n x vector . As an example, for f p x q “ x 2 1 ` x 1 x 2 ´ 4 P R r x 1 , x 2 s : f p x q “ “ ´ 4 0 0 1 1 0 ¨ ¨ ¨ ‰ » — — — — — — — — — – 1 x 1 x 2 x 2 1 x 1 x 2 x 2 2 . . . fi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl – F X (21) where f p x q ’ s coefficient vector F is indexed as: F “ “ f 0 , 0 f 1 , 0 f 0 , 1 f 2 , 0 f 1 , 1 f 0 , 2 ¨ ¨ ¨ ‰ (22) with f α being the coef ficient for x α , α P N k . In this e xample, k “ 2 . For lar ger k , we can index F in a similar pattern. When f p x q is a finite-order polynomial, later terms of F are zero; we can truncate these without losing information. For example, in (21), α “ p α 1 , α 2 q , f α “ 0 for all α 1 ` α 2 ě 3 . T o find real solutions, we also need truncated moment matrices. A moment matrix M p γ q is an infinite matrix defined as: r M p γ qs α,β : “ γ α ` β , α, β P N k (23) A truncation M t p γ q of M p γ q is a submatrix of M p γ q that keeps r M p γ qs α,β with ř n i “ 1 α i ď t and ř n i “ 1 β i ď t . As an example, when n “ 2 : M 1 p γ q “ » – γ 0 , 0 γ 1 , 0 γ 0 , 1 γ 1 , 0 γ 2 , 0 γ 1 , 1 γ 0 , 1 γ 1 , 1 γ 0 , 2 fi fl (24) For polynomials of the form f p x q “ ř β f β x β , we define f γ by: p f γ q α “ ÿ β f β γ α ` β , α P N k (25) Then, moment matrix M p f γ q is pp f γ q α ` β q α,β P N k . Next, define v ariables: d j – r deg p h j q{ 2 s , d – max j “ 1 , ¨¨¨ ,η d j (26) Then, existence of real-v alued solutions to the polynomials can be determined using moment matrices in the follo wing semidefinite program (SDP): min γ 1 s.t. γ 0 “ 1 , M t p γ q ľ 0 , M t ´ d j p h j γ q “ 0 , t ě d (27) When this SDP is infeasible, i.e., D t ě d such that no suitable γ can be found, then V R p I q “ H [11]. When it is feasible, we can find M t p γ q which contains a truncated sequence of γ . Then, for 1 ď s ď t , let M s p γ q be submatrices of M t p γ q with α , β satisfying ř n i “ 1 α i ď s , ř n i “ 1 β i ď s . If D 1 ď s ď t s.t. rank M s p γ q “ rank M s ´ 1 p γ q , we can compute solution candidates for (19) and verify whether they are solutions. The method for computing solution candidates for a giv en s is presented directly in [11] and is quite notationally heavy , so we summarize the general idea here: first, locate a set of linearly independent columns of M s p γ q , then extract sub- matrices of M s p γ q based on indices of linearly independent columns; inv ert one of these submatrices and perform matrix multiplication. Solution candidates of each x i are eigen values of the resulting matrix products. If no solutions are found, we can iteratively increase t and compute more solution candidates. When the polynomial equations hav e a finite number of solutions or no solutions, i.e., | V R p I q | ă 8 , this iteration is guaranteed to terminate, either with a solution or with a report of infeasibility if no solutions exist. Howe ver , termination fails to occur if | V R p I q | “ 8 . T o handle this case, we stop iterations at t “ d ` r for some parameter r . When r is sufficiently large, it is more likely that terminating at t “ r implies | V R p I q | “ 8 ; howe ver , increasing r also increases the time it takes to obtain a solution. This process is summarized in Algorithm 2. Algorithm 2 Computing real-valued solutions to a system of multi variate polynomial equations with finite solutions Input: H – t h 1 , ¨ ¨ ¨ , h η u , r Output: Solution x “ p x 1 , ¨ ¨ ¨ , x k q P R k or Infeasible or Infinite Solutions Likely . 1: for t “ d to d ` r do 2: Run SDP (27) with polynomials from H 3: if SDP Infeasible then 4: retur n Infeasible 5: if D 1 ď s ď t s.t. rank M s p γ q “ rank M s ´ 1 p γ q then 6: Extract solution candidates. 7: if D solution candidate x solving h i “ 0 then 8: retur n x 9: retur n Infinite Solutions Likely When a system of multi v ariate polynomial equations has infinite solutions, then Algorithm 2 cannot retrie ve a solution. W e now work tow ards a solution for this. First, we introduce the follo wing notation: Definition 6. Given V k p I q , k “ R or C , S l k p I q is the set of possible values for x l , written as: S l k p I q – t x l | x “ p x 1 , ¨ ¨ ¨ , x l , ¨ ¨ ¨ , x k q , x P V k p I qu (28) W e also define these sets. G l l was defined in (18) κ p I q – t x l | | S l R p I q | “ 8u (29) ϵ p I q – t x l | G l l “ Hu (30) One can observe that | V R p I q | “ 8 ð ñ | κ p I q | ą 0 . In this case, we want to assign valid v alues to some variables, plug them into h i “ 0 to yield ˜ h i such that ˜ I “ x ˜ h 1 , ¨ ¨ ¨ , ˜ h η y satisfies κ p ˜ I q . In other words, we con vert the system of polynomial equations with infinite solutions to a new system of polynomial equations with finite solutions; then, we can apply Algorithm 2 on this new system. W e now explain how to do this. As V R p I q Ď V C p I q , S l R p I q Ď S l C p I q . x l P S l C p I q ù ñ x l P ϵ p I q . Thus, x P κ p I q ù ñ x P ϵ p I q . W e want to pick an unkno wn in κ p I q . T o do so, we select x l P ϵ p I q , assign a random value ˜ x l and plug it into G l , yielding G l 1 . W e compute the reduced Gr ¨ obner basis of G l 1 , denoted by ˜ G l . At this point, we hav e two options that trade off rigorousness and efficienc y . The more rigorous approach is to run Algorithm 2 with ˜ G l . If the algorithm returns “Infeasible”, we select a new value of ˜ x l and try again. If the algorithm keeps returning “Infeasible” after r attempts, we conclude that x l R κ p I q and start again, picking x τ P ϵ p I q , τ ‰ l . If Algorithm 2 returns “infinite solutions”, we need to assign values to more v ariables; we set G “ ˜ G l 1 and re-compute ϵ p I q . If at any point, Algorithm 2 returns a solution, it is a solution to h i “ 0 and we are done. The do wnside of this Algorithm is that resulting calls to the SDP (27) contain very large moment matrices. T o a void this, we propose an alternativ e approach. A less rigorous but more efficient approach is to assume | V C p I q | “ 8 if and only if | V R p I q | “ 8 . Under this assumption, we can assign random real values to variables and update the Gr ¨ obner basis until complex solutions to the equations of remaining variables is finite. Then, we run Algorithm 2 to find values of the remaining variables. If the algorithm returns “Infeasible”, assign different random values in the first step and try again. Algorithm 3 implements this process. W e note that although this assumption does not generally hold, in all simulations and tests we performed related to this paper , we nev er encountered a scenario in which the assumption was violated. Remark: While conducting numerical experiments, we observed that in all cases where complex solutions existed, real solutions also existed; in these cases, we can very likely extract a real solution using Algorithm 1 directly . Howe ver , to guarantee retriev al of a real solution, we need to run Algorithm 3. Algorithm 3 Computing real-valued solutions to a system of multi variate polynomial equations Input: H “ t h 1 , ¨ ¨ ¨ , h η u , r Output: Solution x “ p x 1 , ¨ ¨ ¨ , x k q P R k or Infeasible . 1: Compute reduced Gr ¨ obner basis G of H 2: if G “ t 1 u then 3: retur n Infeasible 4: a2Out = Algorithm 2( H , r ) 5: if a2Out != Infinite Solutions Likely then 6: retur n a2Out 7: while true do 8: ˆ G Ð G . 9: Compute reduced Gr ¨ obner basis G l l of G . 10: while Exists l s.t. G l l “ H do 11: Assign a random v alue ˜ x l to x l . 12: Substitute ˜ x l into G l which yields G l 1 13: Compute reduced Gr ¨ obner basis ˜ G l of G l 1 14: while ˜ G l “ t 1 u do 15: Repeat Lines 10-14 to generate another ˜ x l 16: G Ð ˜ G l . 17: a2Out2 = Algorithm 2( G , r ). 18: if a2Out2 = x then 19: retur n x 20: else 21: G Ð ˆ G C. F inding sparse r ealizations W ith algorithms to compute solutions to multiv ariate polynomial equations, we preset Algorithm 4 to solve our original problems. Set input v ariable complex to true to solve Problem 1 and false to solve Problem 2. Note that in the pre vious section (Section III), we only deriv ed a sufficient but not necessary criterion for the existence of solutions to Problem 2; this criterion was then translated into a system of multiv ariate polynomial equations. Thus, ev en if we are unable to retriev e real-v alued solutions, this does not definitiv ely mean that no solutions exist; in this case, the algorithm reports No Solutions F ound as opposed to Infeasible . Also, this algorithm nominally iterates over all members T p P P p ; howe ver , in practice, setting T p “ I p often suffices (this is the case for all simulations in the next section). W e can also truncate this iteration for improv ed efficienc y . V . N U M E R I C A L E X A M P L E S W e now apply Algorithm 4 to specific numerical instances of Problems 1 and 2. All code required to reproduce our results can be found at h ttps: //gi thub .com / Yaoz D u/ Sp ar se - L TI - s ys te m- r ea li za ti on - v ia - s ol ving- multivariate- polynomial- equations . Consider a randomly generated unstable plant 9 x “ F x ` Algorithm 4 Computing sparse realizations Input: K p s q , S “ p A , B , C , D q , Boolean variable complex Output: State-space realization p A, B , C , D q of K p s q satis- fying S (real-v alued if complex “ false) or Infeasible or No Solutions Found . 1: Construct modal realization p A p , B p , C p , D p q of K p s q . 2: if D p does not satisfy S then 3: retur n Infeasible . 4: Find P p . 5: while D T p P P p not used do 6: Construct p ˜ A, ˜ B , ˜ C , ˜ D q with new T p . 7: Extract polynomial equation constraints H by (9) from p ˜ A, ˜ B , ˜ C , ˜ D q based on p A , B , C , D q . 8: if complex “ false then 9: Run Algorithm 3. 10: else 11: Run Algorithm 1. 12: if Algorithm returned a solution x then 13: Substitute x into p ˜ A, ˜ B , ˜ C , ˜ D q . 14: retur n p ˜ A, ˜ B , ˜ C , ˜ D q . 15: if complex “ true then 16: retur n Infeasible 17: else 18: retur n No Solutions Found (a) (b) Fig. 1: Controller structures from numerical examples Gu, y “ H x , where x is the plant state and F “ » — — – 0 . 2 0 . 6 0 . 6 0 . 2 0 . 6 0 . 8 0 . 8 0 . 6 0 0 . 4 0 . 6 0 . 2 0 . 6 0 . 40 0 . 4 0 fi ffi ffi fl , G “ » — — – 0 . 8 0 . 2 0 . 4 0 . 6 fi ffi ffi fl H “ “ 0 . 4 0 . 8 0 . 4 0 . 6 ‰ (31) An output feedback controller for this system is: 9 z “ p F ` GK ´ LH q z ` Ly , u “ K z (32) W e obtain an LQG controller by setting Q “ 11 I 4 , R “ 1 , W “ I 4 , V “ 1 where W and V represents the cov ariance of process and sensor noise. Now , we con vert this controller into a transfer function K p s q via K p s q “ K p sI ´ F ´ GK ` LH q ´ 1 L , and attempt to realize it on v arious supports. W e first define the following support, which is depicted in Figure 1a. A s “ » — — – 1 1 0 1 1 1 0 0 0 0 1 1 0 0 1 1 fi ffi ffi fl , B s “ » — — – 1 0 0 1 fi ffi ffi fl C s “ “ 1 0 1 0 ‰ , D s “ 0 (33) W e now wish to find a real-valued realization of K p s q that respects these supports, i.e., solve Problem 2. As previously described, we start with a modal realization A p , B p , C p , D p for this controller then solve for similarity transform ma- trices. This leads to a system of 11 polynomial equations with 12 variables, which has infinite real-valued solutions. Algorithm 4 returns the follo wing realization: A “ » — — – ´ 7 . 10 ´ 5 . 44 0 165 . 23 ´ 27 . 28 ´ 2 . 41 0 0 0 0 1 . 66 ´ 10 . 38 0 0 0 . 36 ´ 2 . 221 fi ffi ffi fl , B “ » — — – ´ 7 . 30 0 0 0 . 17 fi ffi ffi fl C “ “ 23 . 01 0 11 . 34 0 ‰ , (34) Next, we consider a different support, which is depicted in Figure 1b. A s “ » — — – 1 0 0 1 1 1 1 0 1 0 1 0 0 0 0 1 fi ffi ffi fl , B s “ » — — – 1 0 0 1 fi ffi ffi fl C s “ “ 1 1 0 0 ‰ , D s “ 0 (35) W e now wish to find a (possibly) complex-valued realization of K p s q that respects these supports, i.e., solve Problem 1. This leads to a system of 12 polynomial equations with 12 variables. Algorithm 4 returns the following realization: A “ » — — – ´ 17 . 16 0 0 ´ 14 . 94 ´ 17 . 52 7 . 64 ´ 0 . 57 0 0 . 47 0 ´ 0 . 48 0 0 0 0 0 . 08 fi ffi ffi fl , B “ » — — – ´ 8 . 53 0 0 ´ 0 . 76 fi ffi ffi fl C “ “ 19 . 66 ´ 8 . 48 0 0 ‰ (36) Notice that, as pre viously remarked, this algorithm returned a real-v alued solution ev en though it is not guaranteed to. W e consider one final support, which is extremely similar to the first support but with one additional zero: A s “ » — — – 1 0 0 1 1 1 0 0 0 0 1 1 1 0 1 1 fi ffi ffi fl , B s “ » — — – 1 0 0 1 fi ffi ffi fl C s “ “ 1 0 1 0 ‰ (37) W e attempt to solve Problem 1 on this support by again running Algorithm 4. This leads to a system of 12 polynomial equations with 12 variables. This time, the algorithm returns “Infeasible”; no complex or real-valued realizations exist. It’ s interesting to note that many sparse realizations of this LQG controller exist when the original state-space realization of this controller in (32) is dense. V I . C O N C L U S I O N S A N D F U T U R E W O R K W e proposed a novel approach to realize rational proper transfer functions on sparse state space models, and validated our approaches via numerical examples. During the process of formulating and testing our algorithms, we also identified the follo wing directions of future in vestig ation: 1) Deriving necessity and/or suf ficiency results for gen- eral realizing supports. This is related to the notion of canonical microcircuits in neuroscience [3] 2) Improving scalability of Algorithm 1 and Algorithm 3 3) Extending Algorithm 3 to accommodate transfer func- tions with complex conjugate poles Finally , we observed throughout simulations that Algo- rithm 1 always returned real-valued solutions when they ex- isted despite being formulated for complex-v alued solutions. W e are interested in identifying whether this property can be proven; if so, this would offer substantial computational speedups (compared to Algorithm 3). R E F E R E N C E S [1] S ¸ . Sab ˘ au, A. Speril ˘ a, C. Oar ˘ a, and A. Jadbabaie, “Network realization functions for optimal distributed control, ” IEEE T ransactions on Automatic Contr ol , vol. 68, no. 12, pp. 8059–8066, 2023. [2] C. W enk and C. Knapp, “Parameter optimization in linear systems with arbitrarily constrained controller structure, ” IEEE T ransactions on Automatic Control , vol. 25, no. 3, pp. 496–500, 2003. [3] J. S. Li, “T oward neuronal implementations of delayed optimal con- trol, ” in 2025 American Control Conference (ACC) . IEEE, 2025, pp. 2715–2721. [4] J. H ¨ anggi, L. F ¨ ovenyi, F . Liem, M. Meyer , and L. J ¨ ancke, “The hypothesis of neuronal interconnectivity as a function of brain size—a general organization principle of the human connectome, ” Fr ontiers in human neuroscience , vol. 8, p. 915, 2014. [5] J. Anderson and N. Matni, “Structured state space realizations for sls distributed controllers, ” in 55th annual Allerton confer ence on communication, control, and computing . IEEE, 2017, pp. 982–987. [6] C. W . Scherer , “Design of structured controllers with applications, ” in Pr oceedings of the 39th IEEE Confer ence on Decision and Control (Cat. No. 00CH37187) , vol. 5. IEEE, 2000, pp. 5204–5209. [7] F . Lin, M. Fardad, and M. R. Jov anovi ´ c, “Design of optimal sparse feedback gains via the alternating direction method of multipliers, ” IEEE T ransactions on Automatic Control , vol. 58, no. 9, pp. 2426– 2431, 2013. [8] J. S. Li and D. Ho, “Separating controller design from closed-loop design: A new perspectiv e on system-lev el controller synthesis, ” in 2020 American Contr ol Conference (A CC) . IEEE, 2020, pp. 3529– 3534. [9] D. Cox, J. Little, D. O’ shea, and M. Sweedler , Ideals, varieties, and algorithms . Springer , 1997. [10] D. A. Cox, J. Little, and D. O’ shea, Using algebraic geometry . Springer , 1998. [11] J. B. Lasserre, M. Laurent, and P . Rostalski, “Semidefinite charac- terization and computation of zero-dimensional real radical ideals, ” F oundations of Computational Mathematics , vol. 8, no. 5, pp. 607– 647, 2008. [12] B. De Schutter , “Minimal state-space realization in linear system the- ory: an overvie w , ” Journal of computational and applied mathematics , vol. 121, no. 1-2, pp. 331–354, 2000.