Escalonando matrizes em APL

Esta função gera uma matriz-identidade de ordem N. Matrizes-identidade são matrizes quadradas (mesmo número de linhas e colunas) que possuem apenas o número 1 em sua diagonal principal, e 0 em todas as outras posições.

Note que as funções de APL podem possuir uma variável de retorno (normalmente nomeio-a R). No caso de IDENTITY, trata-se de uma função de um parâmetro, que nomeei N.

Veja o código, levando em consideração que APL interpreta cada linha de expressão da direita para a esquerda, diferente do convencional. Esta regra evita alguns problemas de ambiguidade sintática, juntamente com outras que não mencionarei aqui.

∇R←IDENTITY N
  R←(N N)⍴(1,N⍴0)
∇

IDENTITY gera uma matriz N×N, e escreve em cada linha desta matriz a sequência de um número 1 e N números 0. Caso alguma sequência de números "sobre", ela será escrita no início da próxima linha. E quando não houver mais números 1 ou 0 a serem escritos, essa sequência de 1's e 0's se repete. Veja:

      1,3⍴0
1 0 0 0
      3 3⍴(1,3⍴0)
1 0 0
0 1 0
0 0 1

Esta função adiciona uma nova linha com zeros na parte inferior da matriz. Esta função independe do tamanho da matriz, que é deduzido durante a operação.

∇R←EXTENDBOTTOM MAT
  R←(((↑⍴MAT)⍴1),0)\[1]MAT
∇

Veja o que acontece quando aplicamos ela a uma matriz-identidade de dimensões 3×3:

      EXTENDBOTTOM IDENTITY 3
1 0 0
0 1 0
0 0 1
0 0 0

Esta função não possui uso na implementação do escalonamento em si, mas serve para ajudar a gerar algumas matrizes de teste.

Similar a EXTENDBOTTOM, esta função adiciona uma nova coluna de zeros à direita da matriz.

∇R←EXTENDRIGHT MAT
  R←(((2⌷⍴MAT)⍴1),0)\MAT
∇

Aplicando esta operação a uma matriz-identidade 3×3, temos:

      EXTENDRIGHT IDENTITY 3
1 0 0 0
0 1 0 0
0 0 1 0

Também utilizaremos esta função para criar matrizes para testes.

Agora, faremos uma função diádica que coleta o pivô de uma matriz MAT, na linha N. Se a linha não possui pivô, retornamos uma resposta vazia.

Antes de mais nada, vejamos o que é um pivô de uma certa linha. Suponha a seguinte matriz:

O pivô de uma linha da matriz é o primeiro elemento da linha que, quando procuramos da esquerda para direita nela, não seja nulo. Sendo assim, os pivôs de cada linha da matriz serão 5, 2 e 4, respectivamente.

Vejamos como isto pode ser feito em APL. Como já dito, a função PIVOT é uma função diádica (seus parâmetros são N e MAT), que possui os valores locais LST e IDX (que precisam ser indicados no cabeçalho da função, ou serão redeclarados como variáveis globais).

∇R←N PIVOT MAT;LST;IDX
  LST←N⌷[1]MAT
  IDX←(0≠LST)⍳1
  →(IDX>⍴LST)/0
  R←IDX⌷LST
∇

Aqui temos a variável LST armazenando temporariamente a linha onde precisamos encontrar o pivô. Por exemplo, podemos ter LST valendo 0 4 9.

A variável IDX será um valor que armazena o índice do primeiro elemento não-nulo de LST. Por exemplo, para o valor de LST que estamos supondo, IDX valerá 2 (segunda posição na linha).

Se IDX for um valor maior do que o número de elementos naquela linha, então nenhum pivô foi encontrado, e a linha é nula; neste caso, encerramos a função, sem dar uma resposta.

Caso tudo tenha corrido bem, a variável de retorno recebe o elemento de LST que está na posição IDX.

Esta é a primeira função no nosso código que realmente tem alguma espécie de estrutura de repetição.

Esta operação é o algoritmo trivial de Euclides para calcular o máximo divisor comum entre dois números. Convencionaremos que A seja sempre maior que B, para que este algoritmo funcione adequadamente.

Vou me abster de explicá-la totalmente, porque esta função é simples e óbvia (alguns autores dizem que é o algoritmo mais velho do mundo!).

∇R←A GCD B
  →(B=0)/END
  (A B)←(B (|B | A))
  →1
END: R←A
∇

Veja que, sendo uma função diádica, GCD precisa ser escrita entre ambos os números:

      99 GCD 63
9

Esta função monádica simplifica um certo vetor de elementos, dividindo todos pelo máximo divisor comum da porção inteira de todos os números.

Note que o vetor só será realmente simplificado se nenhum dos elementos tornar-se um número fracionário. Decidi fazer isso porque não queria nenhum número "quebrado" em meus resultados.

∇R←SIMPVEC VEC;SIMPL
  SIMPL←VEC÷GCD/⌊|VEC[⍒|VEC]
  →(×+/0<SIMPL-⌊SIMPL)/DONTSIMPL
  R←SIMPL
  →0
  DONTSIMPL:R←VEC
∇

Vejamos a aplicação desta função:

      SIMPVEC 4 2 8
2 1 4

A função GAUSSELIM é uma função diádica, que fixa uma linha de índice L na matriz A, e executa operações nas outras linhas, de forma que o pivô da linha L seja o único elemento não-nulo de sua coluna.

Vejamos o que isto significa com notação matemática. Vamos ver a matriz do problema das maçãs:

Suponhamos que desejamos fixar a primeira linha. Sendo assim, temos L = 1.

O pivô da linha 1 é o número 1. Portanto, os números que estão na mesma coluna (6 e 5) precisam ser zerados.

Para efeito de exemplo, vamos zerar a linha 2. Para tanto, realizaremos o seguinte algoritmo:

• Multiplique os elementos da linha 2 pelo pivô da linha fixada;
• Multiplique os elementos da linha fixada pelo pivô da linha 2;
• Subtraia as linhas geradas, armazene na linha não-fixada (2).

E o mesmo processo seria feito para a linha 3.

Em APL, o primeiro a ser feito é simplificar a linha fixada. Para isso, aplicamos a função SIMPVEC na linha L. Este não é um passo obrigatório, mas garante números menores no processo.

Depois, geramos uma lista de números de linhas sequenciais. Para a nossa matriz das maçãs, teríamos uma lista assim:

      ⍳3
1 2 3

Desta lista gerada, removemos a linha fixada (L = 1), e armazenamos o resultado na variável local ZEROLNS.

Agora, para cada índice de linha em ZEROLNS, aplicamos o processo descrito anteriormente. Isto é feito com uma função anônima declarada inline. O operador Each (¨) garante que esta função seja aplicada a cada elemento de ZEROLNS, em sequência.

A nova matriz é então armazenada na variável de resultado.

∇R←L GAUSSELIM A;ZEROLNS
  →(0=+/|A[L;])/END
  A[L;]←(×L PIVOT A)×SIMPVEC A[L;]
  ZEROLNS←(~L=ZEROLNS)/ZEROLNS←⍳↑⍴A
  R←{A[⍵;]←((L PIVOT A)×A[⍵;])-(A[⍵;L]×A[L;])}¨ZEROLNS
END:R←A
∇

Aplicando eliminação gaussiana, ao tomarmos a primeira linha do problema como fixa, teremos:

      1 GAUSSELIM A
1   3   4   53
0 ¯11 ¯22 ¯209
0 ¯14 ¯17 ¯189

A função GAUSSJORDAN é uma função monádica que realiza o escalonamento completo de uma matriz. Seu algoritmo é simples: para cada linha da matriz, realize uma eliminação gaussiana, tomando-a como fixa.

No final do processo, tudo o que temos que fazer é realizar uma simplificação extra em todas as linhas, e então repassar a nova matriz ao valor de retorno.

∇R←GAUSSJORDAN MAT
  R←{MAT←⍵ GAUSSELIM MAT}¨⍳↑⍴MAT
  R←{MAT[⍵;]←SIMPVEC MAT[⍵;]}¨⍳↑⍴MAT
  R←MAT
∇

Este último passo é, na realidade, redundante. APL possui uma função Dominó (⌹) que, quando usada de forma monádica, calcula a inversa de uma certa matriz. Compare estes dois resultados:

      ⌹F
¯1.5258556E¯16  1 ¯1
 2.0000000E0   ¯2 ¯1
¯1.0000000E0    1  1
      ¯3↑[2]GAUSSJORDAN F,IDENTITY 3
 0  1 ¯1
 2 ¯2 ¯1
¯1  1  1

Veja que o resultado dado pela função Dominó (⌹) não é exatamente o mesmo, porque ela utiliza outros métodos e trabalha com pontos flutuantes. Porém, podemos observar que o valor ¯1.5258556E¯16, que corresponde aproximadamente a -0.00000000000000015258556, é um valor bem próximo do requerido (0). Este tipo de ruído relacionado a pontos flutuantes é esperado, e portanto pode ser mais interessante utilizar o Dominó (⌹) que programar todo o processo de Gauss-Jordan.