Pràctiques de Matlab per àlgebra lineal

Pr`actiques de Matlab per ` Algebra lineal Jaume Amor´os UPC Barcelona 2025/3/27 2 ´ Index 1 Discussi´o i resoluci´o de sistemes lineals. 7 2 Interpolaci´o per splines c´ubics 13 3 Sistemes lineals, i an`alisi del gas natural. 15 4 Tractament d’imatges. 23 5 Xarxa neuronal de tipus MultiLayerPerceptron (MLP) per a aproximar les arrels d’un polinomi 29 3 4´ INDEX Introducci´o Aquest treball ´es un recull de pr`actiques d’ordenador per a fer c`alculs d’` Algebra lineal. Estan pensades per a resoldre amb el programa Matlab o algun dels seus clons de domini p´ublic com ara Octave o Scilab. Busquen ensenyar a l’estudiantat com fer amb l’ordenador els c`alculs m´es b`asics de la materia, i mostrar aplicacions t´ıpiques i/o interessants d’aquests c`alculs, que siguin a m´es realistes. Aquestes pr`actiques s’han oferit als alumnes de l’assignatura d’` Algebra Lineal de primer curs del Grau d’Enginyeria en Tecnologies Industrials de l’ETSEIB, UPC, Barcelona, els anys 20225. S´on adequades per a qualsevol curs d’` Algebra lineal, i tamb´e per a il·lustrar les aplicacions de les Matem`atiques als estudiants d’Enginyeria. Agra¨ıments: Toni Sus´ın ha assesorat durant anys a l’autor en la preparaci´o de pr`actiques de Matlab pels estudiants de l’ETSEIB. Els professors i alumnes de l’assignatura d’` Algebra Lineal de l’ETSEIB dels anys 2022-5 han resolt aquestes pr`actiques i han proporcionat informaci´o valuosa per a millorar-les. Jaume Amor´os Barcelona, 26 de Mar¸c del 2025. 5 6´ INDEX Cap´ıtol 1 Discussi´o i resoluci´o de sistemes lineals. MATLAB t´e una capacitat formidable per a resoldre sistemes d’equacions lineals. En un PC dom`estic pot resoldre en algun minut sistemes 4000 ×4000 amb tots els coeficients no nuls . . . Per`o l’usuari necessita una mica de formaci´o per a entendre els resultats que proporciona el programa. La comanda de MATLAB per a resoldre un sistema lineal ´es x=A\b On A´es la matriu de coeficients del sistema lineal, i b´es el vector columna de termes independents. La resposta del programa ´es un vector columna xamb els valors de les inc`ognites. Mirem exemples per a veure qu`e cont´e exactament aquest vector: Exemple 1: Considerem el sistema lineal 2x−y= 3 x+ 2y= 4 Podem comprovar a m`a que t´e una soluci´o ´unica x= 2, y = 1. Resolem-lo amb Matlab: % matriu de coeficients A=[2 -1 ; 1 2] % terme independent b=[3;4] % resolem el sistema (anomenem xy al vector d’incognites) xy=A\b La resposta que d´ona el programa ´es xy = 2 1 Per tant ens ha retornat la soluci´o exacta, en un vector columna amb components les inc`ognites x, y. Exemple 2: Ara considerem el sistema lineal amb una sola equaci´o x−y+ 2z= 3 7 8CAP´ ITOL 1. DISCUSSI ´ O I RESOLUCI ´ O DE SISTEMES LINEALS. Podem comprovar a m`a que les seves solucions depenen de 2 par`ametres, i formen tot un pla en R3: la soluci´o general ´es  x y z  = 3 0 0  +y 1 1 0  +z −2 0 1  = 3 0 0  + 1−2 1 0 0 1  ·y z, on (3,0,0) ´es una soluci´o particular del sistema (un punt del pla), y, z poden tenir qualsevol valor real, i els dos vectors (1,1,0) i (−2,0,1) s´on vectors directors del pla de les solucions. Resolem el sistema amb Matlab: % matriu de coeficients A=[1 -1 2] % terme independent b=3 % resolem el sistema (anomenem xyz al vector d’incognites) xyz=A\b La resposta que d´ona el programa ´es xyz = 0 0 1.5000 MATLAB ens ha tornat una soluci´o ´unica! Podem comprovar f`acilment que ´es una soluci´o correcta, per`o el programa no ens avisa de que n’hi han infinites m´es. Exemple 3: Considerem el sistema lineal    2x−y= 3 x+ 2y= 4 x+y= 5 Aquest sistema no t´e soluci´o (t´e les dues equacions de l’exemple 1, seguides d’una tercera que la soluci´o de les dues primeres no compleix). Resolem-lo amb Matlab: A=[2 -1 ; 1 2; 1 1] % terme independent b=[3;4;5] % resolem el sistema (anomenem xy al vector d’incognites) xy=A\b La resposta que d´ona el programa ´es xy = 2.2857 1.2857 MATLAB ens ha tornat una soluci´o? El programa sap perfectament que no ho ´es. Si fos una soluci´o exacta es compliria que l’operaci´o matricial % prova de la solucio r=A*xy-b 9 hauria de ser un vector de zeros. Quant val aquest vector r? Aquests exemples ilustren qu`e ´es el que fa MATLAB amb la comanda \quan li demanem que resolgui un sistema d’equacions lineals: - Si el sistema t´e soluci´o ´unica, ens la troba. - Si el sistema t´e soluci´o que dep´en de par`ametres, ens troba una soluci´o particular. Concretament, la m´es propera a l’origen. - Si el sistema no t´e soluci´o, el programa ens ofereix l’anomenada soluci´o aproximada per m´ınims quadrats, que estudiareu en el curs de Geometria, per`o que no ´es una soluci´o exacta del sistema. MATLAB sol fer aix`o sense avisar. L’´unic cas en que emet un av´ıs ´es si detecta que la matriu de coeficients At´e un rang m´es petit que el m`axim que li permet la seva mida. Aix´ı, per aprofitar tota la pot`encia de c`alcul de MATLAB per a resoldre sistemes lineals necessitarem saber fer dues coses: (i) Discutir el sistema lineal: t´e solucions? de quants par`ametres depenen? (ii) En el cas en qu`e la soluci´o no ´es ´unica sin´o que dep´en de mpar`ametres, trobar uns vectors directors u1, . . . , umde manera que la soluci´o general s’obtingui multiplicant-los per qualsevol valor del seu par`ametre, com en l’exemple 2. La discussi´o d’un sistema amb MATLAB es pot fer seguint el teorema de Rouch´e-Frobenius al peu de la lletra. Ho ilustrarem amb el Exemple 4: Discutim el sistema lineal    2x−y−z= 3 x+ 2y+z= 4 4x+ 3y+z= 10 Entrem-lo en MATLAB % matriu de coeficients A=[2 -1 -1; 1 2 1; 4 3 1] % terme independent b=[3;4;10] El teorema de Rouch´e-Frobenius ens diu que si trobem els rangs de la matriu de coeficients A i de la matriu ampliada (A|b) podrem dedu¨ır quantes solucions t´e el sistema. Crear la matriu ampliada ´es senzill % matriu ampliada del sistema Ab=[A b] I la comanda rank de Matlab ens troba el rang de la matriu % discussio del sistema pel teorema RF rank(A) rank(Ab) 16 CAP´ ITOL 3. SISTEMES LINEALS, I AN ` ALISI DEL GAS NATURAL. Podem comprovar a m`a que les seves solucions depenen de 2 par`ametres, i formen tot un pla en R3: la soluci´o general ´es  x y z  = 3 0 0  +y 1 1 0  +z −2 0 1  = 3 0 0  + 1−2 1 0 0 1  ·y z, on (3,0,0) ´es una soluci´o particular del sistema (un punt del pla), y, z poden tenir qualsevol valor real, i els dos vectors (1,1,0) i (−2,0,1) s´on vectors directors del pla de les solucions. Resolem el sistema amb Matlab: % matriu de coeficients A=[1 -1 2] % terme independent b=3 % resolem el sistema (anomenem xyz al vector d’incognites) xyz=A\b La resposta que d´ona el programa ´es xyz = 0 0 1.5000 MATLAB ens ha tornat una soluci´o ´unica! Podem comprovar f`acilment que ´es una soluci´o correcta, per`o el programa no ens avisa de que n’hi han infinites m´es. Exemple 3: Considerem el sistema lineal    2x−y= 3 x+ 2y= 4 x+y= 5 Aquest sistema no t´e soluci´o (t´e les dues equacions de l’exemple 1, seguides d’una tercera que la soluci´o de les dues primeres no compleix). Resolem-lo amb Matlab: A=[2 -1 ; 1 2; 1 1] % terme independent b=[3;4;5] % resolem el sistema (anomenem xy al vector d’incognites) xy=A\b La resposta que d´ona el programa ´es xy = 2.2857 1.2857 MATLAB ens ha tornat una soluci´o? El programa sap perfectament que no ho ´es. Si fos una soluci´o exacta es compliria que l’operaci´o matricial % prova de la solucio r=A*xy-b 17 hauria de ser un vector de zeros. Quant val aquest vector r? Aquests exemples ilustren qu`e ´es el que fa MATLAB amb la comanda \quan li demanem que resolgui un sistema d’equacions lineals: - Si el sistema t´e soluci´o ´unica, ens la troba. - Si el sistema t´e soluci´o que dep´en de par`ametres, ens troba una soluci´o particular. Concretament, la m´es propera a l’origen. - Si el sistema no t´e soluci´o, el programa ens ofereix l’anomenada soluci´o aproximada per m´ınims quadrats, que estudiareu en el curs de Geometria, per`o que no ´es una soluci´o exacta del sistema. MATLAB sol fer aix`o sense avisar. L’´unic cas en que emet un av´ıs ´es si detecta que la matriu de coeficients At´e un rang m´es petit que el m`axim que li permet la seva mida. Aix´ı, per aprofitar tota la pot`encia de c`alcul de MATLAB per a resoldre sistemes lineals necessitarem saber fer dues coses: (i) Discutir el sistema lineal: t´e solucions? de quants par`ametres depenen? (ii) En el cas en qu`e la soluci´o no ´es ´unica sin´o que dep´en de mpar`ametres, trobar uns vectors directors u1, . . . , umde manera que la soluci´o general s’obtingui multiplicant-los per qualsevol valor del seu par`ametre, com en l’exemple 2. La discussi´o d’un sistema amb MATLAB es pot fer seguint el teorema de Rouch´e-Frobenius al peu de la lletra. Ho ilustrarem amb el Exemple 4: Discutim el sistema lineal    2x−y−z= 3 x+ 2y+z= 4 4x+ 3y+z= 10 Entrem-lo en MATLAB % matriu de coeficients A=[2 -1 -1; 1 2 1; 4 3 1] % terme independent b=[3;4;10] El teorema de Rouch´e-Frobenius ens diu que si trobem els rangs de la matriu de coeficients A i de la matriu ampliada (A|b) podrem dedu¨ır quantes solucions t´e el sistema. Crear la matriu ampliada ´es senzill % matriu ampliada del sistema Ab=[A b] I la comanda rank de Matlab ens troba el rang de la matriu % discussio del sistema pel teorema RF rank(A) rank(Ab) 18 CAP´ ITOL 3. SISTEMES LINEALS, I AN ` ALISI DEL GAS NATURAL. En aquest exemple concret, MATLAB ens informa de que la matriu At´e rang 2, i la matriu (A|b) t´e rang 3. Per tant el sistema no t´e solucions. El procediment de l’exemple 4 permet discutir amb MATLAB qualsevol sistema d’equacions lineals. Si el sistema t´e soluci´o ´unica, aquesta ´es la que ens retorna MATLAB executant x=A\b. Ara ens falta veure, en el cas en que la soluci´o no ´es ´unica, com trobar la soluci´o general. Per aix`o hem de recordar l’estructura d’aquesta soluci´o general. Quan tenim un sistema d’equacions lineals en forma matricial Ax =b on A´es una matriu de mida m×n,x´es un vector columna amb ninc`ognites x1,...xn, i b´es un vector columna amb mtermes independents de les equacions, en el que els rangs de Ai de (A|b) valen r < n, de manera que la soluci´o del sistema dep´en de n−rpar`ametres, si trobem una soluci´o particular del sistema x0aleshores totes les solucions s´on de la forma x=x0+z on x0´es la soluci´o particular (en l’exemple 2 seria un punt del pla), i z´es qualsevol vector del subespai director del sistema, que ´es el subespai vectorial format per les solucions del sistema lineal homogene¨ıtzat Az =   0 . . . 0    (homogene¨ıtzar un sistema lineal vol dir canviar el seu terme independent per un de zeros, sense tocar la matriu de coeficients). Quan parlem d’aplicacions lineals, direm que el supespai d’aquests vectors z´es el nucli de la matriu A(identificant la matriu amb l’aplicaci´o lineal que defineix en bases can`oniques). MATLAB disposa d’una comanda per a calcular una base del nucli d’una matriu A. Com el nucli en angl`es es diu nullspace, aquesta comanda s’anomena null. Mostrarem com usar-la per a trobar totes les solucions d’un sistema compatible indeterminat en el Exemple 5: Considerem el sistema lineal 2x−y+z−3t= 3 x+ 2y+ 2t= 4 Entrem el sistema en MATLAB % matriu de coeficients A=[2 -1 1 -3; 1 2 0 2] % terme independent b=[3;4] % matriu ampliada del sistema Ab=[A b] Comencem discutint el sistema % discussio del sistema pel teorema RF rank(A) rank(Ab) Totes dues matrius t´enen rang 2. El sistema t´e doncs soluci´o, i com cont´e 4 inc`ognites la soluci´o dep´en de 4-2=2 par`ametres. Comen¸carem trobant una solucio particular del sistema 19 % solucio particular xyzt0=A\b I ara una base del nucli de la matriu de coeficients del sistema: % nucli de A Z=null(A) La matriu Zt´e per columnes els vectors d’una base del nucli de A. La soluci´o general del sistema lineal es pot escriure com     x y z t    =    x0 y0 z0 t0    +Zλ1 λ2, on (x0, y0, z0, t0) ´es la soluci´o particular que hem trobat fent A\b,Z´es la matriu que cont´e la base del nucli de Acom a columnes, i (λ1, λ2) s´on els par`ametres (un per cada columna de Z) als que es pot donar qualsevol valor per a definir una soluci´o. Notem una difer`encia important entre les solucions generals que troba MATLAB i les que trobem a m`a (com a exercici, podeu aplicar el procediment de l’exemple 5 per a trobar les solucions de l’equaci´o de l’exemple 2): quan treballem a m`a, els par`ametres dels que dep´en la soluci´o del sistema s´on les inc`ognites que no hem pogut despejar al resoldre’l, i els vectors directors que formen la base del nucli de As’obtenen a partir de cadascuna d’aquestes inc`ognites. En canvi, MATLAB busca una base ortonormal del nucli de la matriu A, i la retorna com les columnes de Z, de manera que nomes podem entrar els par`ametres com un vector columna que multiplica per la dreta a Z. Que una base sigui ortonormal vol dir que els vectors que la formen s´on ortogonals 2 a 2 i de longitud 1. Al curs de Geometria treballareu amb aquestes bases, i us explicaran que tenen molts avantatges enfront de les bases que trobem a m`a a l’hora de treballar. Aquests avantatges s´on la ra´o per la que les comandes de MATLAB que retornen una base solen escollir una base ortonormal. Finalment, proposem un problema f´ısic al que es pot trobar una soluci´o aproximada resolent sistemes lineals, de manera que com menys marge d’error vulguem tenir m´es grans seran els sistemes lineals a resoldre. Aplicaci´o: an`alisi del gas natural El gas natural t´e tres components en proporci´o apreciable: el met`a, CH4, que ´es el component principal (70-97% en pes), l’et`a, C2H6, i el prop`a, C3H8. Hi solen haver traces d’altres subst`ancies, per`o en proporci´o ¡¡ 1%. Per a valorar el gas natural que produeix un pou, que transporta un vaixell, . . . cal trobar la proporci´o que cont´e de cadascun dels tres components principals. Un procediment ´es prendre una mostra de 100g de gas i cremar–lo. La combusti´o del gas natural s´on aquestes reaccions qu´ımiques met`a: CH4+ 2O2−→ CO2+ 2H2O et`a: C2H6+7 2O2−→ 2CO2+ 3H2O prop`a: C3H8+ 5O2−→ 3CO2+ 4H2O (3.1) Cremem 100g de gas natural d’una partida que volem analitzar i mesurem: - Consum de 386.4g de O2. 20 CAP´ ITOL 3. SISTEMES LINEALS, I AN ` ALISI DEL GAS NATURAL. - Producci´o de 276.6g de CO2. - Producci´o de 221.07g de H2O. Podem amb aquesta informaci´o trobar els continguts x, y, z (en g) de met`a, et`a, prop`a en el gas cremat? Resoluci´o fallida # 1: Tradu¨ım en equacions la informaci´o que tenim. En primer lloc tenim l’equaci´o de la massa: x+y+z= 100 (els altres component del gas natural s´on inapreciables pesant amb balan¸ca). Usem com a pesos at`omics 1 per l’hidr`ogen, 12 pel carboni, 16 per l’ox´ıgen. Aleshores els pesos moleculars de les subst`ancies que intervenen en la reacci´o (3.1) s´on CH4: 12+4=16, C2H6: 30, C3H8: 44, O2: 32, CO2: 44, H2O: 18. Podem pensar les reaccions (3.1) en unitats de mol`ecules, o en mols (1 mol=6.023 ·1023 mol`ecules, i t´e per pes en grams el pes molecular). Triem la segona opci´o. Aix´ı, si la mostra cont´e xg de met`a, yg d’et`a i zg de prop`a pensarem que t´e x/16 mols de met`a, y/30 mols d’et`a, z/44 mols de prop`a. El consum en mols de O2en la combusti´o ser`a 2x 16 +7 2 y 30 + 5 z 44 =386.4 32 (els termes de l’esquerra s´on els consums de O2de les tres reaccions de (3.1), el de la dreta el consum total). De la mateixa manera, la producci´o de CO2´es x 16 + 2 y 30 + 3 z 44 =276.6 44 i la producci´o d’aigua ´es 2x 16 + 3 y 30 + 4 z 44 =221.07 18 Exercici 1: entreu les equacions de la massa, del consum de O2, de la producci´o de CO2i d’aigua com a sistema lineal en Matlab. Busqueu els rangs per a discutir el sistema. Soluci´o: surt rang A= 2,rang (A|b) = 3. Sistema incompatible! Qu`e ha fallat? rang A < 3 (rang m´es petit del m`axim possible) vol dir que les nostres equacions no s´on independents. Dues d’elles, per exemple la de la massa i la de la producci´o d’aigua, generen les altres dues. D’on ve que rang (A|b) = 3, que fa que el sistema no tingui soluci´o? Ve de que s’han fet arrodoniments en el terme independent (els pesos dels gasos estan arrodonits a la d`ecima de g m´es propera, el de l’aigua a la cent`essima de g), i en els coeficients del sistema (els pesos at`omics han tingut petits arrodoniments). Aquests arrodoniments han fet que el sistema deixi de tenir soluci´o: com rang A= 2, els vectors de terme independent ˜ b∈R4que s´on imatge de A: ˜ b=A x y z   per alguna tria de (x, y, z) formen un subespai de dimensi´o 2 de R4, o sigui un pla. Al canviar el terme independent exacte ˜ bper un arrodoniment b, hem obtingut un terme independent que est`a, per poc, fora d’aquest pla. Aix`o fa que el sistema lineal no tingui soluci´o exacta. 21 Figura 3.1: El problema de l’arrodoniment del terme independent quan el rang del sistema es menor que el nombre d’equacions. S’arreglaria aquest problema obtenint pesos exactes? No. Com rang A= 2, el sistema nom´es permet despejar 2 inc`ognites. Encara que passi a tenir soluci´o, mai ens permetr`a trobar valors ´unics per x, y, z. El que necessitem no ´es pesar amb m´es precissi´o, sin´o noves equacions. Resoluci´o fallida #2: Mesurem el volum del gas. El gas natural t´e mol`ecules petites, i es comporta com un gas ideal. El volum que ocupa segueix la llei del gas ideal (o de Boyle, o de Boyle–Mariotte . . . ) PV =nRT on P= pressi´o, V= volum que ocupa el gas, n= nombre de mols de gas, R= constant dels gasos (universal, i coneguda), T= temperatura (en K). Aix´ı n=PV RT Per tant, si mesurem el volum que ocupa la mostra de 100g de gas natural, i a quina pressi´o i temperatura ´es troba, podem dedu¨ır el nombre de mols nque cont´e. Fent aquesta mesura en la nostra mostra obtenim n= 5.9955, i d’aqu´ı obtenim l’equaci´o del volum x 16 +y 30 +z 44 = 5.9955 Exercici 2: Afegiu l’equaci´o del volum al sistema de l’exercici 1, i discutiu el sistema resultant. Soluci´o: torna a sortir rang A= 2,rang (A|b) = 3. Que ha fallat? rang A= 2 vol dir que l’equaci´o del volum es pot dedu¨ır de les equacions pr`evies (per exemple, de la de la massa i de la de la producci´o d’aigua es continuen dedu¨ınt totes). Que la matriu ampliada surti de rang 3 ha estat causat pels arrodoniments en els pesos, com en l’intent anterior. Resoluci´o amb `exit # 3: Passem a controlar una altra magnitud f´ısica: l’energia produ¨ıda. Fem la combusti´o del gas a pressi´o constant (en obert), deixem circular el O2iCO2perque no ens han aportat res al sistema, i condensem tota l’aigua produ¨ıda per a facilitar la seva mesura. Al cremar en aquestes condicions, es sap que - 1 mol de met`a produeix 212.80 kcal, - 1 mol d’et`a produeix 372.82 kcal, - 1 mol de prop`a produeix 530.60 kcal. Mesurem l’energia produ¨ıda per la combusti´o de la nostra mostra de 100g de gas natural, i resulta ser 1322.3 kcal. Obtenim doncs una equaci´o de l’energia (en kcal): 212.8x 16 + 372.82 y 30 + 530.6z 44 = 1322.3 22 CAP´ ITOL 3. SISTEMES LINEALS, I AN ` ALISI DEL GAS NATURAL. Exercici 3: Formeu un sistema lineal amb nom´es les equacions de la massa, de la producci´o de l’aigua, i de l’energia. Discutiu amb Matlab el sistema resultant, i si t´e soluci´o ´unica, trobeu-la. Soluci´o: Ara s´ı! rang A= 3 = rang (A|b) =. Per tant el sistema t´e soluci´o ´unica, que val x= 91.8289g, y = 6.6250g, z = 1.5461g. Moralitat d’aquesta pr`actica: Que un sistema lineal (A|b) no tingui el rang m`axim possible vol dir que hi han relacions entre les seves equacions. En la vida real, aix`o vol dir relacions entre les magnituds que estudien. En el nostre exemple, les reaccions de combusti´o (3.1) fan que a partir del nombre total de mols d’hidrocarburs i de la producci´o d’aigua es puguin dedu¨ır el consum d’ox´ıgen, la producci´o de CO2, fins i tot el volum que ocupa el gas abans de cremar–lo . . . sense saber quants mols tenim de cada hidrocarbur espec´ıfic. Mirant nom´es el rang de la matriu Ade coeficients del sistema, hem dedu¨ıt que per a analitzar el gas natural necessitem nom´es controlar el volum de gas cremat, el calor que ha em`es, i l’aigua produ¨ıda en el condensador. No fa falta mesurar el consum de O2ni la producci´o de CO2, que s´on bastant m´es costoses de comptar. Cap´ıtol 4 Tractament d’imatges. MATLAB ´es una eina molt potent per al tractament d’imatges. Te menys comandes d’alt nivell per a tractar-les que els programes especialitzats en aquestes tasques, per`o la filosofia de tot el proc`es d’imatges modern ´es digitalitzar-les, ´es a dir convertir-les en matrius num`eriques, i programant Matlab podem executar qualsevol operaci´o amb matrius. En aquesta pr`actica farem alguns tractaments t´ıpics d’imatge, que resulta que s´on aplicacions lineals entre espais de matrius. El men´u de tasques per aquesta pr`actica ´es: (i) Importar una imatge en color a Matlab. (ii) Retallar una imatge. (iii) Canvis de color. (iv) Reescalar la imatge a la baixa. (v) Fer detecci´o de vores (verticals u horizontals) en la imatge. (vi) Guardar la imatge editada en Matlab com un fitxer d’imatge ordinari. (i) Importaci´o d’una imatge a Matlab: Hem de tenir la imatge continguda en un fitxer en el directori de treball de Matlab, i guardar el nom d’aquest fitxer en una cadena de texte per a que la comanda d’importaci´o d’imatges a Matlab la llegeixi. En el nostre exemple infile=’vaca_estirada.jpg’; A0=imread(infile); A0 ´es una taula 3-dimensional de nombres (un tensor 3d). Es pot pensar com que te 3 pisos A0(:,:,1),A0(:,:,2),A0(:,:,3), corresponent respectivament als colors Vermell (R), Verd (G), Blau (B), i que cada pis es una matriu que cont´e la intensitat del color en els p´ıxels de la imatge. Els coeficients de A0 s´on nombres enters, usualment van de 0 a 255 (enters de 8 bits sense signe). Aixi cada coeficient ocupa nomes 1 Byte de memoria. En contrast: Matlab usa 8 Bytes de memoria per a guardar cada nombre real (usa el format doble precissi´o). Com a verificaci´o de que la imatge s’ha llegit correctament, fem que Matlab l’ensenyi com a figura 1 (numerarem les figures per a tenir-ne v`aries obertes simult`aneament): figure(1) imshow(A0) 23 24 CAP´ ITOL 4. TRACTAMENT D’IMATGES. Figura 4.1: Imatge de mostra per aquesta pr`actica: vaca estirada.jpg. (ii) Retallar una imatge: Si volem quedar-nos nom´es amb la cara de la vaca: en la imatge inicial podem esbrinar un rectangle que la contingui amb l’eina ’Data Cursor’ de la barra de sobre de la imatge. Figura 4.2: L’eina Data Cursor en una figura de Matlab. Un cop hem vist que la cara esta entre les files 20 i 200 i les columnes 10 i 210, podem quedar-nos aquesta part de la taula num`erica per a tenir una imatge amb la cara de la vaca. Aret=A0(20:200,10:210,:); figure(2) imshow(Aret) (iii) Conversi´o de color a escala de grisos: Per a operar amb les matrius de color dels pixels, cal passar les taules de nombres de la imatge de nombres enters a reals. De pas canviem l’escala d’intensitats de [0,255] a [0,1]: A=double(A0)/255; Despr`es, separem la taula 3-dimensional Aen les matrius dels 3 colors RGB de la imatge 25 AR=A(:,:,1); AG=A(:,:,2); AB=A(:,:,3); En aquest punt tenim la imatge codificada dins de Matlab com 3 matrius de nombres reals, que son les intensitats del color respectiu. Amb aquestes podrem operar. Per exemple, passar la imatge a escala de grisos ´es canviar les 3 matrius d’intensitats de color R,G,B a una ´unica matriu d’intensitat de gris, que ´es el promig de les matrius originals: Agr=(AR+AG+AB)/3; Verificaci´o: que Matlab ensenyi la imatge en escala de grisos. Primer canviem d’intensitats reals en [0,1] a intensitats enteres en [0,255] (arrodonint a l’enter m´es proper). A0gr=uint8(Agr*255); I ara ja podem demanar a Matlab la imatge: figure(2) imshow(A0gr) Figura 4.3: La imatge de mostra, passada de color a escala de grisos. (iv) Reescalat a la baixa d’una imatge: El cas m´es senzill de reescalat ´es quan volem fer que la imatge tingui menys resoluci´o (perque nom´es volem mostrar-la a tamany petit . . . ). En la nostra imatge d’exemple, podem canviar cada bloc de 4 ×4 p´ıxels per un ´unic pixel, de manera que passem d’una imatge de mida 356 ×476 a una de mida 89 ×119. Quina intensitat de color s’ha de posar al p´ıxel que reempla¸ca a un bloc 4 ×4? La millor tria ´es el promig de les intensitats dels 16 p´ıxels del bloc 4×4. Programem aix`o: trobem la mida de la imatge a reescalar, i fem un bucle doble per a calcular els pixels de la imatge reescalada. [m,n]=size(Agr); f0=1; % fila inicial del bloc for k=1:m/4, c0=1; % columna inicial del bloc for l=1:n/4, bloc=Agr(f0:f0+3,c0:c0+3); prombloc=sum(bloc(:))/16; Ared(k,l)=prombloc; c0=c0+4; end; f0=f0+4; end; 32CAP´ ITOL 5. XARXA NEURONAL DE TIPUS MULTILAYERPERCEPTRON (MLP) PER A APROXIMAR LES ARRELS D’UN POLINOMI Figura 5.3: Flux de treball (pipeline) per a la nostra xarxa MLP. a2x2+a1x+a0. Permetrem que els coeficients del polinomi siguin complexes, de manera que el polinomi p(x) es convertir`a en un vector v0= (Re (a4),Re (a3),...,Re (a0),Im (a4),...,Im (a0)) , que posarem en columna en Matlab. 1) Aplicarem una primera capa FC per a transformar el vector v0en un vector v1∈R50 donat per v1=A1v0+b1. Per a aix`o cal que A1sigui una matriu de mida 50 ×10, i b1un vector de R50. Entre A1ib1 contenen 500+50=550 coeficients, que s´on els par`ametres que necessita tenir definits aquesta capa. 1R) Aplicarem una capa ReLU al vector v1, per a obtenir un vector activat vr1=ReLU(v1), de la mateixa dimensi´o. Aquesta capa no necessita par`ametres. 2) Aplicarem una segona capa FC per a transformar el vector vr1en un vector v2∈R25 donat per v2=A2vr1+b2. Per a aix`o cal que A2sigui una matriu de mida 25 ×50, i b2un vector de R25. Entre A2ib2 contenen 1250+25=1275 coeficients, que s´on els par`ametres que necessita tenir definits aquesta capa. 2R) Aplicarem una capa ReLU al vector v2, per a obtenir un vector activat vr2=ReLU(v2). Aquesta capa no necessita par`ametres. 3) Finalment, aplicarem una tercera capa FC per a transformar el vector vr2en un vector v3∈R10 donat per v3=A3vr2+b3. Per a aix`o cal que A3sigui una matriu de mida 10 ×25, i b3un vector de R10. Entre A3ib3 contenen 250+10=260 coeficients, que s´on els par`ametres que necessita tenir definits aquesta capa. El vector v3que retorna la tecera capa FC de la xarxa ´es la seva resposta. L’entenem com un vector de la forma (Re (z1),Re (z2),...,Re (z5),Im (z1),...,Im (z5)), on z1, z2, . . . , z5s´on els valors que proposa la xarxa per les arrels del polinomi p(x). Validaci´o de la xarxa: Per a poder entrenar una xarxa neuronal amb `exit necessitem quantificar de manera fina l’error que es comet en cada resposta que ens d´ona, perqu`e l’entrenament consisteix en moure lentament els valors dels par`ametres en una direcci´o que provoqui disminucions, tamb´e lentes, de l’error en la resposta. 33 En la nostra xarxa podem aplicar la mesura de l’error m´es senzilla de totes: si la xarxa ha donat un vector v3∈R10 com a posici´o estimada de les arrels de p(x), i les posicions exactes s1, . . . , s5formen un vector s= (Re (s1),Re (s2),...,Re (s5),Im (s1),...,Im (s5)), considerarem que l’error de la resposta ´es la funci´o de p`erdua (o loss function) L(v3) = kv3−sk2, ´es a dir, la norma al quadrat del vector difer`encia v3−s(el producte escalar de v3−sper si mateix). Aquesta mesura de l’error t´e un problema en el cas de la nostra pregunta de trobar les arrels d’un polinomi p(x): encara que el vector de resposta v3contingui els valors correctes de les arrels, si hi apareixen en un ordre diferent de com apareixen a la soluci´o strobarem un error important, quan en realitat hauria de ser 0! Resoldrem aquesta dificultat ordenant les llistes d’arrels del polinomi (tant les pronosticades per la xarxa com les de la soluci´o correcta quan estiguem entrenant-la) abans de calcular la funci´o d’error L. Els nombres complexes formen un pla, i aqu´ı no hi ha cap ordre complet que no tingui algun defecte. Per exemple, si ordenem les arrels del polinomi ordenant la seva part real ens trobarem amb un empat entre les arrels complexes conjugades d’un polinomi real. Si ordenem les arrels del polinomi ordenant la seva part imagin`aria ens trobarem amb l’empat entre les arrels reals del polinomi. Si ordenem les arrels lexicogr`aficament, tenint en compte en primer lloc la part real i usant la part imagin`aria per trencar empats, ens trobarem amb que si per un error d’arrodoniment del c`alcul la part real d’una arrel surt 10−14 enlloc de 0, un error tan petit ja canvia l’ordre de les arrels . . . Per aquestes raons usarem un ordre incomplet en el pla complex, que sigui r`apid de calcular i no provoqui empats entre arrels reals o entre arrels conjugades quan aquestes apareguin: pensarem els nombres complexes com el pla R2i ordenarem cada nombre complex z=a+ib = (a, b) segons el valor del producte escalar π(z) = (1,−0.1) ·(a, b) = a−0.1b. Aquest producte escalar equival a projectar ortogonalment el vector (a, b) sobre la recta amb vector director u= (1,−0.1). Figura 5.4: Projecci´o ortogonal sobre una recta orientada per a ordenar els punts del pla. 3. EL CODI: CONSTRUCCI´ O I PROVA DE LA XARXA MLP AMB PAR` AMETRES DONATS El primer pas que farem ser`a preparar exemples de verificaci´o: una familia de polinomis de grau 5 a coeficients complexes, aleatoris, pels que coneixerem amb gran precissi´o les seves arrels. La manera m´es senzilla de fer aix`o es generar les arrels com punts aleatoris en un rectangle del pla complex, i obtenir el polinomi p(x) que t´e unes arrels z1, . . . , z5prefixades fent les multiplicacions (x−z1)·(x−z2). . . (x−z5): 34CAP´ ITOL 5. XARXA NEURONAL DE TIPUS MULTILAYERPERCEPTRON (MLP) PER A APROXIMAR LES ARRELS D’UN POLINOMI nex=10000; % nombre d’exemples: DEIXAR A 4 O 8 PER VEURE COM FUNCIONA EL CODI % finestra on posarem els zeros xlim=[-2 2]; ylim=[-2 2]; % creem la llista de polinomis (files d’una matriu) i de les seves arrels rng(20241115) % iniciem el generador de nombres aleatoris de Matlab % part real de les arrels: matriu 5 x nex de nombres aleatoris entre els % limits de x rearrels=xlim(1)+rand(5,nex)*(xlim(2)-xlim(1)); % idem per la part imaginaria de les arrels imarrels=ylim(1)+rand(5,nex)*(ylim(2)-ylim(1)); arrels=rearrels+i*imarrels; % cada columna son les arrels d’un exemple % polsex te un polinomi d’exemple en cada fila (incloent el coeficient % inicial que es 1): % de entrada declarem cada polinomi x-arrel(1), i el multipliquem fila a % fila per x-arrel(2),...,x-arrel(5) for k=1:nex, polex=[1 -arrels(1,k)]; for l=2:5, polex=conv(polex,[1 -arrels(l,k)]); end; polsex(k,:)=polex; end; Ja tenim una llista de polinomis de grau 5 (les files de la matriu polsex) amb les arrels conegudes (guardades com a columnes de la matriu arrels). ´ Es el moment de muntar la xarxa neuronal. El primer pas ´es llegir els par`ametres que alg´u ha trobat entrenant la xarxa i ens passa: load parametres_MLP_20241115.mat Apareixen en la mem`oria de Matlab 4 variables: params=vector amb tots els par`ametres, f=nombre de files en cada capa, c=nombre de columnes en cada capa, parcapa=matriu que en cada fila ens diu quines posicions en params ocupen els par`ametres de la seva capa del MLP. A partir de les dades llegides hem de muntar les matrius i vectors de les capes FC. El conveni de Matlab per guardar matrius ´es que la matriu es guarda com una llista dels seus coeficients, llegits columna a columna d’esquerra a dreta, i cada columna llegida de dalt a baix. Per exemple la matriu 11 12 13 21 22 23 es guarda com una llista 11,21,12,22,13,23 . La comanda reshape de Matlab ens permet reconstruir la matriu a partir d’una llista dels seus coeficients, seguint aquest conveni al peu de la lletra. % coeficients de la matriu A1 en la llista de parametres: indexos inicial i % final ind0A1=parcapa(1,1); indfA1=parcapa(1,1)+f(1)*c(1)-1; % organitza els coeficients en una matriu de mida f(1) x c(1) A1=reshape(params(ind0A1:indfA1),f(1),c(1)); % muntem el vector b1 amb els darrers parametres d’aquesta capa 35 ind0b1=parcapa(1,2)-f(1)+1; indfb1=parcapa(1,2); b1=reshape(params(ind0b1:indfb1),f(1),1); % Repetim el procediment identicament per A2,b2,A3,b3 ind0A2=parcapa(2,1); indfA2=parcapa(2,1)+f(2)*c(2)-1; A2=reshape(params(ind0A2:indfA2),f(2),c(2)); ind0b2=parcapa(2,2)-f(2)+1; indfb2=parcapa(2,2); b2=reshape(params(ind0b2:indfb2),f(2),1); ind0A3=parcapa(3,1); indfA3=parcapa(3,1)+f(3)*c(3)-1; A3=reshape(params(ind0A3:indfA3),f(3),c(3)); ind0b3=parcapa(3,2)-f(3)+1; indfb3=parcapa(3,2); b3=reshape(params(ind0b3:indfb3),f(3),1); A continuaci´o, farem un experiment per a comprovar com s’implementa una capa ReLU usant un vector prou petit per a que el vegem c`omodament per pantalla vtest=[3;-1;0] vtestampl=[vtest zeros(size(vtest))] % afegim columna de zeros a la dreta vtestr=max(vtestampl,[],2) % el maxim de cada fila es l’activacio ReLU Ja estem a punt per a aplicar la xarxa MLP a tots els polinomis continguts en les files de la matriu d’exemples polsex. Per`o abans de fer-ho prepararem un darrer detall: l’ordenaci´o de les arrels del polinomi. Com aquesta ordenaci´o necessita algunes l´ınees de codi, i l’haurem d’executar aquest proc`es almenys 2 cops en aquesta pr`actica, el programarem com una funci´o de Matlab. Una funci´o ´es un programa auxiliar, que ´es cridat per un programa principal per a que faci una feina particular i li retorni el resultat. La funci´o queda definida per: - les dades d’entrada que ha de rebre del programa que la crida, - el proc`es que s’ha de fer amb aquestes dades, - les dades de sortida que s’han de retornar al acabar l’execuci´o de la funci´o al programa que l’ha cridat. Nosaltres crearem una funci´o amb nom ordenaR2 que rebi com a dades d’entrada un vector columna v∈R10 (les 5 arrels a ordenar, pensades com a vector v= (Re (z1),Re (z2),...,Re (z5),Im (z1),...,Im (z5)) en columna), i un vector columna u∈R2que sigui director de la recta en la que s’han de projectar les 5 arrels. La funci´o calcular`a aquestes projeccions, reordenar`a les arrels segons la projecci´o usant la comanda sort de Matlab, i retornar`a un vector vord amb les 5 arrels ordenades, per`o escrites en el mateix format que s’ha usat pel vector vd’entrada. La manera m´es senzilla d’usar funcions en Matlab es guardar cada funci´o en un fitxer que es digui exactament com ella, en el nostre cas ordenaR2.m, i que estigui en el directori de treball de Matlab. El contingut del fitxer de la nostra funci´o ser`a 36CAP´ ITOL 5. XARXA NEURONAL DE TIPUS MULTILAYERPERCEPTRON (MLP) PER A APROXIMAR LES ARRELS D’UN POLINOMI function vord=ordenaR2(v,u) % nombre de punts del pla es la mitat del de components de v npunts=length(v)/2; % posa les components de v com npunts vectors columna de R^2, en una matriu vpla=[v(1:npunts).’; v(npunts+1:end).’]; % calcula les projeccions ortogonals de les columnes de la matriu sobre la % recta amb vector director u projs=(u.’)*vpla; % busca l’ordre de les projeccions [~,ordre]=sort(projs); % reordena les columnes de vpla segons el nou ordre trobat vpla=vpla(:,ordre); % a partir de la matriu vpla reordenada torna a formar un vector columna % amb totes les components x, i despres totes les components y vord=[vpla(1,:).’ ; vpla(2,:).’]; Noteu la primera linea: ´es la declaraci´o de la funci´o. Avisa a Matlab de que aquest fitxer cont´e una funci´o i li indica quin ´es el nom de la funci´o i qui s´on les variables d’entrada i de sortida. Al guardar la funci´o en el seu propi fitxer .m, no cal escriure cap comanda tipus return per a indicar que la funci´o s’ha acabat. Quan escribim una funci´o de Matlab que hem d’usar amb dades nombroses, complicades . . .´es prudent provar primer la funci´o en dades senzilles on puguem veure a ull que funciona correctament. Fem aquesta prova amb 4 punts de R2: vtest2=[1 1 2 2 1 -11 2 -3].’ % 4 punts de R^2 com un vector columna de 8 components vtest2ord=ordenaR2(vtest2,[1;-0.1]) Un cop hem comprovat que la funci´o per a ordenar les llistes d’arrels funciona correctament, ja podem fer correr la nostra xarxa MLP i mesurar l’error que comet en cada exemple. Com la nostra xarxa ´es molt simple (nom´es 2085 par`ametres), no aconseguirem una aproximaci´o gaire fina de les 5 arrels de cada polinomi on la provem. Per a veure si hem obtingut m´es informaci´o que inventant 5 punts a l’atzar en el rectangle on hem imposat que estiguin les arrels, farem una comprovaci´o estad´ıstica molt senzilla: - triarem un marge d’error molt tolerant Lmax = 10, - i comptarem per quants polinomis de la nostra llista de verificaci´o l’error global entre les posicions ordenades de les 5 arrels estimades per la xarxa i les 5 arrels correctes ´es m´es petit que Lmax. Ara que ja hem decidit tot el que farem amb cada polinomi de verificaci´o (aplicar-li la xarxa capa a capa, ordenar les arrels que li ha estimat, ordenar les arrels correctes del polinomi, calcular l’error com`es per la xarxa i mirar si cau per sota del marge perm`es), ja podem fer el bucle principal que aplicar`a la xarxa i verificar`a l’error que comet: % marge d’error per sota del que donarem per bo el calcul Lmax=10; % ull! es un marge d’error bastant generos esencert=zeros(1,nex); % aqui guardarem si l’error de la xarxa es admissible o no % aplicacio de la xarxa a cada polinomi de la llista d’exemples, amb % verificacio de l’error comes 37 for k=1:nex, polex=polsex(k,:); % un polinomi de la llista, en fila v0=[real(polex(2:end)) imag(polex(2:end))]; % els coeficients, sense el x^5 inicial v0=v0.’; % posem els 5 coeficients com un vector columna de R^10 % apliquem capa FC1 a la dada v1=A1*v0+b1; % apliquem ReLU al resultat v1ampl=[v1 zeros(size(v1))]; v1r=max(v1ampl,[],2); % apliquem capa FC2 a la dada v2=A2*v1r+b2; % apliquem ReLU al resultat v2ampl=[v2 zeros(size(v2))]; v2r=max(v2ampl,[],2); % apliquem capa FC3 a la dada v3=A3*v2r+b3; % v3 son les arrels que ha pronosticat la xarxa, en format vector de R^10 % Ordenem aquestes arrels v3ord=ordenaR2(v3,[1;-0.1]); % i les guardem, per si serveixen d’alguna cosa arrxarxa(:,k)=v3ord(1:5)+i*v3ord(6:10); % carreguem les arrels correctes, com a vector de R^10 sol=[real(arrels(:,k)); imag(arrels(:,k))]; % les ordenem solord=ordenaR2(sol,[1;-0.1]); % mesurem l’error global comes per la xarxa i el guardem en la llista % d’errors err=norm(v3ord-solord); L(k)=err*err; % mirem si l’error comes per la xarxa era admissible if L(k)<Lmax esencert(k)=1; end; end; Ara comptem en quants polinomis l’estimaci´o de les arrels donada per la xarxa es troba dins del marge d’error que hem adm`es. El vector esencert conte un 1 per cada polinomi on aix`o ha succe¨ıt, i un 0 per cada polinomi on no ha succe¨ıt. Per tant nomes cal sumar les seves components. nencerts=sum(esencert) Pel marge d’error L < 10 la resposta de la xarxa ha funcionat en un 92.94% dels casos. Sembla un marge molt bo, per`o hem de tenir en compte que aquest marge d’error ´es molt ampli. Si dibuixem les 5 arrels estimades per la xarxa (posant ∗en negre) i les 5 arrels correctes (posant o en vermell) per un polinomi dels encertats veurem que el marge d’error era ampli. 38CAP´ ITOL 5. XARXA NEURONAL DE TIPUS MULTILAYERPERCEPTRON (MLP) PER A APROXIMAR LES ARRELS D’UN POLINOMI k=1; figure(1) plot(real(arrxarxa(:,k)),imag(arrxarxa(:,k)),’k*’); hold on plot(real(arrels(:,k)),imag(arrels(:,k)),’ro’); axis([xlim(1) xlim(2) ylim(1) ylim(2)]) axis equal hold off Figura 5.5: Posicions de les arrels correctes (vermell) i pronosticades per la xarxa (negre) pel primer polinomi de la llista de validaci´o. Com ha quedat clar que el marge d’error que permetem ´es ampli, hem de fer una verificaci´o que ´es recomanable en qualsevol cas en que una xarxa neuronal nom´es tingui `exit parcial: Exercici final: la xarxa del mono. (i) Escriviu codi de Matlab per a generar conjunts aleatoris de 5 nombres complexes del mateix rectangle d’on hem seleccionat les arrels pels nostres polinomis de validaci´o. Genereu tants d’aquests conjunts com polinomis de validaci´o, i poseu-los com a columnes d’una matriu. Declarem que cadascuna d’aquestes columnes ´es una nova estimaci´o de la posici´o de les arrels del polinomi corresponent en la llista de validaci´o. Aquest procediment per a estimar les arrels ´es comparable a tenir un mono llen¸cant dards a una diana que sigui el pla complex, i tenint en compte nom´es els dards que es clavin en la regi´o on preveiem les arrels. (ii) Calculeu l’error que ha com`es la xarxa del mono en estimar les arrels dels polinomis de validaci´o que hem emprat per a la nostra xarxa MLP, usant el mateix procediment d’ordenar i mirar la norma al quadrat. Trobeu quants encerts ha fet la xarxa del mono amb el marge d’error que li hem deixat a la nostra xarxa MLP, compareu els encerts de les dues xarxes i discutiu si la xarxa MLP ha aportat informaci´o sobre la posici´o de les arrels dels polinomis. Indicaci´o: Podeu resoldre tot l’exercici copiant i enganxant trossos dels c`alculs que hem fet fins ara en la pr`actica. 39 Figura 5.6: El model biol`ogic que va inspirar la validaci´o definitiva per a xarxes neuronals de precissi´o limitada (imatge de valueofstocks.com).