Kronecker Products_KPpaper

The Kronecker Product and Stochastic Automata Networks Amy N. Langville∗and William J. Stewart† September 29, 2003 Abstract This paper can be thought of as a companion paper to Van Loan’s The Ubiquitous Kronecker Product paper [23]. We collect and catalog the most useful properties of the Kronecker product and present them in one place. We prove several new properties that we discovered in our search for a Stochastic Automata Network preconditioner. We conclude by describing one application of the Kronecker product, omitted from Van Loan’s list of applications, namely Stochastic Automata Networks. Key words: Stochastic automata networks, Kronecker products, Kronecker product properties, Precon- ditioning. ∗Operations Research Program, N. Carolina State University Raleigh, NC 27695-7913, USA [email protected] Phone: (919) 513-1907, Fax: (919) 513-1908 †Department of Computer Science, N. Carolina State University Raleigh, N.C. 27695-8206, USA [email protected] Phone: (919) 515-7824, Fax: (919) 515-7896 Research supported in part by NSF (CCR-9731856). 1 Introduction Stochastic Automata Networks (SANs) have become an increasingly important modeling tool since the 1980s. SANs are used to efficiently model very large Markov chains whose state space is on the order of millions. The key to a SAN’s ability to compactly and efficiently model such large Markov chains lies in their extensive use of the Kronecker product operation. In order to understand SANs and their advantages, one needs some familiarity with the Kronecker product. The first half of this paper (section 2) is meant to provide such familiarity by collecting many of the known names, definitions, and properties of the Kronecker product. In addition, three new properties pertaining to the Kronecker product’s compatibility with generalized inverses are proven in section 2.6. After this theoretical introduction to the Kronecker product, we describe the practical uses of the Kronecker product, listing several applications 1 of the operation, ranging from image processing and generalized spectral analysis to analysis of chess endgames and fast transform algorithms. The number of different uses of the Kronecker product has grown recently, prompting Charlie Van Loan to call the operation ubiquitous . His 2000 paper describes dozens of interesting applications [23]. We use the second half of this paper (section 4) to add one more application of the Kronecker product to Van Loan’s list: SANs. We use examples and a discussion of the solution methods for SANs to show the Kronecker product’s connection to SAN modeling. 2 The Kronecker product The operation defined by the symbol ⊗was first used by Johann Georg Zehfuss in 1858 [5]. It has since been called by various names, including the Zehfuss product, the Producttransformation, the conjunction, the tensor product, the direct product and the Kronecker product. In the end, the Kronecker product stuck as the name for the symbol and operation, ⊗. 2.1 Definition of the Kronecker Product Definition: The Kronecker product ofAmA×nA∈ /RfracturmA×nAandBmB×nB∈ /RfracturmA×nA, written A⊗B, is the tensor algebraic operation defined as A⊗B= a1,1B a 1,2B . . . a 1,nAB a2,1B a 2,2B . . . a 2,nAB ............ amA,1B a mA,2B . . . a mA,nAB . Eachai,jBis a block of size mB×nB.A⊗Bis of size mAmB×nAnB. For example, if A=/parenleftbigg a1,1a1,2a1,3 a2,1a2,2a2,3/parenrightbigg , B= b1,1b1,2 b2,1b2,2 b3,1b3,2 b4,1b4,2 , thenA⊗B= a1,1b1,1a1,1b1,2 a1,1b2,1a1,1b2,2 a1,1b3,1a1,1b3,2 a1,1b4,1a1,1b4,2a1,2b1,1a1,2b1,2 a1,2b2,1a1,2b2,2 a1,2b3,1a1,2b3,2 a1,2b4,1a1,2b4,2a1,3b1,1a1,3b1,2 a1,3b2,1a1,3b2,2 a1,3b3,1a1,3b3,2 a1,3b4,1a1,3b4,2 a2,1b1,1a2,1b1,2 a2,1b2,1a2,1b2,2 a2,1b3,1a2,1b3,2 a2,1b4,1a2,1b4,2a2,2b1,1a2,2b1,2 a2,2b2,1a2,2b2,2 a2,2b3,1a2,2b3,2 a2,2b4,1a2,2b4,2a2,3b1,1a2,3b1,2 a2,3b2,1a2,3b2,2 a2,3b3,1a2,3b3,2 a2,3b4,1a2,3b4,2 . We also mention another Kronecker operation, the Kronecker sum, which is defined as the ordinary sum of Kronecker products. The Kronecker sum, A⊕B, is defined by square matrices AandBand is 2 given by A⊕B/triangle=A⊗InB+InA⊗B, where nAis the size of the square matrix AandnBis the size of the square matrix B. One advantage of Kronecker products is their compact representation. Consider the linear system Cx=din which Ccan be written as the Kronecker product of two much smaller matrices, AandB. The system ( A⊗B)x=dcan be solved quickly without ever forming the full matrix C=A⊗B(as is shown in section 3); only the smaller matrices AandBneed to be stored. An iterative method such as GMRES that uses only matrix-vector multiplications can be used to solve the compact system ( A⊗B)x=dwith the Kronecker product-vector multiplication algorithm [3]. Suppose C10000 ×10000can be expressed as the Kronecker product of A100×100andB100×100. The linear system Cx=donly requires the storage of two 100×100 matrices. In fact, later we will exploit properties of the Kronecker product to solve the special system ( A⊗B)x=dvery fast. 2.2 Properties of the Kronecker Product Before we can discuss some of the interesting applications of the Kronecker product, a complete background of its properties is required. These properties are divided into categories by topic. For example, the first four properties listed are basic Kronecker product properties, while the next three deal with structure. 2.3 Basic Properties Graham’s book [4] lists the following properties (along with proofs) of the Kronecker product such as: 1. Associativity: A⊗(B⊗C) = (A⊗B)⊗C. 2. Distributivity over ordinary matrix addition: (A+B)⊗(C+D) =A⊗C+B⊗C+A⊗D+B⊗D. 3. Compatibility with ordinary matrix multiplication: AB⊗CD= (A⊗C)(B⊗D). 4. Compatibility with ordinary matrix inversion: (A⊗B)−1=A−1⊗B−1. 2.4 Structure and Factorization Properties Van Loan’s paper, [23], lists the following additional properties of the Kronecker product. 3 1. Compatibility with ordinary matrix transposition: (A⊗B)T=AT⊗BT. 2. Structure theorems: (a) If AandBare nonsingular, then A⊗Bis nonsingular. (b) If AandBare square lower (upper) triangular, then A⊗Bis lower (upper) triangular. (c) If AandBare banded, then A⊗Bis banded. (d) If AandBare symmetric, then A⊗Bis symmetric. (e) If AandBare positive definite, then A⊗Bis positive definite. (f) If AandBare stochastic, then A⊗Bis stochastic. (g) If AandBare Toeplitz, then A⊗Bis block Toeplitz. (h) If AandBare orthogonal, then A⊗Bis orthogonal. 3. Factorizations: (a)LU: Let Abe a square nonsingular matrix of order mAwithLUfactorization A=PT ALAUA andBbe a square nonsingular matrix of order mBwithLUfactorization B=PT BLBUB. Then A⊗B= (PT ALAUA)⊗(PT BLBUB) = (PA⊗PB)T(LA⊗LB)(UA⊗UB). (b) Cholesky: Let Abe a positive definite matrix of order mAwith Cholesky factor GAandBbe a positive definite matrix of order mBwith Cholesky factor GB. Then the Cholesky factorization ofA⊗Bis A⊗B= (GT AGA)⊗(GT BGB) = (GA⊗GB)T(GA⊗GB). (c)QR: Let Abe an mA×nAmatrix with linearly independent columns and QRfactorization A=QARA, where Qis anmA×nAmatrix with orthonormal columns and Ris ann×nupper triangular matrix. Bis similarly defined with B=QBRBas its QRfactorization. Then the QRfactorization of A×Bis A⊗B= (QARA)⊗(QBRB) = (QA⊗QB)(RA⊗RB). (d) Schur decomposition: Let Abe a square matrix of order mAwith Schur decomposition A= UATAUT A, where UAis unitary and TAis upper triangular. Let Bbe a square matrix of order mBwith Schur decomposition B=UBTBUT B, where UBis unitary and TBis upper triangular. Then the Schur decomposition of A⊗Bis A⊗B= (UATAUT A)⊗(UBTBUT B) = (UA⊗UB)(TA⊗TB)(UA⊗UB)T. (e) Singular value decomposition: Let Abe an mA×nAmatrix with singular value decomposition UAΣAVT AandBbe an mB×nBmatrix with singular value decomposition UBΣBVT B. Let rank(A) =rAandrank(B) =rB. Then A⊗Bhas rank rArBand singular value decomposition A⊗B= (UAΣAVT A)⊗(UBΣBVT B) = (UA⊗UB)(ΣA⊗ΣB)(VA⊗VB)T. NOTE: All of these factorizations of C=A⊗Bmerely require the factorizations of the small AandBmatrices! Most of the theorems have trivial proofs. Many of the proofs in this section can be found in [23], [4], [19], [10], or [6]. 4 2.5 Measure and Numerical Properties Chapter 4 of the book by Horn and Johnson, [6], contains a wealth of information on Kronecker products and their properties. Some of the more useful ones are listed below. 1. Trace: if AandBare square, then tr(A⊗B) =tr(A)tr(B) =tr(B⊗A). 2. Norms: If AismA×nAandBismB×nB, then for all p-norms /bardblA⊗B/bardbl=/bardblA/bardbl /bardblB/bardbl. 3. Rank: rank(A⊗B) =rank(A)rank(B). 4. Eigenvalues and Eigenvectors: ForAandBsquare, let λbe a member of the spectrum of A. That is, λ∈σ(A). Let xAbe a corresponding eigenvector of λand let µ∈σ(B) and xBbe a corresponding eigenvector. Then λµ∈σ(A⊗B) and xA⊗xBis the corresponding eigenvector of A⊗B. That is, every eigenvalue ofA⊗Barises as a product of eigenvalues of AandB. 5. Singular values: Let the rank(A) =rAandrank(B) =rB. Then the nonzero singular values of A⊗Bare the rArBpositive numbers {σi(A)σj(B) : 1 ≤i≤rA,1≤j≤rB}, where σi(A) is the ithsingular value of A. 6. Determinants: If Aism×mandBisn×nthen det(A⊗B) = [det(A)]n[det(B)]m. 7. Powers: If AandBare square then (A⊗B)n=An⊗Bn. Below we prove or provide references to the proofs of each of the seven theorems above. Proof of 1: LetAbem×mandBben×n. tr(A⊗B) =m/summationdisplay i=1tr(ai,iB) =m/summationdisplay i=1ai,itr(B) =tr(B)m/summationdisplay i=1ai,i=tr(B)tr(A). /square 5 Proof of 2: We begin by proving the Frobenius norm case, /bardblA⊗B/bardblF=/bardblA/bardblF/bardblB/bardblF. /bardblA⊗B/bardbl2 F=tr[(A⊗B)(A⊗B)T] =tr[(A⊗B)(AT⊗BT)] =tr(AAT⊗BBT) =tr(AAT)tr(BBT) =tr(ATA)tr(BTB) =/bardblA/bardbl2 F/bardblB/bardbl2 F= (/bardblA/bardblF/bardblB/bardblF)2. Therefore, /bardblA⊗B/bardblF=/bardblA/bardblF/bardblB/bardblF. Now for the 2-norm. /bardblA/bardbl2/bardblB/bardbl2=/radicalbig λmax(A)λmax(B) =/radicalbig λmax(A⊗B) =/bardblA⊗B/bardbl2. The 1-norm case, /bardblA⊗B/bardbl1= max 1≤jA≤nAmA/summationdisplay iA=1|aiAjAB| = max 1≤jA≤nA,1≤jB≤nBmA/summationdisplay iA=1mB/summationdisplay iB=1|aiAjAbiBjB| = max 1≤jA≤nAmA/summationdisplay iA=1|aiAjA|max 1≤jB≤nBmB/summationdisplay iB=1|biBjB| =/bardblA/bardbl1/bardblB/bardbl1. The∞-norm is similar to the 1-norm except the largest absolute row sum is used rather than the largest absolute column sum. /square Proof of 3: LetAbem×nandBbep×q. IfA=QARAandB=QBRBare the QR factorizations, where QAism×nandQBisp×q, then rank(A⊗B) = rank(QARA⊗QBRB) =rank((QA⊗QB)(RA⊗RB)) =rank(RA⊗RB). Since RAandRBare both upper triangular, then RA⊗RBis upper triangular with upper triangular blocks. Let rank(RA) =rAandrank(RB) =rB. Each row of blocks of size RBhasrB nonzero rows. There are rAnonzero rows of such blocks. Using this and the upper triangular structure of RA⊗RB, we conclude that rank(RA⊗RB) =rArB. Therefore, rank(A⊗B) =rank(RA⊗RB) =rArB=rank(RA)rank(RB) =rank(A)rank(B). /square Proof of 4 and 5: Statement of these theorems and their corresponding proofs can be found in Chapter 4 of the book by Horn and Johnson [6]. /square Proof of 6: LetAbem×mandBben×n. A determinant for an n×nmatrix Gcan be determined from an LU factorization with pivoting [11]. Then det(G) =σGuG1,1uG2,2· · ·uGn,n=σGdiagprod (UG), where σG= +1 if an even number of row interchanges are used to obtain PGand -1 if an odd number of row interchanges are used to obtain PG. With σAB=σAσB. det(A⊗B) =σABdiagprod (UA⊗UB) =σAB(uA1,1)ndiagprod (UB)(uA2,2)ndiagprod (UB)· · ·(uAm,m)ndiagprod (UB) =σAB(diagprod (UA))n(diagprod (UB))m = [det(A)]n[det(B)]m. 6 /square Proof of 7: Proof by induction on n. Base case ( n= 2): (A⊗B)2= (A⊗B)(A⊗B) =A2⊗B2. Induction Step: Assume ( A⊗B)n=An⊗Bn. Show ( A⊗B)n+1=An+1⊗Bn+1. (A⊗B)n+1= (A⊗B)n(A⊗B) = (An⊗Bn)(A⊗B) =An+1⊗Bn+1. /square 2.6 Pseudoinverse Properties We prove three new properties of the Kronecker product in our work with Markov chains and their preconditioners. To our knowledge, these properties have not been stated or proven elsewhere. Before we state our new theorems and proofs, we define some generalized inverses: the Drazin inverse, the group inverse and the Moore-Penrose pseudoinverse [2, 11, 12]. •IfAis ann×nsingular matrix of index ksuch that rank(Ak) =r, then there exists a nonsingular matrix Qsuch that Q−1AQ=/parenleftbigg Cr×r0 0N/parenrightbigg , where Cis nonsingular and Nis nilpotent of index k. The Drazin inverse ofA, denoted by AD, is given as AD=Q/parenleftbigg C−10 0 0/parenrightbigg Q−1. •The group inverse is a special case of the Drazin inverse and applies when the index of the matrix Ais 1. The group inverse is appropriate for singular n×nmatrices of rank n−1 and is denoted byA#. •IfAis anm×nmatrix of rank r, then there exist orthogonal matrices Um×mandVn×nsuch that A=URVT=U/parenleftbigg Cr×r0 0 0/parenrightbigg m×nVT. TheMoore-Penrose pseudoinverse ofA, denoted by A†is the n×mmatrix given as A†=V/parenleftbigg C−10 0 0/parenrightbigg n×mUT. With these definitions we are now ready to state our new theorems. 1. Condition number: For all matrix norms, cond(A⊗B) =cond(A)cond(B). 7 2. Compatibility with the Drazin inverse and the group inverse: (A1⊗A2⊗ · · · ⊗ An)D=AD 1⊗AD 2⊗ · · · ⊗ AD n. (A1⊗A2⊗ · · · ⊗ An)#=A# 1⊗A# 2⊗ · · · ⊗ A# n. 3. Compatibility with the Moore-Penrose pseudoinverse: (A1⊗A2⊗ · · · ⊗ An)†=A† 1⊗A† 2⊗ · · · ⊗ A† n. Proof of 1: Case 1—If AandBare nonsingular, then cond(A⊗B) = /bardblA⊗B/bardbl /bardbl(A⊗B)−1/bardbl =/bardblA⊗B/bardbl /bardblA−1⊗B−1/bardbl =/bardblA/bardbl /bardblB/bardbl /bardblA−1/bardbl /bardblB−1/bardbl =cond(A)cond(B). Case 2—If AandBare singular, then A⊗Bis singular and cond(A⊗B) = /bardblA⊗B/bardbl /bardbl(A⊗B)†/bardbl =/bardblA⊗B/bardbl /bardblA†⊗B†/bardbl =/bardblA/bardbl /bardblB/bardbl /bardblA†/bardbl /bardblB†/bardbl =cond(A)cond(B). Case 3—If Ais nonsingular and Bis singular, then A⊗Bis singular and cond(A⊗B) = /bardblA⊗B/bardbl /bardbl(A⊗B)†/bardbl =/bardblA⊗B/bardbl /bardblA†⊗B†/bardbl =/bardblA/bardbl /bardblB/bardbl /bardblA−1/bardbl /bardblB†/bardbl =cond(A)cond(B), sinceA†=A−1forAnonsingular. /square Proof of 2: (Proof by Induction on n) We begin with a proof of the base case, (A1⊗A2)D=AD 1⊗AD 2. We derive a different expression for the right-hand side, then show that the left-hand side can also be written this way. Every square singular matrix A1of order nA1can be decomposed as A1=PA1/parenleftbigg CA10 0NA1/parenrightbigg P−1 A1, where CA1is nonsingular of size rA1×rA1,NA1is nilpotent of index kA1and rank( AkA1 1) =rA1. Similarly, a square matrix A2of order nA2can be written as A2=PA2/parenleftbigg CA20 0NA2/parenrightbigg P−1 A2, again where CA2is nonsingular of size rA2×rA2,NA2is nilpotent of index kA2and rank(AkA2 2) =rA2. According to the definition of the Drazin inverse, AD 1=PA1/parenleftbigg C−1 A10 0 0/parenrightbigg P−1 A1, 8 and likewise AD 2=PA2/parenleftbigg C−1 A20 0 0/parenrightbigg P−1 A2. Thus, AD 1⊗AD 2=/bracketleftbigg PA1/parenleftbigg C−1 A10 0 0/parenrightbigg P−1 A1/bracketrightbigg ⊗/bracketleftbigg PA2/parenleftbigg C−1 A20 0 0/parenrightbigg P−1 A2/bracketrightbigg = (PA1⊗PA2) C−1 A1⊗C−1 A20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (P−1 A1⊗P−1 A2) = (PA1⊗PA2)/parenleftbigg C−1 A1⊗C−1 A20 0 0/parenrightbigg (P−1 A1⊗P−1 A2). Now we show that ( A1⊗A2)Dis equal to the above expression. (A1⊗A2)D=/braceleftbigg/bracketleftbigg PA1/parenleftbigg CA10 0NA1/parenrightbigg P−1 A1/bracketrightbigg ⊗/bracketleftbigg PA2/parenleftbigg CA20 0NA2/parenrightbigg P−1 A2/bracketrightbigg/bracerightbiggD =/braceleftbigg (PA1⊗PA2)/bracketleftbigg/parenleftbigg CA10 0NA1/parenrightbigg ⊗/parenleftbigg CA20 0NA2/parenrightbigg/bracketrightbigg (P−1 A1⊗P−1 A2)/bracerightbiggD =  (PA1⊗PA2) CA1⊗CA2 0 0 0 0 CA1⊗NA2 0 0 0 0 NA1⊗CA2 0 0 0 0 NA1⊗NA2 (P−1 A1⊗P−1 A2)  D . LetNLbe the 3 ×3 principal submatrix of the middle matrix above. NLis nilpotent of index k=max{kA1, kA2}. That is, Nk L= CA1⊗NA2 0 0 0 NA1⊗CA2 0 0 0 NA1⊗NA2 k = 0. This follows from ( CA1⊗NA2)k=Ck A1⊗Nk A2=Ck A1⊗0 = 0, and similarly for the other two diagonal blocks. Using this fact, (A1⊗A2)D=/braceleftbigg (PA1⊗PA2)/parenleftbigg CA1⊗CA20 0 NL/parenrightbigg (P−1 A1⊗P−1 A2)/bracerightbiggD . Then with the definition of the Drazin inverse and the fact that ( CA1⊗CA2)−1=C−1 A1⊗C−1 A2, the base case is complete. (A1⊗A2)D= (PA1⊗PA2)/parenleftbigg C−1 A1⊗C−1 A20 0 0/parenrightbigg (P−1 A1⊗P−1 A2) =A1D⊗A2D. The base case has been established: ( A1⊗A2)D=A1D⊗A2D. In the induction hypothesis, we assume that ( A1⊗A2⊗. . .⊗An)D=AD 1⊗AD 2⊗. . .⊗AD nand show that (A1⊗A2⊗. . .⊗An⊗An+1)D=AD 1⊗AD 2⊗. . .⊗AD n⊗AD n+1. LetY=A1⊗A2⊗. . .⊗An. Then YD= (A1⊗A2⊗. . .⊗An)D=AD 1⊗AD 2⊗. . .⊗AD nby the induction step. And (A1⊗A2⊗. . .⊗An⊗An+1)D= (Y⊗An+1)D =YD⊗AD n+1 =AD 1⊗AD 2⊗. . .⊗AD n⊗AD n+1. 9 The group inverse of A, denoted by A#, is a special case of the Drazin inverse for singular square matrices with index 1. Thus, (A1⊗A2⊗ · · · ⊗ An)#=A# 1⊗A# 2⊗ · · · ⊗ A# n. /square Proof of 3: (Proof by Induction on n) We begin with a proof of the base case, ( A1⊗A2)†=A† 1⊗A† 2. We derive a different expression for the right-hand side, then show that the left-hand side can also be written this way. Every real mA1×nA1matrix A1has a URV factorization A1=UA1/parenleftbigg CA10 0 0/parenrightbigg VT A1, where the orthogonal matrices UA1andVA1are of order mA1andnA1, respectively, the nonsingular matrix CA1is size rA1×rA1andrA1=rank(A1). Similarly, a real mA2×nA2matrix A2can be written as A2=UA2/parenleftbigg CA20 0 0/parenrightbigg VT A2, where the orthogonal matrices UA2andVA2have size mA2×mA2andnA2×nA2respectively, the nonsingular matrix CA2has size rA2×rA2andrA2=rank(A2). The definition of the Moore-Penrose pseudoinverse of A1, denoted A† 1, gives A† 1=VA1/parenleftbigg C−1 A10 0 0/parenrightbigg UT A1, and likewise, A† 2=VA2/parenleftbigg C−1 A20 0 0/parenrightbigg UT A2. Thus, A† 1⊗A† 2=/bracketleftbigg VA1/parenleftbigg C−1 A10 0 0/parenrightbigg UT A1/bracketrightbigg ⊗/bracketleftbigg VA2/parenleftbigg C−1 A20 0 0/parenrightbigg UT A2/bracketrightbigg = (VA1⊗VA2) C−1 A1⊗C−1 A20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (UT A1⊗UT A2) = (VA1⊗VA2)/parenleftbigg C−1 A1⊗C−1 A20 0 0/parenrightbigg (UT A1⊗UT A2). Now we show that ( A1⊗A2)†is equal to the above expression. (A1⊗A2)†=/braceleftbigg/bracketleftbigg UA1/parenleftbigg CA10 0 0/parenrightbigg VT A1/bracketrightbigg ⊗/bracketleftbigg UA2/parenleftbigg CA20 0 0/parenrightbigg VT A2/bracketrightbigg/bracerightbigg† =/braceleftbigg (UA1⊗UA2)/bracketleftbigg/parenleftbigg CA10 0 0/parenrightbigg ⊗/parenleftbigg CA20 0 0/parenrightbigg/bracketrightbigg (VT A1⊗VT A2)/bracerightbigg† =  (UA1⊗UA2) CA1⊗CA20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (VA1⊗VA2)T  † 10 =/braceleftbigg (VA1⊗VA2)/parenleftbigg (CA1⊗CA2)−10 0 0/parenrightbigg (UA1⊗UA2)T/bracerightbigg =/braceleftbigg (VA1⊗VA2)/parenleftbigg C−1 A1⊗C−1 A20 0 0/parenrightbigg (UA1⊗UA2)T/bracerightbigg . The base case has been established. ( A1⊗A2)†=A† 1⊗A† 2. By the induction hypothesis, we have (A1⊗A2⊗. . .⊗An)†=A† 1⊗A† 2⊗. . .⊗A† nand show that (A1⊗A2⊗. . .⊗An⊗An+1)†=A† 1⊗A† 2⊗. . .⊗A† n⊗A† n+1. Let Y=A1⊗A2⊗. . .⊗An. Then Y†= (A1⊗A2⊗. . .⊗An)†=A† 1⊗A† 2⊗. . .⊗A† nby the induction step. Finally, (A1⊗A2⊗. . .⊗An⊗An+1)†= (Y⊗An+1)† =Y†⊗A† n+1 =A† 1⊗A† 2⊗. . .⊗A† n⊗A† n+1. /square 3 Applications of Properties of Kronecker Products To demonstrate the usefulness of applying these properties of Kronecker products, we return to the linear system problem, ( A⊗B)x=d. Let Abem×mandBbem×m. Property 3(a) of section 2.2.2 regarding LUfactorizations can be exploited. If A⊗Bis nonsingular then a solution exists and the system can be written as (LA⊗LB)(UA⊗UB)x=d, where PAA=LAUAandPBB=LBUB. First the lower triangular system ( LA⊗LB)z= (PA⊗PB)dis solved by forward substitution in O(m3) time. Then ( UA⊗UB)x=zis solved by back substitution in O(m3) time. Without exploiting the structure, Gaussian elimination requires O(m6) arithmetic opera- tions. The Kronecker structure also avoids the formation of m2×m2matrices; only the smaller LA,LB, UA,UBare needed. For example, consider the forward substitution ( LA⊗LB)z=d, where AandBare 3 ×3 matrices andzanddare 9×1 vectors. To simplify the notation, we assume PAandPBare identity matrices.  a110 0 a21a220 a31a32a33 ⊗ b110 0 b21b220 b31b32b33  z1 z2 ... z9 = d1 d2 ... d9 . Then LA⊗LB= a11b11 0 0 0 0 0 0 0 0 a11b21a11b22 0 0 0 0 0 0 0 a11b31a11b32a11b33 0 0 0 0 0 0 a21b11 0 0 a22b11 0 0 0 0 0 a21b21a21b22 0 a22b21a22b22 0 0 0 0 a21b31a21b32a11b33a22b31a22b32a22b33 0 0 0 a31b11 0 0 a32b11 0 0 a33b11 0 0 a31b21a31b22 0 a32b21a32b22 0 a33b21a33b22 0 a31b31a31b32a31b33a32b31a32b32a32b33a33b31a33b32a33b33 , 11 which is a unit lower triangular matrix with lower triangular blocks. The first m= 3 equations of this 9×9 system represent a lower triangular matrix and can be solved in O(m2) arithmetic operations. a11b11z1=d1, a11b21z1+a11b22z2=d2, a11b31z1+a11b32z2+a11b33z3=d3. Now the next three equations are: a21b11z1+a22b11z4=d4, a21b21z1+a21b22z2+a22b21z4+a22b22z5=d5, a21b31z1+a21b32z2+a21b33z3+a22b31z4+a22b32z5+a22b33z6=d6. The boldface expression in the first equation, a21b11z1, can be computed as a21d1/a11. The second bold- face expression, a21b21z1+a21b22z2, is just a21d2/a11, while the third expression, a21b31z1+a21b32z2+a21b33z3, isa21d3/a11. We use the previous expressions for obtaining z1,z2andz3in the first set of equations to simplify the second set of three equations. The simplified second set of equations becomes: a22b11z4=d4−a21d1 a11, a22b21z4+a22b22z5=d5−a21d2 a11, a22b31z4+a22b32z5+a22b33z6=d6−a21d3 a11. Solving the second set of equations takes O(m) arithmetic operations and the forward solve step takes O(m2) operations, so obtaining z4,z5andz6takes O(m2) time. This simplification and using the work from the previous solution step continues so that solving each of the msets of mequations takes O(m2) time, resulting in an overall solution time of O(m3). Exploiting the Kronecker structure reduces the usual, expected O(m4) time to solve ( LA⊗LB)z=dtoO(m3) time. One final note regarding the exploitation of the Kronecker structure of the linear system remains. Suppose the matrices AandBare of different sizes. Then, the time required to solve the linear system (A⊗B)x=disO(mAm2 B), where mAis the size of AandmBis the size of B. Van Loan’s paper [23] provides a thorough catalog of further applications of the Kronecker product. We briefly mention a few here. One application receiving growing interest is semidefinite programming. Due to the surge of work on interior point methods, the solution to systems involving the symmetric Kronecker product has been studied recently. The Kronecker product appears in numerous types of least squares problems; one example is the problem of surface fitting with splines. Kronecker products have also been used to unify the field of fast transforms such as the fast Fourier transform, the Hartley transform, and fast wavelet transforms. The Kronecker product plays an instrumental role in many image restoration algorithms. The Kronecker product has also been used to form approximate inverse preconditioners, an application we emphasize in section 4. One application of the Kronecker product not found in Van Loan’s paper is Stochastic Automata Networks (SANs). The remainder of this paper deals with SANs and their connection to the Kronecker product. First, we define SANs, then we discuss their preconditioning problems. 4 Stochastic Automata Networks Markov chains can be used to model many physical systems. For example, Markov chains are used frequently to answer performance questions about parallel and distributed computer systems. While 12 Markov chains provide accurate measures of the system, the size of the Markov chain can quickly grow to an enormous and even intractable size. Storing the state space and infinitesimal generator matrix Qfor such large Markov chains (on the order of millions) has become a bottleneck. One remedy for this storage problem is Stochastic Automata Networks, which store the infinitesimal generator of the Markov chain in compact form using Kronecker products. Stochastic Automata Networks (SANs) [14] are particularly applicable to parallel and distributed computer systems. The reason for this will become clear after we define SANs. A SAN consists of several individual stochastic automata which act independently for the most part. Occasionally, these automata may need to coordinate their actions, thus connecting the individual au- tomata in a network of automata which depend on one another. Each individual automaton A(i)has a number of states associated with it. A(i)also has a number of rules which determine its movement from one state to the next. The state of any automaton at time tis the state it occupies at time t. The state of the collective SAN at time tis the state of each of its corresponding automata. Figure 1 gives the high-level representation of a SAN. Figure 1: Stochastic Automata Network Automaton A(1)contains 4 states, A(2)contains 4 states and A(3)contains 3 states. The current state of each automaton is denoted by the shaded circle. Thus the current state of the SAN is denoted by all three shaded circles. The line connecting A(1)andA(2)represents the interaction between these two automata. Somehow A(1)andA(2)need to coordinate their actions. Exactly how they might need to do this will be described later. If the lines connecting the automata were not present in the diagram thenA(1),A(2)andA(3)would be completely independent systems. A(1)’s stochastic behavior could be modeled with a separate Markov chain from A(2)and so on. Thus, SANs are only useful for automata which have some interaction. However, too much interaction among the automata can complicate the SAN to the point that its use is questionable. Clearly, SANs should be restricted to systems with appropriate infrequent interaction. 13 There are two general ways in which these automata may interact with one another. One concerns the transitions themselves and the other concerns the actual transition rates. First, the global state may change when a transition occurs. Transitions can be either local or synchronizing. Local transitions only affect the corresponding automaton. When an automaton has a local transition, it moves from one of its states to another of its states. Synchronizing transitions are not local. They affect the global state by changing the state of several automata. A synchronizing transition occurs when one automaton enables a transition to occur in two or more other automata. The second type of interaction draws another distinction between transitions. They can be either constant or functional. A functional transition occurs when an automaton’s transition rate is a function of the state of another automaton. Transitions that are not functional are called constant. Constant or functional transitions, unlike synchronizing transitions, affect only the local automata involved. Note that synchronizing transitions may be constant or functional. This information regarding the automata and their types of transitions provides all the information needed to formally define a SAN, as Atif and Plateau have done [16]. While this infrequent interaction (synchronizing transitions and functional transition rates) does complicate SANs, Plateau and her coworkers have shown that the SAN can still be represented in compact form as a sum of Kronecker products, known as the SAN descriptor [14, 17, 20]. Plateau and Fourneau [17] have shown that SANs with Nautomata and Esynchronizing events should be handled by separating out the local transitions for each automata and writing this local effect as the ordinary sum of NKronecker products with each Kronecker product involving Nsmaller matrices. Then the effect of the synchronizing events is added. Each synchronizing event requires two more Kronecker products of Nmatrices. Thus the infinitesimal generator of a SAN can always be written as Q=2E+N/summationdisplay j=1⊗N i=1Q(i) j. When the infinitesimal generator matrix Qis written and stored in the form/summationtext2E+N j=1⊗N i=1Q(i) j, this is called the SAN descriptor . The reader should note that the SAN descriptor is the sumof Kronecker products. Thus far we have only discussed the effect of synchronizing events on the structure of this SAN descriptor. In summary, we learned that there are two more terms in the descriptor for each synchronizing event. An increasing number of synchronizing events increases the complexity of the SAN model, which is why SANs are restricted to systems with appropriate infrequent interaction. We now mention the effect of functional transitions on the SAN descriptor. By extending the ordinary Kronecker product to the generalized Kronecker product [17], [16], the SAN descriptor can still be written as above, but the elements of the Q(i) jmatrices may now be functions. These functional entries require that the appropriate numerical values be computed and substituted each time the functional rate is needed. Thus, functional transitions do not change the structure of the SAN descriptor, but they do add complexity and a computational burden. In practice, the modeler works from the SAN system and forms and stores each of the Q(i) jmatrices following the rules given in [20, 21, 17]. We emphasize that the global infinitesimal generator matrix Qis never formed or stored. Herein lies the storage-saving capacity of the SAN formalism. Consider a collection of four automata each of size 100. Suppose the infrequent interaction among these four automata is described by two synchronizing events and there are no functional transitions. The global Q is size 108but only 2 E+N= 2·2 + 4 = 8 sparse matrices of size 100 need to be stored, thanks to the Kronecker product! 14 4.1 Stationary Analysis of a Stochastic Automata Network The computation of the stationary solution πof a continuous-time ergodic Markov chain involves solving the linear system πQ= 0 and πeT= 1, where Qis the infinitesimal generator of the Markov chain and erepresents the unit row vector. Qis singular with rank n−1. Thus, finding the stationary solution of a continuous-time Markov chain can be viewed as a linear system problem. Another way to view the same problem is as an eigenvalue problem. Pis the transition probability matrix associated with the same system. In fact, P=I+ ∆tQwhere ∆ t≤1 max|qii|.Pis a stochastic matrix with a unit eigenvalue. Then finding the stationary solution πinvolves solving π=πP, which is an eigenvalue problem. Now this eigenvalue problem can be used to define the power method, an iterative method for finding πby computing iterates with x(k+1)=x(k)P. With a suitable initial iterate x(0),x(k+1)will converge to the eigenvector πwhich can then be normalized so that πcontains the stationary solution. Very large Markov chains are often represented as SANs using the SAN descriptor in place of Q. Namely, Q=/summationtextT j=1⊗N i=1Q(i) j, where T= 2E+N,Eis the number of synchronizing events and Nis the number of automata. Since P=I+ ∆tQ, then in the SAN formalism, P=I+ ∆tQ=⊗N i=1Ini+T/summationdisplay j=1∆t⊗N i=1Q(i) j, and the power method for SANs can be written as x(k+1)=x(k)(I+ ∆tQ) =x(k)+ ∆tx(k)(T/summationdisplay j=1⊗N i=1Q(i) j). The power method is the simplest of all iterative methods for finding the stationary solution vector π. The Jacobi, Gauss-Seidel and SOR method are three more iterative methods used for solving linear systems, such as our homogeneous linear system πQ= 0. Yet these methods are based on splittings of the transition matrix and thus are not easily transferable to the SAN formalism. Another class of iterative methods is that of projection methods. These methods approximate an exact solution (in our case, the stationary solution) by building better and better approximations which are taken from small-dimension subspaces. Some popular projection methods are Arnoldi, GMRES, CGS, BiCGSTAB and QMR. Such projections methods can be and have been applied to SANs [1, 15, 22]. In fact, any iterative or projection method which involves a matrix-vector multiply can be used to find the stationary distribution of a SAN. In place of the matrix-vector multiply, the Kronecker product-vector multiplication algorithm invented by Fernandes and his coworkers [3] can be used. Direct methods for solving linear systems, such as those based on LUdecompositions, are not immediately amenable to SANs because the SAN’s compact descriptor representation of the generator matrix precludes easy access to the L,Ufactors. Furthermore, SANs are used as a compact, alternative representation for very large Markov models. The size of such models makes direct methods impractical [20]. 4.2 Preconditioning for Stochastic Automata Networks It is well known that the iterative methods discussed above perform better when preconditioners are used. The convergence of an iterative method depends on the eigenvalues of the system. Any iterative 15 method can converge slowly if the eigenvalue distribution is undesirable for that method. For example, when the subdominant eigenvalue of the iteration matrix is close to the dominant eigenvalue (which is 1 for our transition matrices P), the power method converges slowly. Thus the goal of preconditioning is to modify the eigenvalue distribution of the iteration matrix so that convergence is improved while the solution remains unchanged. In general, for the linear system Ax=b, we introduce the preconditioning matrix M, so that MAx = Mb. We hope that Mis a good approximation of A−1and thus convergence will be rapid. For Markov chain problems, the preconditioned power method becomes x(k+1)=x(k)(I−(I−P)M). Since the matrix ( I−P) is singular with rank ( n−1), we choose Mto be a good approximation of the group inverse of ( I−P), written as ( I−P)#. Thus for SANs, the preconditioned power method is x(k+1)=x(k)(I−(I−(I+ ∆tQ))M) =x(k)(I+ ∆tQM) =x(k)+ ∆tx(k)QM =x(k)+ ∆tx(k)(T/summationdisplay j=1⊗N i=1Q(i) j)M. The problem now becomes that of finding a suitable preconditioner Mthat fits nicely into the SAN formalism. A popular set of preconditioners, ILUpreconditioners, have largely been dismissed from consideration. The problem with adapting ILUpreconditioners to SANs (for use in an iterative method, like the preconditioned power method) is that they are based on incomplete LUfactorizations of the transition matrix. SANs store the transition matrix information as a sum of Kronecker products. And thus, an LUfactorization of a SAN descriptor is not easily accessible. Numerous other preconditioners have been proposed for SANs but each has been unsuccessful [22, 1, 18]. Recently, we discovered a nearest Kronecker product (NKP) preconditioner for SANs [8]. The initial results for the NKP preconditioner look promising [7, 9]. Our NKP preconditioner is derived from Pitsianis and Van Loan’s work on approximation with Kronecker products. They discovered a method for finding the nearest Kronecker product, A⊗B, for a general matrix R[13]. Since A⊗B≈R, one would hope that A−1⊗B−1≈R−1. They took A−1⊗B−1=Mas the preconditioner and tested this on a small example. Their Kronecker preconditioner compared favorably with many other preconditioners. This sparked us to try to extend this to find a suitable SAN preconditioner. For our case of Markov chains, we want to approximate Q#rather than Q−1. However, the algorithm for finding the AandB almost always results in a nonsingular Aand nonsingular B. Thus, we must use the standard inverses, A−1andB−1, to form the preconditioner. In effect, we are using the ideal preconditioner M=A−1⊗B−1 for a nearby system whose coefficient matrix ˆQis almost Q. Finding the small A−1andB−1matrices is not too difficult and the approximation is good for many matrices with nice structure. The advantage of the Kronecker approximation for SANs is that Mneed never be formed, instead only A−1andB−1 need to be stored and used in the vector-Kronecker product multiplication of the iterative methods. In fact, we were able to extend Pitsianis and Van Loan’s work to find any number of smaller matrices whose Kronecker product approximates the original matrix Q[8]. Thus, we can find A, B, . . . , N such that A⊗B⊗ · · · ⊗ N≈Q. We take M=A−1⊗B−1⊗ · · · ⊗ N−1as our NKP preconditioner for SANs. Our initial battery of tests of the NKP preconditioner on SANs reports good results [9]. In fact, the NKP SAN preconditioner outperforms all other current preconditioners. We would like to remind the reader that the properties and power of the Kronecker product made this discovery of a SAN preconditioner possible. 16 5 Conclusion The use and power of the Kronecker product is indeed ubiquitous as Van Loan [23] has suggested. In this paper, we have gathered and cataloged the most useful properties of the Kronecker product and we also added several new properties to this list. We then used these properties to describe a new application for the Kronecker product, Stochastic Automata Networks. References [1] P. Buchholz, Projection methods for the analysis of stochastic automata networks, in B. Plateau, W. J. Stewart, and M. Silva, eds., Numerical Solution of Markov Chains , Prensas Universitarias de Zaragoza, 1999, pp. 149-168. [2] S. L. Campbell and C. D. Meyer, Generalized Inverses of Linear Transformations , Pitman Publish- ing, London, 1979. [3] P. Fernandes, B. Plateau, and W.J. Stewart, Efficient descriptor-vector multiplications in stochastic automata networks, Journal of Association for Computing Machinery , 45:(1998) 381-414. [4] A. Graham, Kronecker Products and Matrix Calculus with Applications , John Wiley and Sons, New York, 1981. [5] H. V. Henderson, F. Pukelsheim and S. R. Searle, On the history of the Kronecker product, Linear and Multilinear Algebra , 14(1983) 113-120. [6] R. A. Horn and C. R. Johnson, Topics in Matrix Analysis , Cambridge University Press, Cambridge, 1991. [7] A. N. Langville, Preconditioning for Stochastic Automata Networks , Ph. D. Thesis, Operations Re- search Program, North Carolina State University, May 2002. [8] A. N. Langville and W. J. Stewart, A Kronecker product approximate inverse preconditioner for SANs, Journal of Numerical Linear Algebra with Applications , to appear. [9] A. N. Langville and W. J. Stewart, Testing the NKP preconditioner on MCs and SANs, INFORMS Journal on Computing , to appear. [10] C. C. MacDuffee, The Theory of Matrices , Chelsea, New York, 1946. [11] C. D. Meyer, Matrix Analysis and Applied Linear Algebra , SIAM, Philadelphia, 2000. [12] C. D. Meyer, The role of the group generalized inverse in the theory of finite Markov chains, SIAM Review , 17(1975) 443-464. [13] N. Pitsianis and C. Van Loan, Approximation with Kronecker products, in M. S. Moonen and G. H. Golub, eds., Linear Algebra for Large Scale and Real Time Applications , Kluwer Academic Publishers, 1993, pp. 293-314. [14] B. Plateau, On the stochastic structure of parallelism and synchronization models for distributed algorithms, Performance Evaluation Review , 13(1985) 142-154. [15] B. Plateau, PEPS: A package for solving complex Markov models of parallel systems, in R. Puigjaner and D. Potier, eds., Modelling Techniques and Tools for Computer Performance Evaluation , Plenum Press, New York, 1990, pp. 291-306. 17 [16] B. Plateau and K. Atif, Stochastic automata network for modelling parallel systems, IEEE Trans- actions on Software Engineering , 17:10(1991) 1093-1108. [17] B. Plateau and J. M. Fourneau, A methodology for solving Markov models of parallel systems, Journal of Parallel and Distributed Computing , 12(1991) 370-387. [18] B. Plateau, W.J. Stewart, and P. Fernandes, On the benefits of using functional transitions in Kronecker modelling, submitted to Performance Evaluation , Jan. 2001. [19] W-H. Steeb, Kronecker Product of Matrices and Applications , Wissenschaftsverlag, Mannheim, 1991. [20] W. J. Stewart, Introduction to the Numerical Solution of Markov Chains , Princeton University Press, Princeton, 1994. [21] W. J. Stewart, Stochastic automata networks, In W. K. Grassmann, ed., Computational Probability , Kluwer Academic Press, Boston, 2000, pp. 113-152. [22] W. J. Stewart, K. Atif, and B. Plateau, The numerical solution of stochastic automata networks, European Journal of Operational Research , 86(1995) 503-525. [23] C. F. Van Loan, The ubiquitous Kronecker product, Journal of Computational and Applied Mathe- matics , 123(2000) 85-100. 18