ALGORITMI ZA REINVERZIJO BAZNE MATRIKE V LINEARNEM PROGRAMU Janez Barle, Janez Grad Univerza Edvarda Kardelja v Ljubljanl, Ekonomska fakulteta Borisa KidriSa, Ljubljana UDK: 519.852.11 V prlspevku oplsemo nekatere algoritme za invertlranje matrik, ki so uporabni pri reSevanju linearnih programov. Njihova značilnost je specifiSen naSin predstavitve lnverzne matrike in razdelitev postopka v dve fazi. Poseben poudarek dajemo probleraom, ki se nanaSajo na rafiu- nalniško realizacijo algoritmov. BASIS MATRIX REINVERSION ALGORITHMS IN^LINEAR PROGRAMMING. In this paper we describe some of the algorlthms for matrix reinversion whlch can be applied in the procedures for solvlng the linear prograiraning problem. They use a specific presentatlon of the inverse matrix and two-phase computation process. Special attention is paid to those problems associated with the computer processing of these algorithms. UVOD ReSevanje slstema linearnlh ena£b in s tem povezan problem invertiranja matrike spada- ta roed klasične probleme numerične linearne algebre. Običajne metode za reševanje line- arnlh enačb uporabljajo Gaussov eliminaclj- ski postopek, pri katerem se zaradl zagotav- ljanja numerične stabilnosti lzvaja t.i. pi- votlranje (delno ali popolno). Pri delu na programski opremi za linearno programtranje smo prišli do spoznanja, da omenjene metode niso vedno učinkovite. Zaradi lastnosti ma- trik, ki se pojavljajo pri praktlčnih pro- blemih llnearnega programiranja, je treba za njihovo invertlranje uporabiti posebne postopke. Uporaba teh postopkov lahko pre- cej povefia hltrost in natančnost pri re5e- vanju linearnlh programov velikih dimenzij. Zato je podprogram za invertiranje matrike ena od najbolj pomembnih sestavin komercial- nih programsklh paketov za linearno progra- miranje. V nadaljevanju podajamo opis neka- terih najbolj znanih algoritmov za lnverti- ranje matrlke pri llnearnera programiranju. Opisujemo tudi nekatere izkuSnje, kl smo jih pridobili pri programiranju in testira- nju teh algoritmov. Problem linearnega programiranja lahko v matričnl obliki oplSemo takole: T mlnimiziraj c x pri pogojih Ax = b in x > 0 kjer je A realna matrika dimenzije m*n, b, c in x pa realni vektorji ustreznih dlmen- zlj. Optimalna rešitev llnearnega programa je tak vektor x, ki reši gornji problem. UpoStevati moramo, da prl večini praktiSno pomembnih problemov linearnega programira- nja velja: 1) Matrika A je zelo vellkih dimenzij (tudi do nekaj tisoč vrstic in stolpoev) 2) Matrika A je razpršena. To pomeni, da je delez njenih neniCelnih elementov zelo maj- hen (pogosto tudi manj od 1%). Zaradl tega se v praksi za reševanje linear- nih programov skoraj vedno uporablja t.i. revidirana metoda simpleksov, kl jo je možno prllagoditi delu z matrikami velikih dimenzlj. To je iterativna metoda, ki je podrobno opisana v literaturi (npr. |l|). Omenimo le, da se na vsakem koraku te meto- de reSuje zaporedje sistemov llnearnlh enačb irB UB' Bv u in kjer je B nesingularna kvadratna podmatrika matrlke A (imenujemo jo bazna matrika), c^ in xR sta ustrezna podvektorja vektorjev c in b, u pa je stolpec matrike A, ki je doloien z rešitvijo prvega sistema enačb. ReSevanje gornjih sistemov enačb olajšuje dejstvo, da se bazni matrlkl pri dveh zapo- rednih korakih algoritma razllkujeta le v enem stolpcu. PREDSTAVITEV INVERZNE BAZNE MATRIKE Revidirana tnetoda simpleksov ima vefi inačlc, kl izhajajo lz razlifinih načlnov predstavitve ma- trike B (all B-l). Predstavitev matrlke B"1 v eksplicitni obliki lahko uporabimo le pri pro- blemih majhnih dimenzlj. V obstoječlh program- skih paketih je najbolj obifiajna predstavitev lnverzne bazne matrike v obliki produkta t.i. elementarnih matrik: B-1 - EkEk-l ••• ,E1 Elementarne matrike E^ kujejo od enotne matrlke le v enem stolpcu, ta- ko da zadostuje shranlti nenifielne elemente tega stolpca in podatke o tem, kje se le-ti nahajajo. Ena od prednosti produktne predstar vitve inverzne bazne matrlke je tudi enostaven način njenega ažuriranja. Po vsakem koraku re- vidirane metode simpleksov se v produkt dodaja nova eleraentarna matrika, stare elementarne 1, ... k) se razli- 145 matrike pa ostanejo nespremenjene, Tako pred- stavltev matrlke lmenujemo tudi ETA-datoteka, tisti stolpec, v katerem se elementarna matri- ka razlikuje od enotne pa imenujemo ETA-vektor. Znano je, da pri reSevanju linearnih enaSb po Gaussovem elitninacijskem postopku računanje inverzne matrike ni potrebno. BoljSa pot za reSevanje sistema enačb (Bx = b) je izvrSiti razcep permutirane matrike na spodnjo in zgor- njo trikotno matriko {PB = LO) in potem reše- vati dva trikotna sistema enačb (Ly = b in 0x = y). Ta ugotovitev v določeni meri velja tudi za bažno matriko linearnega programa. Za njo lahko uporabimo razcep B = LU, kjer vsa- kega od obeh faktorjev predstavimo z zapored- jem ustreznih trikotnih elementarnih matrik: L = L1L2...LkinU=Ok...O1 Tako predstavitev intenujemo eliminacijska ob- lika inverzne matrike. V zapisu smo zaradi poenostavitve izpustili permutacijske matri- ke, ki so podane z zaporedjem, v katerem ob- ravnavamo vrstice in stolpce matrike B. Iz- kaže se, da v splošnem uporaba eliminacijske oblike matrike omogoča ekonomičnejšo izrabo pomnilnika in Sasovno hltrejSi izrafiun reSit- .ve linearnega programa v primerjavi s postop- kom, ki uporablja produktnp obliko inverzne matrike. Vseeno so trikotni razcep zafiell vgrajevati v komercialne programske pakete Sele potem, ko so bili razviti ufiinkoviti po- stopkl za ažuriranje trikotnih faktorjev (npr. metoda Forresta in Tomlina, glej 111) . Ugotav- ljamo pa, da uporaba trikotnega razcepa zahte- va precej bolj zapletene algoritme in podat- kovne strukture kot postopek, ki uporablja produktno obliko inverzne matrike. REINVERZIJA BAZNE MATRIKE Ker se med izvajanjem simpleksnega algoritma , število elementarnih matrik v predstavitvi B~ ali v trikotnih faktorjih L in U povečuje, se izkaže za nujno občasno izračunati novo pred- stavltev B"1 ali novi trikotnl razcep matrike B. Ta postopek, lmenujemo ga reinverzija baz- ne matrike, sproSSa pomnllnik ln omogofia hi- trejše in natančnejše izvajanje algoritraa na preostallh korakih revidirane metode simplek- sov. Nazlv relnverzlja je primeren tudi za po- stopek dolo&anja novih elementarnlh matrik trikotnega razcepa matrike, ker je ta naloga v tesnl povezavi z invertiranjem matrlke. Ta- ko lahko iz znane eliminacijske oblike matrike takoj dobimo produktno obliko inverzne matrike: B -1 š1 .Za produktno in eliminacijsko obllko matrike uporabljamo tudl skupen naziv ETA-faktoriza- cija matrike. Omenlmo še, da je inverzura ele- mentarne matrlke zopet elementarna matrika. Vzemimo, da se prvotna elementarna matrika raz- llkuje od enotne matrike v stolpcu, kl ga ozna- Slmo z y. Pri inverzumu se ta stolpec spremeni po naslednjih pravilih: diagonalni element yk postane l/v^r vsi ostali elementi y^ pa se spremenijo v "V Pri relnverzlji lzhajamo iz prvotne matrike B, ki je praviloma razpršena. Lahko se zgodi, da ima matrika B"1 mnogo veS neničelnih elementov kot B. To velja tudi za skupno Stevllo neniSel- nih elementov v ETA-faktorizaciji inverzne ma- trike. Ta pojav,' imenujemo ga napolnitev, je posebno izrazit takrat, ko pri izvajanju Gaus- sovega eliminacijskega postopka uporabimo del- no ali popolno pivotiranje. To sta postopka, pri katerih se teži k izbiri ključnega elemen- ta (pivota) s Cim vefijo absolutno vrednostjo ("pivoting for size"). Zato je pri relnverziji bazne matrlke treba uporabiti take kriterije za izbor pivotov, ki vodijo k ETA-faktorizaci- ji s Sim manjSim številom nenlSelnih elementov ("plvoting for sparslty"). Poskusi optimalne reSltve tega problema, ki je znan tudi pod ime- nom "problem minimalne napolnitve", so pripe- ljall do ugotovitve, da je ta problem NP-poln. To pomeni, da je malo upanja, da bo za njegovo reSevanje kdaj oblikovan ekonomičen algoritem (glej |1O|). Zato so blli razvlti številni al- goritml, osnovani na razliSnih heurističnih pravilih za izbor pivota, ki so bolj ali manj uspešno reševali zastavljeno nalogo. Pomemben dogodek v raziskovanju tega problema je bil, ko je leta 1957 H. Markowitz objavil članek o tej problematiki (glej |8|). V njem je posta- vil nekaj trdltev, ki so spodbudile Stevilne poznejše raziskave: 1) Postopek invertiranja je smiselno razdeliti v dve fazi. V prvl fazi, imenujemo jo ana- litična faza, se vnaprej doloSi zaporedje, po katerem se bo vršilo pivotiranje. Pri tem se uporabljajo le informacije o tem, katerl elementi so neničelni, ne oziraje se.na njihovp vrednost. V drugi fazi, lmenujemo jo numerična faza, se na podlagi vnaprej do- ločenega zaporedja pivotov lzvrSi dejanski izračun nove ETA-faktorizacije. 2) Dobre rezultate pri minlmizaciji nenifielnih elementov v ETA-faktorizaciji lahko doseže- mo tako, da na posaraeznem koraku izbiramo tlsti pivot, prl katerem lma "Markowitzevo Stevilo" mlnimalno vrednost. Z v, ("vrstlčnl števec") ln 3j ("stolpifini Stevec") sta oznaSeni 5te-- vili neniSelnih eleraentov v 1-tl vrstici in j-tem stolpcu, pri Cemer upoštevamo le tiste vrstice in stolpce, v katerih še ni bilo iz- vršeno pivotiranje. Markovitz je bil tudi prvi, ki je ugotovil pri- mernost elimlnacijske oblike za predstavitev bazne matrike pri linearnera programiranju. ANALITIČNA FAZA REINVERZIJE Ideja o ločenem izvajanju analitiSne faze se ni uveljavila le pri reinverzlji bazne matrike, temveS tiidi pri vseh drugih problemih numerične materaatike, ki zahtevajo reSevanje razpršenih sistemov linearnih enaSb. Praksa je potrdila tudi to, da uporaba Markowitzevega pravlla pri izbiri pivota vodi k zelo razprSeni ETA-fakto- rizacijl inverzne matrike. Vseeno je dlrektna uporaba tega pravila otežkoSena, ker zahteva hkraten dostop do vrstic in stolpcev matrike. Kolikor nam je znano, uporablja Markowitzevo pravilo v Slsti oblikl le podprogram za reln- verzijo pri IBM-ovem programskem paketu MPSX/370 (glej |2|). Izkazalo se je, da je mož- no uporabiti različne kriterije, kl na posreden naSin merijo Markoviritzevo števllo, hkrati pa zahtevajo veliko manj rafiunanja. Sodobni programski paketl v analitifini fazi re- inverzije veSinoma uporabljajo algoritem Hel- lermana in Rarlcka (glej |6| in |7|). Ta algo- ritem je prirejen strukturi matrik, kl so zna- fiilne za praktifine probleme linearnega progra- miranja, s čimer lahko pojasnimo njegovo učin- kovitost. Cilj tega algorltma je, da se z zame- njavo vrstlo ln stolpcev matrika B prevede v t.i. HR matrlko, ki iraa strukturo kot na sllki 1, 146 B1 /////// /////// iiiniiiimn (IIIIIIIIIIIII B3 Slika 1 kjer so B in B spodnje trikotne raatrike z ne- ničelnimi elementi na diagonali. Lahko privza- memo, da ima B , če ni prazna, v vsaki vrsticl ln stolpcu vsaj dva nenifielna elementa, ker bi v nasprotnem primeru lahko razširili B ali B . Matriko B imenujemo jedro all izboklina ("bump"). Algoritem Hellermana in Rarlcka se nadaljuje z ustreznimi operacijaml na jedru, dokler jedro ne postane strukturirano kot na sliki 2, kjer so G ali prazne ali spodnje tri- kotne matrike z neničelnimi elementi na diago- nali: F1 //// G1 ////////// ///// ///// ///// ///// /////; /////> /////> /////> F2 'ii/i '//// i/iir, I/I/I, G2 '//// >//// ////////////////////77777 ///////////////////////// FP Slika 2 Matrike F^ lmenujemo zunanje izbbkline ("exter- nal bumps"). Imajo vsaj po en stolpec z neničel- nimi elementi nad diagonalo. Take stolpce lme- nujemo konice ("spikes"). Zunanje lzbokline se morajo odlikovati z naslednjimi lastnostmi: i) zadnji stolpec v izboklini je konica z ne- ničelnim elementom v najvišji vrstlcl, ii) stolpoi, ki niso konice, morajo imeti neni- čelne diagonalne elemente. Aigpritma ne bomo podrobneje opisovall. Omenimo le, da se konice izbirajo na temelju stolpiSne števfine funkcije tk(j), kjer tj,(j) = število nenl8elnih elementov j-tega stolpca v vrsticah z vrstičnim štev- cem manjšlm ali enakim k Algoritem Hellermana in Raricka pokaže tudi za- nimivo povezavo med teorijo grafov in linearnim programiranjem. Binarno matriko, s katero dela- mo v analitični fazi reinverzije, lahko inter- pretiramo kot matriko sosednosti U9merjenega grafa. V jeziku teorije grafov pomeni iskanje bločne trikotne obllke matrike (to je oblika matrike na sliki 2) določanje močno povezanih komponent usraerjenega grafa. To je problem, ki je spodbudil precejšnjo pozornost razlskoval- cev s področja teorlje grafov. Mi smo za reše- vanje te naloge sprogramlrali eno varianto t. i. Tarjanovega algoritma, ki smo jo povzeli po |3|. Zanimivo je, da je uporaba metod iz teo- rlje grafov pripeljala do izboljšane variante algoritma Hellermana in Rarioka, ki je v lite- raturi dobila ime P4 ("Partitioned Preassigned Pivot Procedure"). S to varlanto se praviloma dobi večje Stevllo zunanjlh izboklin ln manjSe število konic kot s starejSo varianto algorit- ma, ki je znana pod nazivom P3 ("Preassigned Pivot Procedure"). Na koncu naj omenlmo, da so bili na temelju . algoritma Hellermana in Raricka razviti še ne- kateri zanimlvi algoritmi, kot je npr. algori- tem opisan v |9|. NUMERIČNA FAZA REINVERZIJE V numerični fazi reinverzije se, na temelju zaporedja plvotov doloSenega v analitiSni fa- zi, računa nova ETA-faktorizacija inverzne bazne matrike. V tej fazi moramo razllSno ob- ravnavati produktno in ellminacijsko obliko inverzne matrike. kar nl bilo potrebno v ana- litiSni fazi. V |4| ln |5| je podan zelo lz- črpen pregled različnih algoritmov za izvaja- nje numerične faze reinverzije. V nadaljevanju podajamo opis algoritma za produktno obliko inverzne bazne matrike, ki smo ga spogramira- 11 v |ll| . Predhodno določeno zaporedje plvotov je poda- no z vektorjem C, ki vsebuje lndekse stolpcev, in z vektorjem R, kl vsebuje lndekse vrstic. Razen tega potrebujemo še indekse konic ln podatke o tem, kje se začnejo posamezne zuna- nje lzbokllne. Zaradi prihranka pomnilniškega prostora je te lnforraaclje koristno vključitl v vektorja R ln C. To lahko storimo tako, da z negativnim predznakora shranlmo v C indekse stolpcev, ki so konlce, in v R lndekse vrstic, kjer se začnejo posamezne zunanje izbokllne. Nl: |lnlcializacija| a) Vpeljemo naslednje parametre: (I) 1=1 za lndekslranje pivotov, (II) 1=0^ za stolpični lndeks pivota, (ili)r=|RjJ za vrstični indeks plvota, (iv) p=l za indekslranje začetka zuna- nje izbokline. b) Rezerviramo prostor za (I) y - realni vektor dimenzlje m, (II) ETA-datoteko (predstavltev matrik i) c) Deflnlramo y ke B). (1-tl stolpec matrl- N2: |Določanje novega ETA-vektorja| ( l/yr za k=r Vektor z je r-tl stolpec elementarne matrlke E^. N3: |Testlranje konca algorltma| a) če je i=m nadaljujemo na N7 (konec algo- ritina) b) 6e je i 0 deflniramo y = B,^ ln nadalju- jemo na N2 N5: |Transformacija konlce| Izračunamo V (E ustreza prvemu stolpcu v izboklini, kl prlpada Bt^) N6: ]zamenjava konice, Se je plvot enak nlč| a) če je y ?*0 postopek nadaljujemo na ko- raku N2 b) PrelšSemo vse konice s stolpičniml in- deksl 1 |C. ki se nahajajo 147 N7: v lsti zunanji izboklini, dokler ne naj- demo takšne za katero je r-ta komponenta E. ,...E^B* 1—1 k w neničelna = E. , . . .E B, i-1 p ' j c) Izračunamo y d) Zamenjamo C^ in C. e) Postopek nadaljujemo na koraku N2 |Konec| Določanje ETA-datoteke je konSano. Kot vidlmo, je korak N5 potreben le za tiste stolpce, ki so konice. Edino v teh stolpcih se lahko pri relnverziji poveča število neničelnih elementov. To je razlog, da je algoritem Heller-. roana in Raricka tako pomemben. Preverjanje ali je pivot enak niS, ki je potrebno na koraku N6, se v praksl zamenjuje z testi, kl zagotavljajo numerično stabilnost postopka. Mathematical Prograinming, 12, pp. 260-278, 1977. 10) PetkovSek M., "NP-polni probleml", v: M. PetkovSek, Pisanski.T., Izbrana poglavja iz računalništva 1. del, Ljubljana: DMFA SRS, 1982. 11) Barle J., J. Grad, V. Rupnik, "Programska oprema za rešitev linearnega programa z razSiritvami", Razlskovalni projekt po na- ročilu Iskre-Delta, Ljubljana: RCEF, Eko- nomska fakulteta Borisa Kidriča v Ljublja- ni, 1984. ZAKLJUČEK Uporaba matematike v praksi često zahteva poz- navanje metod, ki jih ni možno nau&itl iz stan- dardnih učbenikov. Ta ugotovitev velja tudi za metode, ki smo jlh morali vkljufiiti v podpro- gram za invertiranje bažne matrike linearnega prograraa. Naša raziskovalna skupina nadaljuje s spremljanjera tega raziskovalnega področja in razvojem ustreznih rešitev v sklopu programske oprenie za llnearno programiranje. LITERATORA 1) Bastian M., Lineare Optimierung grosser Systeme: Compact Inverse Verfahren und Basisfaktoriesirung, Konigstain: Verlag Anton Hain, 1980. 2) Benichou M., J.M. Gauthier, G. Hentges, G. Ribiere, "The Efficient Solution of Large- Scale Linear Prograiraning Problems - Some Algorithmic Techniques and Computational Results", Mathematical Progranming, 13, pp. 280-322, 1977. 3) Gustavson F., "Finding the Block Lower Triangular Form of a Sparse Matrix", v: Bunch J.R., D.J. Rose (eds.), Sparse Matrix Computatiorts, New York: Academic Press, 1976. 4) Helgason R.V., J.L. Kennington, "Spike Swapping in Basis Reinversion", Naval Re- search Logistic Quarterly, 27, pp. 697- 701, 1980. . 5) Helgason R.V., J.L. Kennington, "A note on Splltting the Bump in an Elimination Fa- ctorisation", Naval Research Logistics Quarterly, 29, pp. 169-178, 1982. 6) Hellerman E., D. Rarick, "Reinversion with the Preassigned Pivot Procedure", Mathema- tical Programning, 1, pp. 195-216, 1971. 7) Hellerman E., D. Rarlck, "The Partitioned Preassigned Pivot Procedure", v: Rose D.J., R.A. Willoughby (eds.), Sparse Matrlces and Thelr Applicatlons, New York-London: Plenum Press, 1972. 8) Markovitz H., "The Elimination Form of the Inverse and its Applicatlon to Linear Pro- granuning", Management Science, 3, pp. 255- 269, 1957. 9) Lin T.D., R.S.H. Mah, "Hierarchical Parti- tion - a New Optimal Pivoting Algorithm",