IE172 – Algorithms in Systems EngineeringLecture 35Pietro BelottiDepartment of Industrial & Systems EngineeringLehigh UniversityApril 13, 2009Matrices: sparse an d dense◮Density of a matrix: percentage of entries that are 6= 0◮Dense vectors can simply be stored in an array.◮Dense matrices can be stored in a 2-dimensional array:class SparseMatrix {double**rows_; // rows_ [3] points to vector// containing the fourth rowint nRows_;int nCols_;};◮In practical applications, matrices are usually s parse.e.g. in Linear Programming, instances with ≥ 10 nonzeros percolumn are rare, even if tens of thousands of rows◮Space complexity: Θ(mn): one double per cellExampleA =4 1 0 0 20 3 0 0 05 0 0 9 08 0 0 0 4double**A = new double*[4];for (int i=0; i<4; i++)A [i] = new double [5];A[0][0]=4; A[0][1]=1; A[0][2]=0; A[0][3]=0; A[0][4]=2;A[1][0]=0; A[1][1]=3; A[1][2]=0; A[1][3]=0; A[1][4]=0;A[2][0]=5; A[2][1]=0; A[2][2]=0; A[2][3]=9; A[2][4]=0;A[3][0]=8; A[3][1]=0; A[3][2]=0; A[3][3]=0; A[3][4]=4;Sparse matrix storage◮Sparse matrices can be stored using several strategies◮We will learn threeThree Sparse Matrix Formats◮Yale Spars e Matrix Format (Triples)◮(Compressed) Sparse Column Format◮(Compressed) Sparse Row For matYale sparse matrix formatStores the matrix in three arrays:◮A: double array holding the element values◮IA: int array with the row indices of the non-zero values◮JA: int array with the column indices of the non-zeros◮The arrays A, IA, and JA all have lengt h equal t o thenumber of non-zero elements in the matrix◮Is this the “sparsest” way to h old a matrix?◮Does it support efficient operations that we w ould like todo (like inner p roduct of two columns?)◮Space complexity: if z is # of nonzeros, Θ(z) and O(mn); atmost one double per cell, but in general ≪ kmnExampleA =4 1 0 0 20 3 0 0 05 0 0 9 08 0 0 0 4double*A = new double [8];double*IA = new int [8];double*JA = new int [8];A[0]=4; IA[0]=0; JA[0]=0; A[1]=1; IA[1]=0; JA[1]=1;A[2]=2; IA[2]=0; JA[2]=4;A[3]=3; IA[3]=1; JA[3]=1;A[4]=5; IA[4]=2; JA[4]=0; A[5]=9; IA[5]=2; JA[5]=3;A[6]=8; IA[6]=3; JA[6]=0; A[7]=4; IA[7]=3; JA[7]=4;Compressed sparse column format◮Again stores the matrix in three arrays, but the arrays herehave different meanings:◮matval: double array holding the element va lues◮matind: int array holding the row indices of the nonzeroentries in each column.◮matbeg: int array with the index, in the matval andmatind arrays, of thefirst element of ea ch column◮matval and matind: each h ave length equal to number ofnon-zeros◮matbeg: has length one more than the number of columnsin t h e matrixExampleA =4 1 0 0 20 3 0 0 05 0 0 9 08 0 0 0 4Let’s use A for matval, I for matind, and B for matbeg.double*A = new double [8];double*I = new int [8];double*B = new int [6];A[0]=4; I[0]=0; B[0]=0;A[1]=5; I[1]=2;A[2]=8; I[2]=3;A[3]=1; I[3]=0; B[1]=3; B[2]=5; // empty column (==B[3])A[4]=3; I[4]=1;A[5]=9; I[5]=2; B[3]=5;A[6]=2; I[6]=0; B[4]=6;A[7]=4; I[7]=3; B[5]=8; // to mark the last elementCompressed sparse row format◮Like Compressed Column Format, except “row-wise”◮matval: double array holding the element va lues◮matind: int array holding the columns indices of thenonzero entries in each column.◮matbeg: int array holding the index, in the matval andmatind arrays, for the first element of each row◮matval and matind: each h ave length equal to number ofnon-zeros◮matbeg: has length one more than the number of rows inthe matrixExampleA =4 1 0 0 20 3 0 0 05 0 0 9 08 0 0 0 4Let’s use A for matval, I for matind, and B for matbeg.double*A = new double [8];double*I = new int [8];double*B = new int [5];A[0]=4; I[0]=0; B[0]=0;A[1]=5; I[1]=1;A[2]=8; I[2]=4;A[3]=1; I[3]=1; B[1]=3;A[4]=3; I[4]=0; B[2]=4;A[5]=9; I[5]=3;A[6]=2; I[6]=0; B[3]=6;A[7]=4; I[7]=4; B[4]=8; // to mark the last elementAnother look at matrix multiplica tionImportant NotationIf A ∈ Rm×n, then Ajis the jthcolumn, and ajis the jthrow.◮If A ∈ Rm×k, B ∈ Rk×n, then [AB]ij= aTiBj.i.e. the product AB is a matrix; its (i, j) element is equal to theinner product of the i-th row of A with the j-th column of B.◮Also: [AB]ij=ℓXk=1aikbkj.◮Naturally this is only defined if A ∈ Rm×ℓand B ∈ Rℓ×nNote: If two matrices A and B are composed of blocks of the rightsizes, we can compute AB by multiplying the blocks:A BC DE FG H=AE + BG AF + BHCE + DG CF + DH#cols of A, C = #rows of E, F, and#cols of B, D = #rows of G, HMatrix multiplication: linear combinations ofcolumnsLooking at it another way, write B as its columns:B = (B1B2· · · Bn)Then the jthcolumn of AB is ABj, orAB = A (B1B2· · · Bn) = (AB1AB2· · · ABn)Each column of AB is alinear combination of the columns of A.The multipliers for the linear combination are g iven in thecolumn Bj.Matrix Multiplication: Linear Combi nations ofRowsWe can also express the relationship in terms of the rows of AA =a1a2...am, AB =a1Ba2B...amBThe ithrow of AB is a linear combination of the rows of B, withthe weights in the combination coming from t h e weights in ai.Another nice formula:(AB)T= BTATSome definit i ons you already knewVectors {A1, A2. . . An} are linearly dep endent if the zero vectorcan be written as a non-zero linear combination of the vectors.That is, ∃α1, α2, . . . , αnnot all equal to zero such thatnXj=1αjAj= 0.Alternatively, if Ajare columns of a matrix A, then the Ajarelinearly dependent if Az = 0 for some z 6= 0More definitions you a lready knewVectors {A1, A2. . . An} are linearly independent if the zerovector cannot be written as a non -zero linear combination ofthe ve ctors:nXj=1αjAj= 0 ⇒ α1= α2= · · · = αn= 0Alternatively, if Ajare columns of a matrix A, then the Ajarelinearly independent if Az = 0 ⇒ z = 0 (i.e., 0 is the onlysolution to Az = 0).More DefinitionsThe Range of a matrix A ∈ Rm×n, denoted R(A), is the set of alllinear combinations of the columns of A. Thus, R(A) ⊆ Rm.◮b ∈ R(A) ⇔ ∃x ∈ Rnwith Ax = b◮The rang e is sometimes called the column space of A◮R(AT) ⊆ Rnis sometimes called the row space of A.◮Two vectors x, y are orth ogonal if xTy = 0◮The set of all (m-dimens ional) vectors or thogonal t ovectors in R(A) is thenull space of AT◮z ∈ N (AT) ⇔ ATz = 0◮zTb = zT(Ax) = xTATz = 0 if b ∈ R(A), z ∈