Apostila Proc Iml

Introdução ao SAS/IML para Álgebra de

Matrizes

Escola Superior de Agricultura “Luiz de Queiroz”Departamento de Ciências Exatas - ESALQ/USP

Fevereiro de 2002( ve r s ã o 3 , 14 )

Antonio Augusto Franco GarciaAntonio Augusto Franco GarciaAntonio Augusto Franco GarciaAntonio Augusto Franco GarciaAntonio Francisco IemmaAntonio Francisco IemmaAntonio Francisco IemmaAntonio Francisco Iemma

Introdução ao SAS/IML para Álgebra de Matrizes

Esse texto destina-se aos alunos do curso de verão de Álgebra de Matrizes doDepartamento de Ciências Exatas da ESALQ/USP, visando auxiliar naresolução dos exercícios propostos e fornecendo ao mesmo tempo subsídiospara a utilização de um dos programas computacionais mais empregados naEstatística: o SAS (Statistical Analysis System).

Optou-se por apresentar as possíveis aplicações do módulo IML mantendo amesma seqüência do texto usado em aula, de modo a complementá-lo.Leitores interessados em maiores detalhes poderão consultar o livro SAS/IMLSoftware: Usage and Reference, Version 6, First Edition, Cary, NC: SASInstitute Inc., 1989. 501pp. (6th printing, December 1995), além de outraspublicações do referido instituto.

• Por ter finalidades didáticas, os programas aqui citados foram elaboradosde forma a facilitar o entendimento do assunto, não visandoexclusivamente aplicações práticas ou análise de dados. Supõe-se contudoque esta forma de abordagem propiciará aos leitores com interesseaplicado os subsídios necessários para o desenvolvimento de seuspróprios programas.

• Críticas e sugestões, bem como correção de eventuais erros no material,serão bem recebidas.

Cons ide raçõe s sobre o SAS/IMLCons ide raçõe s sobre o SAS/IMLCons ide raçõe s sobre o SAS/IMLCons ide raçõe s sobre o SAS/IML

1 Considerações sobre o SAS/IMLO software SAS/IML (Interactive Matrix Language) é uma linguagem de programação muitoflexível, tendo como principais vantagens a interatividade com o usuário, a possibilidadede trabalhar com matrizes de forma simples, além da possibilidade do uso de funções pré-programadas e outras procedures do próprio SAS. Os resultados dos comandos podem servistos imediatamente ou armazenados para posterior execução, a critério do usuário. Nessetexto, por razões óbvias, será dada ênfase as operações matriciais. No entanto, o SAS/IML éuma linguagem de programação, permitindo que o usuário escreva seus própriosprogramas.

Operações como inversão de matrizes, obtenção de autovalores e autovetores, etc, podemser realizadas usando comandos pré-programados, facilitando o trabalho e evitando o usode algoritmos sofisticados que utilizam várias linhas de programação. O fato dos algoritmosempregados serem amplamente testados e discutidos por usuários do mundo todo permiteque os resultados obtidos sejam muito confiáveis. Além disso, o SAS/IML faz uso muitoeficiente dos recursos do computador, facilitando o trabalho com matrizes de grandesdimensõ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 asmudanças de linhas. Um primeiro programa:

proc iml;a={3 -1 2, 2 -2 3, 4 1 -4};x=inv(a);print x;quit;

X 0.3333333 -0.133333 0.0666667 1.3333333 -1.333333 -0.333333 0.6666667 -0.466667 -0.266667

Esse programa permite a obtenção da matriz X, que é a inversa de A. O comando PRINTsolicita a visualização da matriz X. Vale ressaltar que, ao usar o comando QUIT, a sessão doIML é encerrada e as matrizes e variáveis criadas são apagadas da memória RAM.

I n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze s 1111

Defin indo matr i z e sDef i n indo matr i z e sDef i n indo matr i z e sDef i n indo matr i z e s

2 Definindo matrizesO objeto fundamental sobre o qual operam os comandos do IML é uma matriz com linhas ecolunas, podendo ser numérica ou alfanumérica. O número máximo de caracteres que umelemento de uma matriz alfanumérica pode apresentar é 32676, sendo que todos oselementos serão armazenados com o mesmo número de caracteres que o maior deles. Osnomes 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 partirdo 2º caractere (versões 6.12 e anteriores).

Em versões mais recentes, como a 8.0 e 8.2, os nomes das variáveis podem exceder oitocaracteres. É possível criar matrizes com valores perdidos, sendo que os mesmos devem serindicados por um ponto ("."). Contudo, nem todos os comandos funcionam com matrizescom pontos.

Uma vez definido um nome para a matriz, ele pode ser usado diretamente ou substituídoem ações futuras, a critério do usuário, sendo possível inclusive substituir apenasdeterminados elementos, conforme será visto a seguir. Ao se referir a uma dada matriz,pode ser simplesmente utilizado seu nome, ou ainda especificando seus elementos.

Matrizes com múltiplos elementos são especificadas com o uso de chaves ({ }), sendo queas linhas são separadas por vírgulas. Evidentemente, todas as linhas devem conter omesmo número de elementos.

Os valores a serem adicionados podem ser:

• Números reais inteiros ou com decimais (notação científica, como 1E-8, é permitida).

• Caracteres alfanuméricos ("strings"), que podem estar entre aspas simples (´ ´) ouduplas (" "), ou até mesmo sem elas, sendo que nesse caso todos os valores serãoarmazenados com letras maiúsculas.

• Números entre colchetes ( [ ] ) indicam o número de vezes que um dado valor de repete,facilitando a digitação em muitas situações.

Matrizes também podem ser criadas a partir do uso de funções pré-programadas, docomando CALL ou de um programa específico. A seguir, são apresentados exemplos dematrizes com um único elemento (escalares):

proc iml;a=32;b=.;c='oi';d='olá pessoal';e="ciao";print a b c d e;

A B C D E 32 . oi olá pessoal ciao

2222 Introdução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze s

Def in i ndo matr i ze sDef i n i ndo matr i ze sDef i n i ndo matr i ze sDef i n i ndo matr i ze s

Já as matrizes com múltiplos elementos são escritas da seguinte maneira:

proc iml;x={10 20 30 40};y={10, 20, 30, 40};z={1 2, 3 4, 5 6};print x y z;

X Y Z 10 20 30 40 10 1 2 20 3 4 30 5 6 40

Ou ainda, para caracteres alfanuméricos: este

proc iml;aa={ abcd Efgh};bb={'abcd' 'Efgh'};print aa bb;

AA BB ABCD EFGH abcd Efgh

Pelo fato da matriz AA não apresentar seus elementos entre aspas, eles automaticamentesão substituídos por letras maiúsculas. Já para BB, foi respeitada a forma original. Caso sejanecessário definir elementos com espaços entre eles, deve-se sempre usar a segundaforma.

Uma forma mais sintética de escrever uma matriz com elementos repetidos é possívelindicando o número de vezes que o elemento deve ser repetido, entre colchetes. Porexemplo:

proc iml;resposta={[3] 'Sim', [3] 'Não', [3] 'Talvez'};print resposta;

RESPOSTA Sim Sim Sim Não Não Não Talvez Talvez Talvez



Isso é muito conveniente para construir matrizes com elementos repetidos, como ocorremcom as matrizes de delineamentos experimentais.

Conforme mencionado, uma matriz pode ser redefinida no momento em que o usuáriojulgar conveniente. A segunda matriz denominada RESPOSTA sobrepõe-se à primeira.

proc iml;resposta={[3] 'Sim', [3] 'Não', [3] 'Talvez'};resposta={[3] 1, [3] 2, [3] 0};print resposta;

RESPOSTA 1 1 1 2 2 2 0 0 0

Matrizes também podem ser obtidas como resultado de funções ou expressões, como porexemplo x=inv(a); em que os elementos de X são obtidos pela inversão de A. Há ainda umaopçã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 oincremento.

proc iml;r=1:5;s=10:5;t='xyz1':'xyz4';u=do(-2,1,1);v=do(-2,1,.5);print r s t u v;

R 1 2 3 4 5

S10 9 8 7 6 5

T xyz1 xyz2 xyz3 xyz4

U-2 -1 0 1

V-2 -1.5 -1 -0.5 0 0.5 1


Def in i ndo matr i ze sDef i n i ndo matr i ze sDef i n i ndo matr i ze sDef i n i ndo matr i ze s

Um comando também útil é o SHOW NAMES, que permite visualizar os atributos de todasas matrizes armazenadas. No caso, tem-se a dimensão, o tipo - numérica ou alfanumérica -e o número de bits para armazenamento.

proc iml;x={1 1 0, 1 1 0, 1 0 1, 1 0 1};y={1, 2, 3, 4};show names;print x y;

X 4 rows 3 cols num 8 Y 4 rows 1 col num 8 Number of symbols = 2 (includes those without values)

X Y 1 1 0 1 1 1 0 2 1 0 1 3 1 0 1 4

Matrizes especiais, como a matriz identidade, possuem comandos específicos para suaobtenção, sendo que o número entre parênteses indica sua dimensão:

proc iml;identid=i(5);print identid;

IDENTID 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á matrizes com elementos iguais, que podem ser iguais a um (matriz de uns) ou a qualqueroutro valor de interesse, são obtidas com a função J.



proc iml;uns1_5=j(1,5,1);uns5_5=j(5,5,1);quatros=j(3,4,4);print uns1_5 uns5_5 quatros;

UNS1_51 1 1 1 1

UNS5_5 QUATROS1 1 1 1 1 4 4 4 41 1 1 1 1 4 4 4 41 1 1 1 1 4 4 4 41 1 1 1 11 1 1 1 1

Para essa função, o primeiro valor dentro dos parênteses indica o número de linhas damatriz a ser criada, o segundo, o de colunas, e o terceiro, o valor que terão todos oselementos.

Outros comandos úteis são os que fornecem a dimensão das matrizes armazenadas, NROW(número de linhas) e NCOL (número de colunas):

proc iml;x={1 2 3, 4 5 6, [3] 9, [2] 1 0};m=nrow(x);n=ncol(x);print m n x;

M N X4 3 1 2 3 4 5 6 9 9 9 1 1 0

Finalizando esta sessão, mencionaremos ainda o comando STORAGE, que permitearmazenar (gravar) as matrizes numa biblioteca, permitindo seu uso posterior. Isso podeser útil nos casos em que os programas são extensos e/ou há limitação de memória docomputador. O comando FREE elimina as matrizes da memória RAM.

proc iml;


Operaçõe s Bás i ca sOperaçõ e s Bás i ca sOperaçõ e s Bás i ca sOperaçõ e s Bás i ca s

x={1 1 0, 1 1 0, 1 0 1, 1 0 1};y={1, 2, 3, 4};store x y;free x y;show storage;

Contents of storage library = WORK.IMLSTOR

Matrices: X Y

Modules:

3 Operações BásicasA partir desse ponto, em função de estarmos sempre apresentando comandos do SAS/IML,não apresentaremos mais os comandos PROC IML e QUIT, necessários para iniciar eencerrar a sessão.

3.1. Adição

A={2 2, 3 4};B={4 5, 1 0};SOMA=A+B;print A B soma;

A B SOMA2 2 4 5 6 73 4 1 0 4 4


Operaçõe s Bás i ca sOperaçõe s Bás i ca sOperaçõe s Bás i ca sOperaçõe s Bás i ca s

3.2. Subtração

DIF=A-B;print dif;

DIF-2 -3 2 4

3.3. Multiplicação por Escalar

C=4;AxC=A*C;print AxC;

AXC 8 812 16

3.4. Multiplicação de Matrizes

AxB=A*B;print AxB;

AXB10 1016 15

3.5. Soma Direta

Para tanto, deve-se utililizar a função BLOCK, que permite realizar a soma direta de umgrande número de matrizes.

a={2 2, 4 4};b={6 6, 8 8};c=block(a,b);



print a b c;

A B C2 2 6 6 2 2 0 04 4 8 8 4 4 0 0 0 0 6 6 0 0 8 8

3.6. Produto Direto (Kronecker)

a={1 0, 0 1};b={1, 2, 3};c1=a@b;c2=b@a;c3=a@a;print a b c1 c2 c3;

A B C1 C21 0 1 1 0 1 00 1 2 2 0 0 1 3 3 0 2 0 0 1 0 2 0 2 3 0 0 3 0 3C31 0 0 00 1 0 00 0 1 00 0 0 1

3.7. Produto de Hadamard

a={1 2, 3 4};b={4 8, 0 5};c=a#b;print c;



C4 160 20

3.8. Potência de Matriz Quadrada

É importante ressaltar que matrizes elevadas a números muito grandes podem apresentarcomo resultado valores com grandes imprecisões numéricas. A potência é obtida comosegue.

a={1 2, 3 4};a2=a**2;a3=a**3;inversa=a**(-1);print a2 a3 inversa;

A2 A3 INVERSA 7 10 37 54 -2 115 22 81 118 1.5 -0.5

3.9. Traço

x={1 2 3, 4 5 6, 6 2 4};traco=trace(x);print traco;

TRACO 10

3.10. Transposição

x={ 1.9 2 3, 40 5.2 6, 60 20 4.3};transp1=t(x);transp2=x`;



print transp1 transp2;

TRANSP1 TRANSP2 1.9 40 60 1.9 40 60 2 5.2 20 2 5.2 20 3 6 4.3 3 6 4.3

Note que há duas formas diferentes de transpor uma matriz, uma com o argumento pós-fixado e uma com o argumento pré-fixado. Evidentemente, ambas produzem o mesmoresultado.

3.11. Comparações entre Matrizes

O SAS/IML possui operadores que permitem realizar a comparação entre os elementos deduas matrizes, fornecendo como resultado uma nova matriz de uns e zeros, conforme acomparação seja verdadeira ou falsa, respectivamente. A sintaxe dos operadores é: matriz1< matriz2, matriz1 <= matriz2, matriz1 > matriz2, matriz1 >= matriz2, matriz1 =matriz2, matriz1 ^= matriz2. Os dois últimos operadores permitem identificar matriziguais, mas devemos sempre observar que isso pode ocorrer devido a erros de aproximação.Alguns exemplos:

a={1 7 3, 4 2 4};b={0 8 2, 4 1 3};menor= a<b;menorig= a<=b;maior= a>b;maiorig= a>=b;igual= a=b;diferent= a^=b;print menor menorig maior maiorig igual diferent;

MENOR MENORIG 0 1 0 0 1 0 0 0 0 1 0 0

MAIOR MAIORIG 1 0 1 1 0 1 0 1 1 1 1 1

IGUAL DIFERENT 0 0 0 1 1 1 1 0 0 0 1 1



3.12. Concatenação de Matrizes e Vetores

A concatenação pode ser feita horizontalmente ou verticalmente, fornecendo umaexcelente ferramenta para análise de dados em situações reais. Vejamos alguns exemplos:

A={4 3, -2 5};b={10, 8};C={7 -10, 4 13};A_b=A||b;A_C=A//C;print A_b A_C;

A_B A_C 4 3 10 4 3-2 5 8 -2 5 7 -10 4 13

A matriz A e o vetor b foram concatenados horizontalmente usando o operador ||. Já asmatrizes A e C foram concatenadas verticalmente usando //, obviamente observando se asdimensões são compatíveis.

3.13. Escalonamento de Matrizes

O SAS/IML possui duas funções para obtenção de formas escalonadas, HERMITE eECHELON, resultando matrizes na forma de Hermite e Echelon, respectivamente.

x={3 6 9, 1 2 5, 2 4 10};h=hermite(x);e=echelon(x);print h e;

H E1 2 0 1 2 00 0 0 0 0 10 0 1 0 0 0

Note que, para essas matrizes, as formas escalonadas resultantes são diferentes, sendo asegunda delas a forma escalonada canônica.

1 21 21 21 2 In trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze s


3.14. Posto de Matrizes

As formas escalonadas permitem determinar o posto de matrizes:

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=echelon(A);Besc=hermite(B);Cesc=echelon(C);print Aesc Besc Cesc;

AESC BESC1 0 1 1 0 00 1 -1 0 1 00 0 0 0 0 1

CESC1 2 3 40 0 0 0

Verifica-se que o posto de A, B e C é 2, 3 e 1, respectivamente.

Embora o SAS/IML não possua uma função que forneça diretamente o posto de uma matriz(a função RANK fornece o rank dos elementos, e não da matriz, como desejado), isso podeser conseguido com o auxílio de um programa auxiliar, apresentado a seguir. Ele bas

eia-se nos teoremas 55 (r[A]=r[A+A]) e 50 (A+A é simétrica e idempotente), lembrandoainda que para matrizes idempotentes o rank é igual ao traço. A função ROUND arredondao valor encontrado, que é sujeito a pequenos erros de aproximaçã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=round(trace(ginv(A)*A));RANK_B=round(trace(ginv(B)*B));RANK_C=round(trace(ginv(C)*C));


I n ve r sa de Matr i z Não Singu larIn ve r sa de Matr i z Não Singu larIn ve r sa de Matr i z Não Singu larIn ve r sa de Matr i z Não Singu lar

print RANK_A RANK_B RANK_C;

RANK_A RANK_B RANK_C 2 3 1

4 Inversa de Matriz Não SingularApresentaremos a seguir três formas diferentes de calcular a inversa de matriz nãosingular, sendo que a primeira delas (INV1) foi incluída apenas por razões didáticas, usandoo algoritmo de Gauss.

B={3 1 2, 1 4 1, 2 1 5};INV1=hermite(B||I(3));INV2=inv(B);INV3=B**(-1);print INV1 INV2 INV3;

INV11 0 0 0.475 -0.075 -0.1750 1 0 -0.075 0.275 -0.0250 0 1 -0.175 -0.025 0.275

INV2 INV3 0.475 -0.075 -0.175 0.475 -0.075 -0.175-0.075 0.275 -0.025 -0.075 0.275 -0.025-0.175 -0.025 0.275 -0.175 -0.025 0.275

Aproveitaremos a matriz INV1, que contém a inversa de B à direita da identidade, paramostrarmos como selecionar submatrizes de uma dada matriz. É possível selecionar linhas,colunas ou até mesmo elementos de interesse. Para tanto, utilizam-se subscritos,que sãooperadores localizados dentros de colchetes, após o operando. O primeiro número refere-sea(s) linha(s) a ser(em) retirada(s) e o segundo a(s) coluna(s), sendo que eles devem serseparados por uma vírgula. Para selecionar um único elemento, como por exemplo da 3alinha e 4a coluna (lembre-se que INV1 tem dimensão 3 x 6):

B={3 1 2, 1 4 1,

1 4141414 Introdução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze sI n trodução ao SAS/IML para Álgebra de Matr i ze s

I n ver sa de Matr i z Não Singu larIn ver sa de Matr i z Não Singu larIn ver sa de Matr i z Não Singu larIn ver sa de Matr i z Não Singu lar

2 1 5};INV1=hermite(B||I(3));elem3_4=INV1[3,4];print elem3_4;

ELEM3_4 -0.175

Para selecionar a linha 3 ou a coluna 4:

B={3 1 2, 1 4 1, 2 1 5};INV1=hermite(B||I(3));linha3=INV1[3,];coluna4=INV1[,4];print linha3 coluna4;

LINHA3 COLUNA40 0 1 -0.175 -0.025 0.275 0.475 -0.075 -0.175

Para selecionar uma submatriz, como por exemplo a inversa sem a parte referente à matrizidentidade:

B={3 1 2, 1 4 1, 2 1 5};INV1=hermite(B||I(3));INVERSA=INV1[, {4 5 6}];print INVERSA;

INVERSA 0.475 -0.075 -0.175-0.075 0.275 -0.025-0.175 -0.025 0.275

Note que é necessário utilizar chaves para correta especificação da sub-matriz.


Determ inante sDeterm inante sDeterm inante sDeterm inante s

5 DeterminantesO SAS/IML calcula o determinante de matrizes quadradas de forma simples, através dafunção DET, considerando o fato de que o determinante pode ser obtido através do produtodos autovalores da matriz em questão, que por sua vez são obtidos pela decomposição LU(matrizes triangulares). O resultado fornecido é uma matriz 1x1 (escalar) contendo odeterminante.

proc iml;a={1 1 1, 1 2 4, 1 3 9};determ=det(a);print determ;

DETERM 2

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

6 Sistemas Consistentes de Equações LinearesO SAS possui a função SOLVE para a solução de sistemas de equações lineares (consistentese determinados), que é muito eficiente em termos computacionais, devendo ser usado emsituações práticas. Contudo, sempre com finalidades didáticas, serão apresentadas outrasformas de resolver sistemas.

Para usar a função SOLVE para resolver o sistema Ax=b, pode-se proceder comoapresentado a seguir.

A={4 3, -2 5};b={10, 8};x=solve(A,b);print x;

X


Sistemas Cons i s t ente s de Equaçõe s L ineare sSi st emas Cons i s t ente s de Equaçõe s L ineare sSi st emas Cons i s t ente s de Equaçõe s L ineare sSi st emas Cons i s t ente s de Equaçõe s L ineare s

1 2

As formas escalonadas apresentadas também podem ser usadas para solução de sistemasde equações lineares. Considerando o mesmo sistema com duas equações e duasincógnitas, escrito na forma matricial Ax=b, um programa para obter sua solução pode ser:

A={ 4 3, -2 5};b={10, 8};A_b=A||b;SOL1=echelon(A_b);SOL2=hermite(A_b);print SOL1 SOL2;

O escalonamento forneceu a solução do sistema para x, com resultados 1 e 2.

SOL1 SOL2 1 0 1 1 0 1 0 1 2 0 1 2

Caso desejemos resolver o sistema usando o método de Kramer, podemos usar o seguinteprograma:

reset print;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 maissintéticos e eficientes. Nesse contexto, foram definidas as matrizes A, B, A1, A2 e A3, com


Sistemas Cons i s tente s de Equaçõe s Lineare sSi st emas Cons i s tente s de Equaçõe s Lineare sSi st emas Cons i s tente s de Equaçõe s Lineare sSi st emas Cons i s tente s de Equaçõe s Lineare s

seus respectivos determinantes. Em seguida, foi obtida a solução do sistema, concatenadano vetor SOLUCAO. O comando RESET PRINT solicita a visualização de todas as matrizescriadas, facilitando a visualização dos resultados, embora com excesso de informações. Paradesativá-lo, deve-se usar RESET NOPRINT. A saída do programa acima é

A 3 rows 3 cols (numeric) 1 1 2 2 0 3 1 1 0

B 3 rows 1 col (numeric) 9 11 3

A1 3 rows 3 cols (numeric) 9 1 2 11 0 3 3 1 0



DET_A 1 row 1 col (numeric) 4

DET_A1 1 row 1 col (numeric) 4



X1 1 row 1 col (numeric) 1




Sistemas Cons i s t ente s de Equaçõe s L ineare sSi st emas Cons i s t ente s de Equaçõe s L ineare sSi st emas Cons i s t ente s de Equaçõe s L ineare sSi st emas Cons i s t ente s de Equaçõe s L ineare s

SOLUCAO 3 rows 1 col (numeric) 1 2 3

Um programa mais sintético pode ser:

reset noprint;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);print solucao;

SOLUCAO 1 2 3

No caso de obter soluções para sistemas de equações lineares homogêneos do tipo Ax=(,sendo A de dimensão m x n, mn, é possível utilizar a função HOMOGEN, que retornasoluções como colunas ortonormais de uma matriz X com as propriedades Ax=( (P1) eX'X=I (P2). Para tanto, o posto de A deve ser r<n, ou seja, A deve ser de posto incompleto.A dimensão dessa matriz X é n x (n-r), e portanto deve ser dada especial atenção a matrizesmal-condicionadas, já que o correto posto de A pode não ser obtido em função de erros deaproximação.

A={4 2 2, 2 2 0, 2 0 2};X=HOMOGEN(a);P1=A*X;P2=T(X)*X;print X P1 P2;

X P1 P2-0.57735 2.22E-16 1 0.5773503 2.22E-16 0.5773503 0

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


Raíze s Caracter í s t i ca s (aut ova l o re s )Ra í ze s Caracter í s t i ca s (aut ova l o re s )Ra í ze s Caracter í s t i ca s (aut ova l o re s )Ra í ze s Caracter í s t i ca s (aut ova l o re s )

Um outro fato que pode ser mencionado é que a dimensão da matriz com as soluções dosistema homogêneo podem ser usadas para calcular o posto da matriz A, já que não há umafunção específica para tanto:

POSTO_A=NROW(X)-NCOL(X);print POSTO_A;

POSTO_A 2

7 Raízes Características (autovalores)Para obter os autovalores ("eingenvalues") de matrizes quadradas simétricas, usa-se afunção EIGVAL, que fornece um vetor com os autovalores ordenados em ordemdecrescente.

A={3 1 1, 1 3 1, 1 1 3};AUTOVAL=EIGVAL(A);print AUTOVAL;

AUTOVAL 5 2 2

Caso necessitemos obter os autovalores a partir da equação característica, podemos usar afunçã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 umpolinômio de grau n-1, sendo que eles devem estar em ordem decrescente, com relação aspotências do polinômio. O resultado é uma matriz de dimensão (n-1) x 2, com as raízes dopolinômio; a primeira coluna contém a parte real e a segunda, a parte imaginária (quedeverá ser zero, no contexto desse curso). No exemplo em questão, a equação característicaé -x3+9x2-24x+20=0, e assim

p={-1 9 -24 20};raizes=POLYROOT(p);


Raíze s Caracter í s t i ca s (au to va lo re s)Ra í ze s Caracter í s t i ca s (au to va lo re s)Ra í ze s Caracter í s t i ca s (au to va lo re s)Ra í ze s Caracter í s t i ca s (au to va lo re s)

print raizes;

RAIZES 2.0000003 0 1.9999997 0 5 0

É interessante notar que as soluções estão sujeitas a erros de aproximação, em função doalgoritmo empregado. Um outro aspecto que deve ser lembrado diz respeito ao fato do SAS/IML usar um algoritmo que não garante que todas as raízes do polinômio sejamencontradas; nesses casos, há uma mensagem para o usuário na janela LOG. Nas situaçõesem que desejamos obter matrizes diagonais com as raízes características de uma dadamatriz, combinamos a função DIAG, que cria matrizes diagonais, com a função EIGVAL.

M={ 2 -1 -1, -1 2 -1, -1 -1 2};AUTOV=EIGVAL(M);L=DIAG(AUTOV);PRINT AUTOV L;

AUTOV L 3 3 0 0 3 0 3 0-9.11E-17 0 0 -9.11E-17

É fácil verificar que, exceto devido a erros de aproximação, M é positiva semi-definida.

Caso apliquemos o comando DIAG em uma matriz (e não num vetor, como no exemploacima), teremos como resultado uma matriz formada pelos elementos da diagonal dessamatriz, sendo os demais substituídos por zero.

Para facilitar a visualização, poderíamos ainda usar a opção FUZZ, associada ao comandoRESET, de modo que os números muito pequenos (por volta de 1E-12) sejam impressoscomo zero.

reset fuzz;PRINT AUTOV L;

AUTOV L 3 3 0 0 3 0 3 0


Vetore s Caracter í s t i c o s (aut ovet or e s )Vetore s Caracter í s t i c o s (aut ovet or e s )Vetore s Caracter í s t i c o s (aut ovet or e s )Vetore s Caracter í s t i c o s (aut ovet or e s )

0 0 0 0

Para voltarmos ao default, dentro da mesma seção do IML, devemos usar o comando RESETNOFUZZ.

8 Vetores Característicos (autovetores)A obtenção dos autovetores ("eingenvectors") de matrizes quadradas simétricas pode serfeita pelo uso da função EIGVEC, que cria uma matriz cujas colunas são autovetoresortonormais.

A={3 1 1, 1 3 1, 1 1 3};AUTOVEC=EIGVEC(A);print AUTOVEC;

AUTOVEC0.5773503 -0.408248 0.70710680.5773503 -0.408248 -0.7071070.5773503 0.8164966 0

É fácil verificar que esses vetores são ortogonais entre si, pois o produto entre autovaloresdiferentes é 0, apresentando como resultado de seu produto consigo mesmo o valor 1(norma 1). Em outras palavras, o produto entre a transposta da matriz com os autovalores eela própria é a identidade:

RESET FUZZ;AUTOVEC=EIGVEC(A);ORTONORM=AUTOVEC`*AUTOVEC;print ORTONORM;

ORTONORM1 0 00 1 00 0 1

Em situações práticas, podemos obter os autovalores e autovetores diretamente, com o usoda sub-rotina EIGEN, acionada pelo comando CALL, que invoca determinadas estruturas


Formas Equ iva l en te sFormas Equ iva l en te sFormas Equ iva l en te sFormas Equ iva l en te s

(operações) pré-armazenadas. O primeiro termo dentro dos parênteses especifica o nomeda matriz que conterá os autovalores; o segundo, da matriz que conterá os autovetores;oterceiro indica sobre qual matriz o comando será aplicado.

CALL EIGEN(AUTOVAL,AUTOVET,A);print AUTOVAL AUTOVET;

AUTOVAL AUTOVET 5 0.5773503 -0.408248 0.7071068 2 0.5773503 -0.408248 -0.707107 2 0.5773503 0.8164966 0

9 Formas EquivalentesAs matrizes A e B são equivalentes se existem M e N não singulares tais que MAN=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 é triangularsuperior e idempotente. Se A for não singular, o resultado será uma matriz identidade.

A={4 2 2, 2 2 0, 2 0 2};B=hermite(A);print B;

B1 0 10 1 -10 0 0

É possível também obter a forma normal de Echelon que, nesse caso, coincide com a formanormal de Hermite:

B=echelon(A);print B;

B


Formas Equ i va l en te sFormas Equ i va l en te sFormas Equ i va l en te sFormas Equ i va l en te s

1 0 10 1 -10 0 0

Se a matriz em questão for quadrada, as formas normais de Hermite e Echelon podem serobtidas uma a partir da outra, apenas rearranjando as linhas de zeros. A forma de Echelonsempre fornece uma forma escalonada canônica.

Para a obtenção de uma matriz D diagonal congruente (e, portanto, equivalente) à matriz Areal e simétrica considerada, tal que TAT '=D, utiliza-se a sub-rotina EIGEN, pelo comandoCALL EIGEN(D, T, A), já descrito anteriormente. Para tanto, utilizam-se as seguintespropriedades das matrizes obtidas com esse comando: T' = T-1 (ortogonal) eA=T*diag(D)*T', ou seja, diag(D) é uma matriz diagonal congruente à matriz A.

A={4 2 2, 2 2 0, 2 0 2};CALL EIGEN (D,T,A);CONG=DIAG(D);A2=T*DIAG(D)*T`;reset FUZZ;print A A2;

A A24 2 2 4 2 22 2 0 2 2 02 0 2 2 0 2

Dessas propriedades, conseguimos também verificar que diag(D)= T'*A*T, em que diag(D)é uma matriz diagonal que exibe as raízes características de A, e T é uma matriz ortogonalcujas linhas são vetores característicos normalizados associados as raízes características deA.

L1=DIAG(D);L2=T`*A*T;print L1 L2;

L1 L26 0 0 6 0 00 2 0 0 2 00 0 0 0 0 0

Nesse caso, A e diag(D) são também similares.


FatoraçãoFatoraçãoFatoraçãoFatoração

10 FatoraçãoO SAS/IML não realiza diretamente fatorações do tipo A=RR', fatorações de posto completodo tipo A=BC e operações de diagonalização simultânea. No primeiro caso, é possível aobtenção de fatorações do tipo raiz quadrada de outra forma, apresentada abaixo.

Uma vez que o principal interesse na fatoração de posto completo é para a obtenção dainversa generalizada de Moore-Penrose, a ausência da fatoração não representa empecilho,uma vez que essa inversa generalizada é calculada pelo SAS (itens seguintes). No manual,temos diversas outras decomposições que podem interessar ao leitor.

Fatorações do tipo A=U'U são obtidas com a função ROOT, que retorna uma matriztriangular superior U que satisfaz a igualdade, sendo que a matriz A deve ser simétrica epositiva definida. Essa decomposição é conhecida como decomposição de Cholesky, e amatriz U é dita raiz quadrada de A.

U=ROOT(A);A1=U`*U;print U A1;

U A12 1 1 4 2 20 1 -1 2 2 00 0 0 2 0 2

Vale ressaltar que alguns autores denominam a fatoração A=RR' como raiz quadrada, e nãoA=U'U, como o SAS considera. É interessante notar que essa fatoração também é do tipotriangular (A=LU), uma vez que U' é triangular inferior e U é triangular superior.

Para realizar decomposição por valores singulares do tipo X=USV´, o SAS/IML possui a sub-rotina SVD, acessada pelo comando CALL. Sua sintaxe é CALL SVD(U, S, V, X), em que Xrepresenta a matriz a ser decomposta, U, S e V, o nome das matrizes a serem criadas com oresultado da decomposição, de forma que X=U*diag(S)*V´ e S contenha os valoressingulares de A. Note que o algoritmo apresentado no texto das aulas teóricas produzmatrizes diferentes, mostrando que a decomposição não é única.

reset print;X={1 1 0, 1 1 0, 1 0 1, 1 0 1};call svd(U,S,V,X);X1=U*diag(S)*V`;


I n ve r sa s Genera l i zada sIn ve r sa s Genera l i zada sIn ve r sa s Genera l i zada sIn ve r sa s Genera l i zada s

X 4 rows 3 cols (numeric)1 1 01 1 01 0 11 0 1

U 4 rows 3 cols (numeric)0.5 0.5 0.70710680.5 0.5 -0.7071070.5 -0.5 00.5 -0.5 0

S 3 rows 1 col (numeric)2.44948971.4142136 0

V 3 rows 3 cols (numeric)0.8164966 0 -0.577350.4082483 0.7071068 0.57735030.4082483 -0.707107 0.5773503

X1 4 rows 3 cols (numeric)1 1 01 1 01 0 11 0 1

11 Inversas GeneralizadasÉ possível calcular a inversa generalizada de Moore-Penrose usando o SAS/IML, com ocomando GINV(X), em que X é a matriz a ser invertida. A inversa obtida satisfaz as quatrocondições da inversa, conforme apresentado abaixo.

X={1 1 0, 1 1 0, 1 0 1, 1 0 1};MP1=ginv(X);P1=MP1*X*MP1;P2=X*MP1*X;P3a=X*MP1; P3b=t(X*MP1);P4a=MP1*X; P4b=t(MP1*X);print MP1;


I nv er sas Genera l i zada sInv er sas Genera l i zada sInv er sas Genera l i zada sInv er sas Genera l i zada s

print P1 P2 P3a P3b P4a P4b;

MP10.1666667 0.1666667 0.1666667 0.16666670.3333333 0.3333333 -0.166667 -0.166667-0.166667 -0.166667 0.3333333 0.3333333

P1 P20.1666667 0.1666667 0.1666667 0.1666667 1 1 00.3333333 0.3333333 -0.166667 -0.166667 1 1 0-0.166667 -0.166667 0.3333333 0.3333333 1 0 1 1 0 1P3A P3B0.5 0.5 0 0 0.5 0.5 0 00.5 0.5 0 0 0.5 0.5 0 0 0 0 0.5 0.5 0 0 0.5 0.5 0 0 0.5 0.5 0 0 0.5 0.5

P4A P4B0.6666667 0.3333333 0.3333333 0.6666667 0.3333333 0.33333330.3333333 0.6666667 -0.333333 0.3333333 0.6666667 -0.3333330.3333333 -0.333333 0.6666667 0.3333333 -0.333333 0.6666667

As inversas generalizadas reflexivas g2 para matrizes de delineamento podem ser obtidascom a função SWEEP. Tomando o modelo linear de Gauss-Markov , podemosobter o sistema

e, usando o operador SWEEP sobre esse sistema, calcular uma matriz

A partição superior esquerda dessa matriz contêm a inversa g2 de X'X, a partição superiordireita as soluções para B e apartição inferior direita a soma de quadrados do resíduo. Esseé o procedimento usado pelo PROC GLM do SAS para realizar análises de variância. Umexemplo é apresentado abaixo, sendo que o segundo argumento da função SWEEP é umvetor com índices variando de 1 até o posto de X.

reset print;X={1 1 0, 1 1 0, 1 0 1, 1 0 1};y={1,

y XB e+=

X'X X'yy 'x y 'y

X'X( )g2 X'X( )g2X'y

y 'X X'X( )g2– y ' I X X'X( )g2X '–( )y


I n ve r sa s Genera l i zada sIn ve r sa s Genera l i zada sIn ve r sa s Genera l i zada sIn ve r sa s Genera l i zada s

2, 4, 3};A=(t(X)*X||t(X)*y)//(t(y)*X||t(y)*y);t=1:ncol(X);G=sweep(A,t);reset noprint;

X 4 rows 3 cols (numeric)

1 1 01 1 01 0 11 0 1

Y 4 rows 1 col (numeric)

1243

A 4 rows 4 cols (numeric)

4 2 2 10 2 2 0 3 2 0 2 710 3 7 30

T 1 row 3 cols (numeric)

1 2 3

G 4 rows 4 cols (numeric)

0.5 -0.5 0 3.5-0.5 1 0 -2 0 0 0 0-3.5 2 0 1

Leitores já familiarizados com o PROC GLM podem verificar as afirmativas acima com oseguinte programa:

data;input trat Y@@;datalines;1 1 1 2 2 4 2 3;


I nv er sas Genera l i zada sInv er sas Genera l i zada sInv er sas Genera l i zada sInv er sas Genera l i zada s

proc glm; class trat; model Y=trat/inverse;run;

The GLM Procedure Class Level Information Class Levels Values trat 2 1 2

Number of observations 4

X'X Generalized Inverse (g2)

Intercept trat 1 trat 2 YIntercept 0.5 -0.5 0 3.5trat 1 -0.5 1 0 -2trat 2 0 0 0 0Y 3.5 -2 0 1

Dependent Variable: Y

Sum ofSource DF Squares Mean Square F Value Pr > F

Model 1 4.00000000 4.00000000 8.00 0.1056Error 2 1.00000000 0.50000000Corrected Total 3 5.00000000

R-Square Coeff Var Root MSE Y Mean0.800000 28.28427 0.707107 2.500000

Source DF Type I SS Mean Square F Value Pr > F

trat 1 4.00000000 4.00000000 8.00 0.1056

Source DF Type III SS Mean Square F Value Pr > F

trat 1 4.00000000 4.00000000 8.00 0.1056


Documents

Apostila Proc Iml