Apostila Proc Iml

I n trtrod oduu ção ao SA SAS/I ML p araraa Álg Álgebra ebra de de Matrizes Escola Superior de Agric ultu ra “Luiz de Queiroz” Departamento de Ciências Exatas - ESALQ/USP Fevereiro de 2002 (versão 3,14)

Antonio Augusto Franco Garcia Anto nio Fra Francisco ncisco Iem m a

I n t rod In oduu çã çãoo ao a o SA S AS /I/IM ML p paa ra Álge lgebr braa d e Ma M a t ri zes Esse texto de stinastina-se se aos alun os do curso de verão d e Álg Álgebra ebra de Matrizes do Departamento de Ciências Exatas da ESALQ/USP, visando auxiliar na resolução dos exercícios propostos e fornecendo ao mesmo tempo subsídios para a utiliza utilizaçã çãoo de u m dos program program as co com m putaci putacionais onais m ai aiss em prega pregados dos n a Estatística: o SAS (Statistical Analysis System). Optou-se por apresentar as possíveis aplicações do módulo IML mantendo a mesma seqüência do texto usado em aula, de modo a complementá-lo. Leitores Lei tores int eressados em m aiores detalhes p oderão consu ltar o livro livro SSAS AS/I /IML ML Software: Usage and Reference, Version 6, First Edition, Cary, NC: SAS Institute Inc., 1989. 501pp. (6th printing, December 1995), além de outras pu blicaç blicações ões do referido instituto. • Por ter finali finalidades dades didátic didáticas, os programas programas aqui cita citados dos foram foram elaborados elaborados de forma a facilitar o entendimento do assunto, não visando exclusivam excl usivam en te aplicações aplicações práticas práticas ou an áli álise se de d ados. Sup Sup õeõe-se se contu do que esta forma de abordagem propiciará aos leitores com interesse aplicado os subsídios necessários para o desenvolvimento de seus próprios programas. • Críti Crítica cass e sugestões, sugestões, bem com o correç correção ão de eventu eventu ai aiss erros no m ateri aterial al,, serão bem re cebi cebidas. das.

C o n si de d e r a çõ ç õ es so b r e o S A S / I M L

1

Co n s i d e raç õ es s o b re o SAS/IM L O softw are SAS/ SAS/IML IML ((Interactiv Interactivee M atrix Language) Language) é um a linguagem linguagem de p rogra rogram m ação m uito flexível flexível,, tend o com o prin cipais cipais vantagens a in terativi teratividade dade com o usu ário, a possibilidade possibilidade de trabalhar com com m atriz atrizes es de form a sim sim ples, além da p ossibi ossibililidade dade do u so de fun ções préprogramadas prog ramadas e outras procedures do próprio SA SAS. Os resultados dos coman dos podem ser vistos vistos ime diatam ent e ou arm azenados p ara posterior execução, execução, a cri critério tério do usu ário. Nesse texto, por razões óbvias, será dada ê n fase as ope rações m atric atriciai iais. s. No en tan to, o SAS/ SAS/IML IML é uma linguagem de programação, permitindo que o usuário escreva seus próprios programas. Operações como inversão de matrizes, obtenção de autovalores e autovetores, etc, podem ser realizadas usando comandos pré-programados, facilitando o trabalho e evitando o uso de alg algoritm oritm os sofisti sofisticados cados que u til tiliz izam am vári várias as linhas de p rogra rogram m ação. O fato dos algoritmos algoritmos em preg pregados ados serem am plamen te testados e discutidos discutidos por usuários usuários do m un do todo perm ite que os resultados obtidos sejam muito confiáveis. Além disso, o SAS/IML faz uso muito eficiente dos recursos do computador, facilitando o trabalho com matrizes de grandes di mensões. As sessões do IML iniciam-se com o comando "proc iml;" e finalizam-se com o comando "quit;". Uma matriz é definida entre chaves ({}), sendo que as vírgulas indicam as m udan ças de linh linh as. Um p rimeiro programa: programa: p r o c i ml ; a ={ ={ 3 - 1 2 , 2 - 2 3, 4 1 - 4} ; x =i =i n v ( a ) ; pr i n t x ; qui t ;

X 0. 333 333333 3333 3 - 0. 13333 133333 3 0. 066 066666 6667 7 1. 3333 333333 333 3 - 1. 3333 333333 33 - 0. 3333 333333 33 0. 6666 666666 667 7 - 0. 4666 466667 67 - 0. 2666 266667 67

Esse programa permite a obtenção da matriz X, que é a inversa de A. O comando PRINT solicita solic ita a visualização visualização da m atriz X. X. Vale Vale ressaltar qu e, ao usar o com an do Q UIT, UIT, a sessã o do IML é en cerrada e as m atriz atrizes es e vari variávei áveiss criadas criadas são apagadas da m em ória RA RAM.

I nntt r o d u ç ã o a o S A S / I M L pa p a r a Á l g e b r a ddee M at a t r i zzee s

1

D eeff i n iinndd o m at r i zes

2

Definindo m atrizes O objeto fundam ental sobre o qual operam os com andos do IML é um a m atriz com linh as e colunas, podendo ser n um érica ou alfanu m érica. O núm ero m áxim o de caracteres que um elemento de uma matriz alfanumérica pode apresentar é 32676, sendo que todos os elemen tos serão armazenados com o m esm o nú m ero de caracteres que o m aior deles. Os nomes atribuídos as matrizes devem ser válidos no ambiente SAS, seguindo suas regras: possuir de 1 a 8 caracteres sem espaços e começar com letras, aceitando números a partir do 2º caractere (versões 6.12 e an teriores). Em versões mais recentes, como a 8.0 e 8.2, os nomes das variáveis podem exceder oito caracteres. É possível criar m atrizes com valores perdidos, sendo qu e os m esm os devem ser indicados por um ponto ("."). Contu do, nem todos os comand os funcionam com m atrizes com p ontos. Uma vez definido um nom e para a matriz, ele pode ser usado diretamen te ou su bstituído em ações futuras, a critério do usuário, sendo possível inclusive substituir apenas determinados elementos, conforme será visto a seguir. Ao se referir a uma dada matriz, pode ser simp lesm ente u tilizado seu n om e, ou ainda esp ecificand o seus elemen tos. Matrizes com m últiplos elem en tos são especificadas com o u so de chaves ({ }), sen do que as linhas são separadas por vírgulas. Evidentemente, todas as linhas devem conter o m esmo núm er o de elemen t os. Os valores a serem adicionados podem ser: • Núm eros reais inteiros ou com decim ais (notação científica, com o 1E-8, é perm itida). • Caracteres alfanum éricos ("strings"), que podem estar entre aspas sim ples (´ ´ ) ou duplas (" "), ou até mesmo sem elas, sendo que nesse caso todos os valores serão armazenados com letras m aiúsculas. • Núm er os ent r e colchet es ( [ ] ) indi cam o nú m er o de vezes que um dado valor de r epet e, facilitando a digitação em m uitas situações. Matrizes também podem ser criadas a partir do uso de funções pré-programadas, do comando CALL ou de um programa específico. A seguir, são apresentados exemplos de m atrizes com um ún ico elem ento (escalares): pr oc a=32; b=. ; c=' oi d=' ol e=" c i pr i nt

i ml ;

' ; á pes soal ' ; ao" ; a b c d e;

A 32

2

B C . oi

D ol á pes soal

E c i ao

I nt r oduçã oduçãoo ao S A S / I M L par zes pa r a Á lgebr a de M aatt r iizes

D ef i nniinn d o m at r i ze s

Já as matrizes com m últiplos elem en tos são escritas da seguin te m aneira: pr oc i ml ; x={ 10 20 30 40} ; y ={ 10, 20, 30, 40} ; z ={ 1 2, 3 4, 5 6} ; pr i nt x y z;

X 10

Y 20 20 30 40

Z 30 3 5

40 4 6

10

1

2

Ou ainda, para caracteres alfanum éricos: este pr oc aa={ bb={ ' pr i nt

i ml ; abcd Ef gh} ; abc d' ' Ef gh' } ; aa bb;

AA BB ABCD EFGH ab cd Ef gh

Pelo fato da matriz AA não apresentar seus elementos entre aspas, eles automaticamente são su bstitu ídos po r letras m aiúscu las. Já para BB, foi resp eitad a a form a original. Caso seja necessário definir elementos com espaços entre eles, deve-se sempre usar a segunda forma. Uma forma mais sintética de escrever uma matriz com elementos repetidos é possível indicando o número de vezes que o elemento deve ser repetido, entre colchetes. Por exemplo: pr oc i ml ; r e spos t a={ [ 3] ' Si m' , [ 3] ' Não' , [ 3] ' Tal v ez' } ; pr i nt r es pos t a;

RESPOSTA Si m Si m Si m Não Não Não Ta l v e z Ta l v e z Ta l v e z

I nt r odução oduçã o ao S A S / I M L par pa r a Á lgebr a de M at a t r iizes zes

3


Isso é m uito conveniente para constru ir matrizes com elemen tos repetidos, com o ocorrem com as matrizes de delineam entos experimen tais. Conforme mencionado, uma matriz pode ser redefinida no momento em que o usuário ju lga r co n ve n ie n te. A s egu n d a m at riz den om in ada RESPOSTA s obrep õe-se à p rim eir a. pr oc i ml ; r e spos t a={ [ 3] ' Si m' , [ 3] ' Não' , [ 3] ' Tal v ez' } ; r espos t a ={ [ 3] 1, [ 3] 2, [ 3] 0} ; pr i nt r es pos t a;

RESPOSTA 1 2 0

1 2 0

1 2 0

Matrizes também podem ser obtidas como resultado de fun ções ou e xpressões, como por exemp lo x= inv(a); em que os elem entos de X são obtidos pela inversão de A. Há ainda um a opção muito prática, que é a construção de vetores indexados, como apresentado abaixo. Note que é possível inclusive usar a função DO, sendo o último número da função o incremento. pr oc i ml ; r =1: 5; s=10: 5; t =' x yz 1' : ' u=do( - 2, 1, v =do( - 2, 1, pr i nt r s

xyz 4' ; 1) ; . 5) ; t u v;

R 1

2

3

4

5

S 10

9

8

7

6

5

0

0. 5

T xyz1 xyz2 xyz3 xyz4

4

U -2

-1

0

1

V -2

- 1. 5

-1

- 0. 5


1

D ef i nniinn d o m at r i ze s

Um com and o tam bém útil é o SHOW NAMES, que perm ite visualizar os atributos de tod as as m atrizes arm azenadas. No caso, tem -se a dim en são, o tipo - n um érica ou alfanum érica e o n úm er o de bit s par a arm azenam ent o. pr oc i ml ; x={ 1 1 0, 1 1 0, 1 0 1, 1 0 1} ; y ={ 1 , 2, 3, 4} ; show names ; pr i nt x y ;

X 4 r ows Y 4 r ows Nu mb er o f s y mb ol s = 2

X 1 1 1 1

1 1 0 0

3 col s num 8 1 col num 8 ( i n c l u de s t h os e wi t h ou t v al u es )

Y 1 2 3 4

0 0 1 1

Matrizes especiais, como a matriz identidade, possuem comandos específicos para sua obtenção, send o que o nú m ero entre parênteses indica sua dimen são: pr oc i ml ; i dent i d=i ( 5 ) ; pr i nt i dent i d;

I DENTI D 1 0 0 0 0

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

Já m atrizes com elem en tos iguais, que pod em ser iguais a um (m atriz de un s) ou a qualquer outro valor de int eresse, são obtidas com a fun ção J.


55


pr oc i ml ; uns 1_5=j ( 1, 5, 1) ; uns 5_5=j ( 5, 5, 1) ; quat r os =j ( 3, 4, 4) ; p r i n t u ns 1 _5 un s5 _5 qu at r o s ;

1

UNS1_5 1

1

1

1

1 1 1 1 1

UNS5_5 1 1 1 1 1

1 1 1 1 1

1 1 1 1 1

1 1 1 1 1

QUATROS 4 4 4 4 4 4

4 4 4

4 4 4

Para essa função, o primeiro valor dentro dos parênteses indica o número de linhas da matriz a ser criada, o segundo, o de colunas, e o terceiro, o valor que terão todos os elementos. Outros coman dos úteis são os que fornecem a dimen são das matrizes armazenadas, NROW (núm ero de lin has) e NCOL (n úm ero de colun as): pr oc i ml ; x={ 1 2 3, 4 5 6, [ 3] 9, [ 2] 1 0} ; m=nr ow( x) ; n=nc ol ( x ) ; pr i nt m n x;

M 4

N 3

X 1 4 9 1

2 5 9 1

3 6 9 0

Finalizando esta sessão, mencionaremos ainda o comando STORAGE, que permite armazenar (gravar) as matrizes numa biblioteca, permitindo seu uso posterior. Isso pode ser útil nos casos em que os programas são extensos e/ou há limitação de memória do comp utad or. O coman do FREE elim ina as m atrizes da m em ória RAM. pr oc i ml ;

6


O perações B ásicas

x={ 1 1 0, 1 1 0, 1 0 1, 1 0 1} ; y ={ 1 , 2, 3, 4} ; s t or e x y; f r ee x y; show st or age;

Cont ent s of st or age l i br ar y = WORK. I MLSTOR Mat r i c es : X Y Modul es:

3

Operações Básicas A partir desse ponto, em função de estarm os sem pre apresen tando com andos d o SAS/IML, não apresentaremos mais os comandos PROC IML e QUIT, necessários para iniciar e en cerrar a sessão.

3.1. Adição A={ 2 2, 3 4} ; B={ 4 5, 1 0} ; SOMA=A+B; pr i nt A B soma;

A 2 3

2 4

B 4 1

5 0

SOMA 6 4

7 4


7


3.2. Subtração DI F=A- B; pr i nt di f ;

DI F -2 2

-3 4

3 . 3 . M u l t i p l i c a çã o p o r Esc al a r C=4; AxC=A* C; pr i n t Ax C;

AXC 8 12

8 16

3 . 4 . M u l t i p l i c aç ão d e Ma t r i z e s AxB=A* B; pr i n t Ax B;

AXB 10 16

10 15

3 . 5 . So m a Di r e t a Para tanto, deve-se utililizar a função BLOCK, que permite realizar a soma direta de um grande nú m ero de m atrizes. a={ 2 2, 4 4} ; b={ 6 6, 8 8} ; c =bl oc k( a, b) ;

8



pr i nt a b c ;

A 2 4

B 6 8

2 4

C 2 4 0 0

6 8

2 4 0 0

0 0 6 8

C2 1 0 2 0 3 0

0 1 0 2 0 3

0 0 6 8

3 . 6 . Pr o d u t o Di r e t o ( Kr o n e c k er ) a={ 1 0, 0 1} ; b={ 1, 2, 3} ; c1=a@b; c2=b@a; c3=a@a; pr i n t a b c1 c2 c 3;

A 1 0

C3 1 0 0 0

B 1 2 3

0 1

0 1 0 0

C1 1 2 3 0 0 0 0 0 1 0

0 0 0 1 2 3 0 0 0 1

3 . 7 . Pr o d u t o d e Ha d am a r d a={ 1 2, 3 4} ; b={ 4 8, 0 5} ; c=a#b; pr i nt c;


9


C 4 0

16 20

3 . 8 . Po t ê n ci a d e M a t r i z Qu a d r a d a É im portante ressaltar que m atrizes elevadas a nú m eros m uito grandes p odem apresentar como resultado valores com grandes imprecisões numéricas. A potência é obtida como segue. a={ 1 2, 3 4} ; a2=a* * 2; a3=a* * 3; i nv er s a=a* * ( - 1) ; pr i nt a2 a3 i nv er s a;

A2 7 15

10 22

A3 37 81

54 118

I NVERSA -2 1. 5

1 - 0. 5

3.9. Traço x={ 1 2 4 5 6 2 t r ac o=t pr i nt t

3, 6, 4} ; r ac e( x ) ; r aco;

TRACO 10

3.10. Transposição x ={ 1 . 9 2 3, 40 5. 2 6, 60 20 4. 3} ; t r ans p1=t ( x ) ; t r an sp 2=x ` ;

10



pr i nt t r ans p1 t r ans p2;

TRANSP1 1. 9 2 3

40 5. 2 6

60 20 4. 3

TRANSP2 1. 9 2 3

40 5. 2 6

60 20 4. 3

Note que há duas formas diferentes de tran spor um a matriz, um a com o argum ento pósfixado e um a com o argum ento pré-fixado. Evidentem ente, ambas produzem o m esm o resultado.

3.11. Comparações entre Matrizes O SAS/IML possui operadores que permitem realizar a comparação entre os elementos de duas m atrizes, fornecendo como resultado um a nova m atriz de un s e zeros, conform e a comp aração seja verdadeira ou falsa, respectivam en te. A sintaxe dos op eradores é: m atriz1 < m at r iz2, mat r iz1 < = m at ri z2, m at ri z1 > m at ri z2, m at r iz1 > = m at ri z2, m at riz1 = m atriz2, m atriz1 ^ = m atriz2. Os dois últim os operadores permitem iden tificar matriz iguais, mas devem os sem pre obser var que isso pode ocorrer devido a erros de aproxim ação. Algun s exem plos: a={ 1 7 3, 4 2 4} ; b={ 0 8 2, 4 1 3} ; menor = ab; mai or i g= a>=b; i g ua l = a =b ; di f er ent = a^ =b; pr i nt menor menor i g mai or mai or i g i gual

MENOR 0 0 MAI OR 1 0 I GUAL 0 1

1 0

0 1

0 0

di f er ent ;

0 0

MENORI G 0 1

1 0

0 0

1 1

MAI ORI G 1 1

0 1

1 1

0 0

DI FERENT 1 0

1 1

1 1

I nt r odução oduçã o ao S A S / I M L par pa r a Á lgebr a de M aatt r iizes zes

11


3 . 1 2 . Co n ca t e n a çã o d e M a t r i z e s e Ve t o r e s A concatenação pode ser feita horizontalmente ou verticalmente, fornecendo uma excelen te ferram en ta para análise de dados em situações reais. Vejam os algun s exem plos: A={ 4 3, -2 5}; b={ 10, 8} ; C={ 7 - 10, 4 13 } ; A_b=A| | b; A_C=A/ / C; pr i nt A_b A_C;

A_B 4 -2

3 5

10 8

A_C 4 -2 7 4

3 5 -10 13

A m atriz A e o vetor b foram concatenados horizontalmen te usand o o operador | | . Já as matrizes A e C foram concatenadas verticalmente usando //, obviamente observando se as dimen sões são comp atíveis.

3 . 1 3 . Esc al o n a m e n t o d e Ma t r i z e s O SAS/IML possui duas funções para obtenção de formas escalonadas, HERMITE e ECHELON, resultando m atrizes na form a de Herm ite e Echelon, respectivam en te. x={ 3 6 9, 1 2 5, 2 4 10} ; h=he r mi t e ( x ) ; e=ec he l on ( x ) ; pr i nt h e;

H 1 0 0

2 0 0

0 0 1

E 1 0 0

2 0 0

0 1 0

Note que, para essas matrizes, as formas escalonadas resultantes são diferentes, sendo a segund a delas a form a escalonada canôn ica.

12



3 . 1 4 . P o st o d e M a t r i z e s As formas escalonadas perm item d eterminar o posto de m atrizes: A={ 4 2 2, 2 2 0, 2 0 2} ; B={ 3 1 2, 1 4 1, 2 1 5} ; C={ 1 2 3 4, 2 4 6 8} ; Aesc=echel on( A) ; Besc=her mi t e( B) ; Cesc=echel on( C) ; pr i nt Aesc Besc Cesc;

AESC 1 0 0

0 1 0

1 -1 0

CESC 1 0

2 0

3 0

BESC 1 0 0

0 1 0

0 0 1

4 0

Verifica-se qu e o p osto d e A, B e C é 2, 3 e 1, resp ectivam en te. Embora o SAS/IML não possua um a função que forneça diretamen te o posto de um a m atriz (a fun ção RANK fornece o ran k dos elem en tos, e não d a m atriz, com o desejado), isso pod e ser conseguido com o au xílio de u m programa auxiliar, apresen tado a se guir. Ele bas eia-se n os teorem as 55 (r[A]= r[A+ A]) e 50 (A+ A é simé trica e idem poten te), lembran do ainda qu e para m atrizes idem poten tes o ran k é igual ao traço. A fun ção ROUND arredond a o valor encon trado, que é sujeito a pequen os erros de aproxim ação. A={ 4 2 2, 2 2 0, 2 0 2} ; B={ 3 1 2, 1 4 1, 2 1 5} ; C={ 1 2 3 4, 2 4 6 8} ; RANK_ A=r o un d( t r a ce ( g i n v( A) * A) ) ; RANK_ B=r o un d( t r a ce ( g i n v( B) * B) ) ; RANK_C=r ound( t r ace( gi nv( C) * C) ) ;


13

I n ver sa de d e M a t r i z N ã o S i ng ul a r pr i nt

RANK_A RANK_B RANK_C;

RANK_A 2

4

RANK_B 3

RANK_C 1

Inversa de M at riz Não Singular Apresentaremos a seguir três formas diferentes de calcular a inversa de matriz não singular, send o que a prim eira delas (INV1) foi incluída ap en as por razões didáticas, usan do o algoritm o de Gauss. B={ 3 1 2, 1 4 1, 2 1 5} ; I NV1=her mi t e( B| | I ( 3) ) ; I NV2 =i n v( B) ; I NV3 =B* * ( - 1 ) ; pr i nt I NV1 I NV2 I NV3;

I NV1 1 0 0

0 1 0

I NV2 0. 4 75 - 0. 0 75 - 0. 1 75

- 0. 0 75 0. 2 75 - 0. 0 25

0 0 1

0 . 4 75 - 0 . 0 75 - 0 . 1 75

- 0. 1 75 - 0. 0 25 0. 2 75

- 0 . 0 75 0 . 2 75 - 0 . 0 25

I NV3 0. 4 75 - 0. 0 75 - 0. 1 75

- 0 . 1 75 - 0 . 0 25 0 . 2 75

- 0. 0 75 0. 2 75 - 0. 0 25

- 0. 1 75 - 0. 0 25 0. 2 75

Aproveitaremos a matriz INV1, que contém a inversa de B à direita da identidade, para m ostrarm os com o selecionar subm atrizes de um a dada m atriz. É possível selecionar linh as, colunas ou até mesmo elementos de interesse. Para tanto, utilizam-se subscritos,que são operadore s localizados den tros de colchetes, após o operand o. O primeiro nú m ero refere-se a(s) linha(s) a ser(em) retirada(s) e o segundo a(s) coluna(s), sendo que eles devem ser separados por uma vírgula. Para selecionar um único elemento, como por exemplo da 3a linh a e 4a colun a (lem bre-se que INV1 tem dim en são 3 x 6): B={ 3 1 2, 1 4 1,

14


I n ver ssaa ddee M a t r i z N ã o S i ng ul a r

2 1 5} ; I NV1=her mi t e( B| | I ( 3) ) ; el em3_4=I NV1[ 3, 4] ; pr i n t el e m3_4 ;

ELEM3_4 - 0. 175

Para selecionar a linh a 3 ou a colun a 4: B={ 3 1 2, 1 4 1, 2 1 5} ; I NV1=her mi t e( B| | I ( 3) ) ; l i nha3=I NV1[ 3, ] ; c ol u na 4=I NV1[ , 4] ; pr i nt l i nha3 c ol una4;

LI NHA3 0

0

1

- 0 . 1 75

- 0 . 0 25

0 . 2 75

COLUNA4 0 . 4 75 - 0. 075 - 0. 175

Para selecionar u m a subm atriz, com o por e xemplo a inversa sem a parte referente à m atriz identidade: B={ 3 1 2, 1 4 1, 2 1 5} ; I NV1=her mi t e( B| | I ( 3) ) ; I NVERSA=I NV1[ , { 4 5 6} ] ; pr i nt I NVERSA;

I NVERSA 0. 4 75 - 0. 0 75 - 0. 1 75

- 0. 0 75 0. 2 75 - 0. 0 25

- 0. 1 75 - 0. 0 25 0. 2 75

Note qu e é ne cessário u tilizar chaves p ara correta especificação da sub-matriz.


15

D et er m i na nt es

5

Determ inantes O SAS/IML calcula o determinante de matrizes quadradas de forma simples, através da função DET, considerando o fato de que o determ inante pode ser obtido através do produto dos au tovalores da m atriz em questão, que p or sua vez são obtidos pela decomp osição LU (matrizes triangulares). O resultado fornecido é uma matriz 1x1 (escalar) contendo o det er mi nant e. pr oc i ml ; a={ 1 1 1, 1 2 4, 1 3 9} ; d et e r m=d et ( a ) ; pr i nt det er m;

DETERM 2

Não há uma função específica para se obter matrizes adjuntas e dos cofatores de uma m atriz de interesse.

6

Sist em as Consist ent es de Equações Lineares O SAS possui a fu n ção SOLVE para a solução de sistem as de equações lineares (consisten tes e determ inados), que é m uito eficiente em termos com putacionais, devendo ser u sado em situações práticas. Contudo, sempre com finalidades didáticas, serão apresentadas outras form as de resolver sistem as. Para usar a fu nção SOLVE para resolver o sistem a Ax= b, pode-se proceder como apresen tado a seguir. A={ 4 3, - 2 5} ; b={ 10, 8} ; x =s ol v e( A, b) ; pr i nt x;

X

16


S ist emas C onsist ent es de E quações L inear es

1 2

As formas escalonadas apresen tadas também podem ser usadas para solução de sistemas de equações lineares. Considerando o mesmo sistema com duas equações e duas incógnitas, escrito na form a m atricial Ax= b, um programa p ara obter sua solução pode ser: A={ 4 3, - 2 5} ; b={ 10, 8} ; A_b=A| | b; SOL1=echel on( A_b) ; SOL2=her mi t e( A_b) ; pr i nt SOL1 SOL2;

O escalonam en to forneceu a solução do sistem a para x, com resu ltados 1 e 2. SOL1 1 0

0 1

1 2

SOL2 1 0

0 1

1 2

Caso desejem os resolver o sistem a usan do o m étodo de Kramer, podem os usar o seguinte programa: r eset pr i nt ; A={ 1 1 2, 2 0 3, 1 1 0} ; B={ 9, 11, 3} ; A1=B | | A[ , { 2 3} ] ; A2=A[ , 1] | | B | | A[ , 3] ; A3=A[ , { 1 2} ] | | B; DET_A=det ( A) ; DET_A1=det ( A1) ; DET_A2=det ( A2) ; DET_A3=det ( A3) ; X1=DET_A1/ DET_A; X2=DET_A2/ DET_A; X3=DET_A3/ DET_A; SOLUCAO=X1/ / X2/ / X3;

Vale lembrar que não houve necessariamente a preocupação em elaborar programas mais sint éticos e eficiente s. Nesse con text o, foram defin idas as m atrizes A, B, A1, A2 e A3, com


17


seus re spectivos determ inan tes. Em seguida, foi obtida a solução do sistem a, concaten ada no vetor SOLUCAO. O comando RESET PRINT solicita a visualização de todas as matrizes criadas, facilitando a visualização dos resultados, em bora com excesso de in form ações. Para des ativá-lo, deve-se u sar RESET NOPRINT. A saída do pro gram a acim a é

18

A

3 r ows 1 2 1

3 col s 1 0 1

( numer i c) 2 3 0

B

3 r ows 9 11 3

1 col

( numer i c)

A1

3 r ows 9 11 3

3 col s 1 0 1

( numer i c) 2 3 0

A2

3 r ows 1 2 1

3 col s 9 11 3

( numer i c) 2 3 0

A3

3 r ows 1 2 1

3 col s 1 0 1

( numer i c) 9 11 3

DET_A

1 r ow 4

1 col

( numer i c)

DET_A1

1 r ow 4

1 col

( numer i c)

DET_A2

1 r ow 8

1 col

( numer i c)

DET_A3

1 r ow 12

1 col

( numer i c)

X1

1 r ow 1

1 col

( numer i c)

X2

1 r ow 2

1 col

( numer i c)

X3

1 r ow 3

1 col

( numer i c)



SOLUCAO

3 r ows 1 2 3

1 col

( numer i c)

Um programa m ais sintético pode ser: r es et nopr i nt ; A={ 1 1 2, 2 0 3, 1 1 0} ; B={ 9, 11, 3} ; A1=B | | A[ , { 2 3} ] ; A2=A[ , 1] | | B | | A[ , 3] ; A3=A[ , { 1 2} ] | | B; SOLUCAO=DET( A1) / DET( A) / / DET( A2) / DET( A) / / DET( A3) / DET( A) ; pr i nt s ol uc ao;

SOLUCAO 1 2 3

No caso de obter soluções para sistem as de equações lineares h om ogêneos do tipo Ax= (, sendo A de dimensão m x n, mn, é possível utilizar a função HOMOGEN, que retorna soluções com o colun as ortonorm ais de u m a mat riz X com as prop riedades Ax= ( (P1) e X'X= I (P2). Para tanto, o posto d e A deve ser r< n , ou seja, A deve ser de p osto incom pleto. A dim en são dessa m atriz X é n x (n -r), e portan to deve ser dad a especial aten ção a matrizes m al-condicion adas, já que o correto posto de A pode n ão ser obtido em fun ção de erros de aproximação. A={ 4 2 2, 2 2 0, 2 0 2} ; X=HOMOGEN( a) ; P1=A* X; P2=T( X) * X; pr i nt X P1 P2;

X - 0. 57735 0. 5773503 0. 5773503

P1 2. 22E- 16 2. 22E- 16 0

P2 1

Nesse caso, x é a solução não trivial, e P1 e P2 a verificação das propriedades da solução. Note que pequen os erros de aproxim ação podem ocorrer.


19

R aízes C ar act eríst icas (aut ovalores)

Um ou tro fato que pode ser m encionado é que a dim ensão da m atriz com as soluções do sistem a hom ogêneo podem ser usadas para calcular o posto da m atriz A, já que n ão há u m a fun ção específica para tanto : POSTO_A=NROW( X) - NCOL( X) ; pr i nt POSTO_A;

POSTO_A 2

7

Raízes Características (autovalores) Para obter os autovalores ("eingenvalues") de matrizes quadradas simétricas, usa-se a função EIGVAL, que fornece um vetor com os autovalores ordenados em ordem decrescente. A={ 3 1 1, 1 3 1, 1 1 3} ; AUTOVAL=EI GVAL( A) ; pr i nt AUTOVAL;

AUTOVAL 5 2 2

Caso necessitem os obter os au tovalores a partir da e quação característica, podem os usar a função POLYROOT, que encontra as raízes de um polinômio real. A sintaxe é POLYROOT(vetor), em que o vetor de dimensão 1 x n contém os coeficientes de um polinômio de grau n-1, sendo que eles devem estar em ordem decrescente, com relação as potên cias do polinôm io. O resultado é u m a m atriz de dim ens ão (n -1) x 2, com as raízes do polinômio; a primeira coluna contém a parte real e a segunda, a parte imaginária (que deverá ser zero, no cont exto des se curso). No exem plo em questão, a equ ação característica é -x3+ 9x2-24x+ 20= 0, e assim p={ - 1 9 - 24 20 } ; r ai zes =POLYROOT( p) ;

20

I nt r oduçã oduçãoo ao S A S / I M L par zes pa r a Á lgebr a de M a t r iizes

R aízes C ar act eríst es) er íst icas (aut ovalor ovalores)

pr i nt r ai zes;

RAI ZES 2. 0000003 1. 9999997 5

0 0 0

É interessante notar que as soluções estão sujeitas a erros de aproximação, em função do algoritm o em pregado. Um outro asp ecto que deve ser lem brado diz respeito ao fato do SAS/ IML usar um algoritmo que não garante que todas as raízes do polinômio sejam en contradas; nesses casos, há u m a m en sagem para o us uário na jane la LOG. Nas situações em que desejamos obter matrizes diagonais com as raízes características de uma dada m atriz, com bina m os a fun ção DIAG, que cria m atrizes diagona is, com a fun ção EIGVAL. M={ 2 - 1 - 1, - 1 2 - 1, - 1 - 1 2} ; AUTOV=EI GVAL( M) ; L=DI AG( AUTOV) ; PRI NT AUTOV L;

AUTOV

L

3 3 - 9 . 1 1E- 1 7

3 0 0

0 0 3 0 0 - 9 . 1 1E- 1 7

É fácil verificar que, exceto de vido a erros de ap roxim ação, M é po sitiva sem i-defin ida. Caso apliquem os o coman do DIAG em um a m atriz (e n ão nu m vetor, com o n o exemp lo acima), terem os como resultado um a m atriz form ada pelos elem entos da diagonal dessa m atriz, sendo os dem ais substituídos por zero. Para facilitar a visualização, poderíamos ainda usar a opção FUZZ, associada ao comando RESET, de modo que os números muito pequenos (por volta de 1E-12) sejam impressos como zero. r es et f uz z; PRI NT AUTOV L;

AUTOV 3 3

L 3 0

0 3

0 0


21 21

V et or es C a r a c t er í st i c os ( a ut ov et or es) 0

0

0

0

Para voltarm os ao default, dent ro da m esm a seção do IML, devem os usar o com and o RESET NOFUZZ.

8

Vetores Característicos (autovetores) A obtenção dos autovetores ("eingenvectors") de matrizes quadradas simétricas pode ser feita pelo uso da função EIGVEC, que cria uma matriz cujas colunas são autovetores ortonormais. A={ 3 1 1, 1 3 1, 1 1 3} ; AUTOVEC=EI GVEC( A) ; pr i nt AUTOVEC;

AUTOVEC 0. 5773503 0. 5773503 0. 5773503

- 0. 408248 - 0. 408248 0. 8164966

0. 7071068 - 0. 707107 0

É fácil verificar que esses vetores são ortogonais en tre si, pois o produ to en tre au tovalores diferentes é 0, apresentando como resultado de seu produto consigo mesmo o valor 1 (n orm a 1). Em outras palavras, o produ to ent re a transposta da m atriz com os aut ovalores e ela própria é a iden tidade: RESET FUZZ; AUTOVEC=EI GVEC( A) ; ORTONORM=AUTOVEC` * AUTOVEC; pr i nt ORTONORM;

ORTONORM 1 0 0 1 0 0

0 0 1

Em situações práticas, pode m os obter os aut ovalores e autovetores diretam en te, com o uso da sub-rotina EIGEN, acionada pelo comando CALL, que invoca determinadas estruturas

22


F or m maass E q ui v al en t e s

(operações) pré-armazenadas. O primeiro termo dentro dos parênteses especifica o nome da matriz que conterá os autovalores; o segundo, da matriz que conterá os autovetores; oterceiro ind ica sobre qual m atriz o coman do será ap licado. CALL EI GEN( AUTOVAL, AUTOVET, A) ; pr i nt AUTOVAL AUTOVET;

AUTOVAL 5 2 2

9

AUTOVET 0. 5773503 - 0. 408248 0. 7071068 0. 5773503 - 0. 408248 - 0. 707107 0. 5773503 0. 8164966 0

Form as Equivalentes As m atrizes A e B são equivalentes se existem M e N n ão singulares tais que M AN= B. B é dita forma normal de A. O SAS/IML fornece a forma normal de Hermite da matriz A, usando função já vista anteriormente. Para matrizes quadradas, essa forma é triangular sup erior e idem poten te. Se A for não singular, o resultado será um a m atriz identidad e. A={ 4 2 2, 2 2 0, 2 0 2} ; B=h er mi t e ( A) ; pr i nt B;

B 1 0 0

0 1 0

1 -1 0

É possível tam bém obter a forma n orm al de Echelon que, nesse caso, coincide com a forma norm al de Hermite: B=echel on( A) ; pr i nt B;

B


23 23

F or m maass E q uuiivv aallen en t e s 1 0 0

0 1 0

1 -1 0

Se a m atriz em qu estão for quadrada, as form as norm ais de Herm ite e Echelon p odem ser obtidas um a a partir da outra, apen as rearranjando as linh as de zeros. A form a de Eche lon sem pre fornece um a forma escalonada canônica. Para a obtenção de u m a m atriz D diagon al congruen te (e, portan to, equivalente) à m atriz A real e sim étrica con siderad a, tal que TAT '= D, utiliza-se a su b-rotin a EIGEN, pelo coma n do CALL EIGEN(D, T, A), já descrito anteriormente. Para tanto, utilizam-se as seguintes propriedad es das m atrizes obtidas com esse com and o: T' = T-1 (ortogonal) e A= T*diag(D)*T', ou se ja, diag(D) é um a m atriz diagon al con gruen te à m atriz A. A={ 4 2 2, 2 2 0, 2 0 2} ; CALL EI GEN ( D, T, A) ; CONG=DI AG( D) ; A2=T*DI AG( D) * T` ; r eset FUZZ; pr i n t A A2;

A 4 2 2

2 2 0

2 0 2

A2 4 2 2

2 2 0

2 0 2

Dessas prop riedades, conseguimos tam bém verificar que d iag(D)= T'*A*T, em que d iag(D) é u m a m atriz diagonal que exibe as raízes características de A, e T é um a m atriz ortogonal cujas linhas são vetores característicos normalizados associados as raízes características de A. L1=DI AG( D) ; L2=T` * A* T; pr i nt L1 L2;

L1 6 0 0

0 2 0

0 0 0

L2 6 0 0

0 2 0

0 0 0

Nesse caso, A e diag(D) são tam bém sim ilares.

24


F a t or a çã o

10

Fatoração O SAS/IML n ão realiza diretam en te fatorações do tipo A= RR', fatorações de posto comp leto do tipo A= BC e operações de d iagonalização sim ultân ea. No prime iro caso, é possível a obten ção de fatorações do tipo raiz quadrad a de out ra forma, apresen tada abaixo. Uma vez que o principal interesse na fatoração de posto completo é para a obtenção da inversa gen eralizada de Moore-Pen rose, a ausência da fatoração não rep resen ta em pecilho, um a vez que essa in versa generalizada é calculada pelo SAS (itens s eguintes). No m anu al, tem os diversas outras decom posições que pode m interessar ao leitor. Fatorações do tipo A= U'U são obtidas com a fu n ção ROOT, que re torna u m a m atriz triangular superior U que satisfaz a igualdade, sendo que a matriz A deve ser simétrica e positiva definida. Essa decomposição é conhecida como decomposição de Cholesky, e a m atriz U é dita raiz quadrada de A. U=ROOT( A) ; A1=U` * U; p r i n t U A1 ;

U 2 0 0

1 1 0

1 -1 0

A1 4 2 2

2 2 0

2 0 2

Vale ressaltar que algun s autores d en om inam a fatoração A= RR' como raiz quadrada, e não A= U'U, como o SAS con sidera. É interessan te n otar que essa fatoração também é do tipo triangular (A= LU), um a vez que U' é triangular inferior e U é trian gular superior. Para realizar de comp osição p or valores singulares do tipo X= USV´ , o SAS/IML possu i a su brotina SVD, acessada pelo comando CALL. Sua sintaxe é CALL SVD(U, S, V, X), em que X represen ta a m atriz a ser decom posta, U, S e V, o nom e das m atrizes a serem criadas com o resultado da d ecomp osição, de forma qu e X= U*diag(S)*V´ e S cont en ha os valores singulares de A. Note que o algoritmo apresentado no texto das aulas teóricas produz m atrizes diferentes, m ostrando que a decomposição não é ún ica. r eset pr i nt ; X={ 1 1 0, 1 1 0, 1 0 1, 1 0 1} ; c al l s vd( U, S, V, X) ; X1=U* di ag( S) * V` ;


25 25

Inver sas G eneralizadas

X 1 1 1 1

4 r ows 1 1 0 0

3 col s 0 0 1 1

( numer i c)

U 0. 5 0. 5 0. 5 0. 5

4 r ows 3 col s 0 . 5 0 . 7 07 10 68 0. 5 - 0. 70 7107 - 0. 5 0 - 0. 5 0

( numer i c)

S 3 r ows 2. 4494897 1. 4142136 0

1 col

( numer i c)

V 3 r ows 3 col s ( numer i c) 0. 8164966 0 - 0. 57735 0. 4082483 0. 7071068 0. 5773503 0. 4082483 - 0. 707107 0. 5773503 X1 1 1 1 1

11

4 r ows 1 1 0 0

3 col s 0 0 1 1

( numer i c)

Inversas Generalizadas É possível calcular a inversa generalizada de Moore-Penrose usando o SAS/IML, com o comando GINV(X), em que X é a matriz a ser invertida. A inversa obtida satisfaz as quatro condições da inversa, conform e apresen tado abaixo. X={ 1 1 0, 1 1 0, 1 0 1, 1 0 1} ; MP1=gi nv( X) ; P1=MP1* X* MP1; P2=X* MP1* X; P3a=X* MP1; P3b=t ( X* MP1) ; P4a=MP1* X; P4b=t ( MP1* X) ; pr i nt MP1;

26


Inver sas G eneralizadas

pr i nt P1 P2 P3a P3b P4a P4b;

MP1 0. 1666667 0. 1666667 0. 1666667 0. 1666667 0. 3333333 0. 3333333 - 0. 166667 - 0. 166667 - 0. 166667 - 0. 166667 0. 3333333 0. 3333333 P1 0. 1666667 0. 1666667 0. 1666667 0. 1666667 0. 3333333 0. 3333333 - 0. 166667 - 0. 166667 - 0. 166667 - 0. 166667 0. 3333333 0. 3333333

P2 1 1 1

1 1 0 1

P3A 0. 5 0. 5 0 0

0. 5 0. 5 0 0

0 0 0. 5 0. 5

P3B 0. 5 0. 5 0 0

0 0 0. 5 0. 5

P4A 0. 6666667 0. 3333333 0. 3333333 0. 3333333 0. 6666667 - 0. 333333 0. 3333333 - 0. 333333 0. 6666667

0. 5 0. 5 0 0

0 0 1 0

0 0 0. 5 0. 5

1 0 0 0. 5 0. 5

P4B 0. 6666667 0. 3333333 0. 3333333 0. 3333333 0. 6666667 - 0. 333333 0. 3333333 - 0. 333333 0. 6666667

As inversas generalizadas reflexivas g2 para matrizes de delineamento podem ser obtidas com a fun ção SWEEP. Tom and o o m odelo line ar de Gauss-Markov y = XB + e, pod em os obter o sistem a X'X X'y y 'x y 'y

e, usand o o operad or SWEEP sobre esse sistem a, calcular um a m atriz ( X 'X )

g2

– y ' X( X ' X )

( X'X ) g2

g2

X'y g2

y ' ( I – X ( X 'X ) X ' ) y

A partição superior esqu erda d essa m atriz contêm a inversa g2 de X'X, a partição sup erior direita as soluções para B e apartição inferior direita a som a de qu adrados d o resídu o. Esse é o procedimento usado pelo PROC GLM do SAS para realizar análises de variância. Um exemplo é apresentado abaixo, sendo que o segundo argumento da função SWEEP é um vetor com índices variand o de 1 até o p osto de X. r eset pr i nt ; X={ 1 1 0, 1 1 0, 1 0 1, 1 0 1} ; y ={ 1 ,


27 27

Inver sas G eneralizadas 2, 4, 3} ; A=( t ( X) * X| | t ( X) * y ) / / ( t ( y ) * X| | t ( y ) * y ) ; t =1: nc ol ( X) ; G=sweep( A, t ) ; r es et nopr i nt ;

X 1 1 1 1 Y

4 r ows 1 1 0 0

3 col s

( numer i c)

0 0 1 1

4 r ows

1 col

( numer i c)

4 r ows

4 col s

( numer i c)

1 2 4 3 A 4 2 2 10 T 1 G 0. 5 - 0. 5 0 - 3. 5

2 2 0 3 1 r ow 2 4 r ows - 0. 5 1 0 2

2 0 2 7

10 3 7 30

3 col s

( numer i c)

3 4 col s 0 0 0 0

( numer i c) 3. 5 -2 0 1

Leitores já familiarizados com o PROC GLM podem verificar as afirmativas acima com o seguinte programa: d at a ; i nput t r at Y@@; dat al i nes ; 1 1 1 2 2 4 2 3 ;

28


Apostila Proc Iml

Recommend Documents