Primeira vez aqui? Seja bem vindo e cheque o FAQ!
x

Precisão Simples ou Dupla: Implicações na Solução de Sistemas Lineares

0 votos
37 visitas
perguntada Out 26 em Ciência da Computação por Carlos Piccioni (21 pontos)  
editado Out 26 por Carlos Piccioni

O problema descrito a seguir encontra-se no livro Scientific Computing: An Introductory Survey, de Michael T. Heath, segunda edição (2018), capítulo 2, seção Computer problems, página 100, Problema 2.5:

  • item a)

Use uma rotina de precisão simples de Eliminação Gaussiana para resolver o sistema \(Ax=b\), em que:

\[ A= \begin{bmatrix} 21.0 & 67.0 & 88.0 & 73.0 \\ 76.0 & 63.0 & 7.0 & 20.0 \\ 0.0 & 85.0 & 56.0 & 54.0 \\ 19.3 & 43.0 & 30.2 & 29.4 \end{bmatrix} \]
\[ b= \begin{bmatrix} 141.0 \\ 109.0 \\ 218.0 \\ 93.7 \end{bmatrix} \]

  • item b)

Calcule o residual \({r=b-Ax}\) usando aritmética de precisão dupla, se disponível (mas armazenando o resultado final em um vetor \({r}\) de precisão simples). (Se apenas um tipo de precisão estiver disponível em seu ambiente computacional, então faça todo o problema utilizando o tipo disponível.)

  • item c)

Resolva o sistema linear \({Az=r}\) para obter a solução "aprimorada", \({x+z}\). Note que \({A}\) não precisa ser refatorada.

  • item d)

Repita os passos b e c até que nenhuma melhora no resultado seja observada.

Compartilhe

2 Respostas

0 votos
respondida Out 26 por Carlos Piccioni (21 pontos)  
editado Nov 11 por Carlos Piccioni

Item a)

Inicialmente, é interessante apresentar os conceitos de precisão simples (ou única) e precisão dupla (ou single precision e double precision). Para isso, primeiro devemos responder a questão: como computadores armazenam números reais? Sabemos que toda informação é armazenada de forma binária: zeros ou uns. Então, a princípio, uma simples conversão de base poderia parecer adequada. Mas como representar números irracionais, por exemplo? Ou frações exatas na base decimal, mas que geram sequências infinitas na base binária? (como será visto mais a frente). Como padronizar a quantidade de memória que será alocada para armazenar cada número? (além de outras questões, como otimizar as operações aritméticas realizadas pelo processador.)

Como solução, adota-se o chamado sistema de ponto flutuante (floating-point system). Nele, cada número real \(N\) é representado por uma mantissa (\(m\)), um expoente (\(e\)), uma base (\(b\)), e (\(s\)), que define o sinal no número real, tal que:

\[ N=(-1)^s \times m \times b^e \]

Em sistemas digitais, \(b=2\), e a mantissa e o expoente são obviamente armazenados na forma binária. Por exemplo, para representarmos o número real 5.75:

  • Conversão da base decimal para a base binária: \(5.75_{10} = 1\times2^2+0\times2^1+1\times2^0+1\times2^{-1}+1\times2^{-2} = 101.11_{2}\)

  • "Normalizando": \(101.11_{2} = 1.0111_{2}\times2^2\)

  • Ou seja, uma possível representação de \(5.75_{10}\) em um sistema de ponto flutuante binário é:

\[ 5.75_{10} = (-1)^0\times1.0111_{2}\times2^{{10}_2} \]

Em que \(s = 0_{2}, m = 1.0111_{2}, e = 10_{2} (=2_{10})\).

E em casos de números com uma quantidade bem maior de dígitos, ou mesmo infinita, como definir quantos bits devemos alocar para cada um desses elementos?

Aqui entra o conceito de precisão simples e precisão dupla. O Institute of Electrical and Electronics Engineers, por meio da norma IEEE 754, definiu requisitos mínimos para esses dois formatos, sendo que o formato de precisão simples utiliza 32 bits para armazenar um número real, enquanto o formato de precisão dupla, 64 bits. E a quantidade de bits para cada elemento da representação de ponto flutuante é divida da seguinte forma em cada padrão:

\[ \begin{array}{r|ccc|c} Precisão & Sinal (s) & Expoente (e) & Mantissa (m) & Total(bits)\\ \hline Simples & 1 & 8 & 23 & 32\\ Dupla & 1 & 11 & 52 & 64 \end{array} \]

Existem outros detalhes da implementação desses formatos, como rotação do expoente para permitir expoentes negativos, e omissão do primeiro bit da mantissa. Esses detalhes podem ser conferidos nas referências ao final dessa resposta. Também é importante salientar que a IEEE 754 define requisitos mínimos, logo algumas implementações podem, internamente, utilizar precisões maiores, mas armazenar os resultados finais na precisão indicada.

Como curiosidade, no caso do formato de precisão simples, você pode testar como funciona a conversão de números reais para a representação de ponto flutuante de forma rápida em:

www.h-schmidt.net/FloatConverter/IEEE754.html

Note, por exemplo, que não é possível armazenar de forma precisa um valor real como \(0.1_{10}\) (base decimal) em uma representação de base binária, incluindo a de ponto flutuante. No caso de precisão simples, o valor real armazenado será:

0.100000001490116119384765625

ou seja, com erro na ordem de \(10^{-9}\). Isso ocorre pois \(0.1_{10}\) é representado por uma sequência infinita na forma \(0.00011001100110011...\) na base binária. Como temos uma limitação de bits para representar o número, haverá esse erro devido ao truncamento do valor até o limite de bits disponível. O problema é o mesmo da fração \(1/3=0.333...\) na base decimal: precisamos de uma sequência infinita de \(3\)'s para representá-la.

E \(0.1_{10}\), em precisão dupla, na realidade será armazenado como:

0.1000000000000000055511151231257827021181583404541015625

ou seja, um erro na ordem de \(10^{-18}\). Para verificar esse valor, no Python, basta executar:

from decimal import Decimal
Decimal.from_float(0.1)

Mas então por que, no Python, ao realizarmos uma operação como \(1/10\) no console, o resultado padrão exibido é exatamente 0.1? Pois a exibição do resultado geralmente é arredondada até certas casas decimais. Mas, para verificar, por exemplo, o problema que apontamos acima com relação as representações binárias, tente executar \(.1 + .1 + .1 == .3\) em um console Python, ou verifique o resultado da operação \(0.8 * 12\) no console padrão do Spyder, e surpreenda-se com a resposta...

Apesar do erro ser pequeno nessas representações, no caso de um grande número de operações complexas ou consecutivas, estes erros podem ser propagados e amplificados, gerando, em alguns casos, resultados finais com imprecisões em um nível indesejado. E assim, sem entrar em mais detalhes sobre as representações e operações aritméticas envolvendo ponto flutuante, já pode-se notar que a escolha de um ou de outro formato, simples ou duplo, pode determinar a precisão do resultado de operações envolvendo números reais.

Para o problema em questão, se resolvido por qualquer técnica não computacional em que carregamos as frações, quando existentes, até o final, encontramos como solução o vetor \({x} = [-1, 2, -3, 4]\).

Porém, para explorar as questões envolvendo precisão de ponto flutuante, utilizamos a implementação em Python descrita a seguir, que é um solver simples de sistemas lineares, em que dados \({A}\), \({b}\) e uma determinada precisão, retorna a solução \({x}\).

É importante salientar que, para simular precisão simples, alimentamos todas as operações algébricas do numpy com dados em precisão simples ("float32"). Dessa forma, o resultado de cada operação também será em precisão simples. Na codificação, em cada operação, ainda utilizamos o parâmetro dtype=precisao, mas não é necessário para nosso objetivo, visto que o numpy realizará a operação na precisão desejada caso os operandos estejam na mesma precisão.

Observação 1: Porém, se em um operação algébrica, determinada entrada estiver um precisão dupla, enquanto as demais em precisão simples, o numpy "promove" a operação para a maior precisão.

Observação 2: Os resultados apresentados foram obtidos ao executar os scripts no Python 3.8.3, numpy 1.18.5, sobre o Windows 10 em uma arquitetura x64.

Basicamente, o algoritmo desenvolvido:

  1. Armazena \(A\) e \(b\) com precisão simples ou dupla, conforme o usuário define a variável precisão, para float32 ou float64, respectivamente. Para o exercício em questão, utilizaremos precisão simples conforme solicitado.

  2. Concatena \(A\) e \(b\) na matriz \(G\) e realiza o processo de Eliminação Gaussiana, com pivoteamento parcial (em que apenas a coluna atual é avaliada na busca do maior pivô. Esse feature foi implementado de forma opcional, e será desabilitado em alguns ensaios mais a frente). \({G}\) é retornada então no formato triangular superior (com relação a parte quadrada referente a \({A}\))

  3. Resolve o sistema linear a partir da última linha de \({G}\), em que é possível determinar diretamente o valor de \(x_k\), \(k\) índice da última linha, e a partir deste solucionar a linha imediatamente anterior, e assim consecutivamente.

Que é implementado por:

import numpy as np
import sys


def eliminacao_gaussiana(A, b, precisao, pivoteamento_parcial):
    ''' Dados A e b, realiza processo de Eliminação Gaussiana, 
    com ou sem pivoteamento parcial
    '''
    G = np.concatenate((A, b), axis=1)
    n = G.shape[0]

    # inicializa array de multiplicadores para as operações lineares da
    # Eliminação Gaussiana:
    m = np.ones(n, dtype=precisao)

    for k in range(n):

        if(pivoteamento_parcial == "on"):

            # obtem índice da linha que contém elemento da maior valor para a
            # coluna em questão.
            # note que consideramos só as linhas abaixo das já processadas:
            indice_linha_pivo = abs(G[k:, k]).argmax() + k

            if(indice_linha_pivo != k):

                linha_pivo = G[indice_linha_pivo, :].copy()
                G[indice_linha_pivo] = G[k]
                G[k] = linha_pivo

        if(G[k, k] != 0):

            for i in range(k + 1, n):

                m[i] = np.divide(G[i, k], G[k, k], dtype=precisao)

            for j in range(k, n + 1):

                for i in range(k + 1, n):

                    G[i, j] = np.subtract(
                        G[i, j], np.multiply(m[i], G[k, j], dtype=precisao),
                        dtype=precisao)

        else:

            sys.exit("Matriz A é singular")

    return G


def resolva(A, b, precisao, pivoteamento_parcial="on"):
    ''' Dados A e b, resolve o Sistema Linear, com pivoteamento parcial
    por default, ou sem caso definido na chamada
    '''
    G = eliminacao_gaussiana(A, b, precisao, pivoteamento_parcial)

    n = G.shape[0]
    x = np.zeros((n, 1), dtype=precisao)

    for i in range(n-1, -1, -1):

        x[i, 0] = np.divide(G[i, n], G[i, i], dtype=precisao)

        for j in range(i):

            G[j, n] = np.subtract(G[j, n], np.multiply(G[j, i], x[i, 0],
                                                       dtype=precisao),
                                  dtype=precisao)

    return x


def imprima(desc, array):
    ''' imprime o resultado, com todos os dígitos armazenados
    conforme precisão utilizada
    '''    
    print(desc, ":")
    [print("{}{} = {:0.99g}".format(desc, i[0]+1, p))
         for i, p in np.ndenumerate(array)]


#%% Main

if __name__ == '__main__':

    #%% a)

    # define precisão simples:
    precisao = "float32"

    A = np.array(
        [[21.0, 67.0, 88.0, 73.0],
         [76.0, 63.0,  7.0, 20.0],
         [ 0.0, 85.0, 56.0, 54.0],
         [19.3, 43.0, 30.2, 29.4]], dtype=precisao)

    b = np.array([[141.0],
                  [109.0],
                  [218.0],
                  [ 93.7]], dtype=precisao)

    x = resolva(A, b, precisao)

    print("\nitem a)\n")
    imprima("x", x)

Utilizando essa rotina com precisão simples, note que obtemos o seguinte resultado:

x1 = -1.00000035762786865234375
x2 = 2.000000476837158203125
x3 = -3.000000476837158203125
x4 = 4

E assim podemos notar a imprecisão com relação a resposta esperada (\(x=[-1, 2, -3, 4]\)).

Item b)

O residual \({r=b-Ax}\) é então calculado por meio de aritimética de precisão dupla, com o resultado sendo armazenado inicialmente na variável \({r_double}\), e então convertido em um vetor de precisão simples, \({r}\), conforme solicitado pela questão, pelo segmento de código abaixo:

    #%% b)

    # calcula o resíduo r = b - Ax, com precisão dupla:
    r_double = np.subtract(b, np.matmul(A, x, dtype="float64"), \
                           dtype="float64")

    # converte o resultado, r_double, para precisão simples, e armazena em r:
    r = r_double.astype(precisao)

    # imprime os desvios encontrados e a norma 1 de r:
    print("\nitem b)\n")
    imprima("r", r)
    print("||r||_1 =", np.linalg.norm((r), ord=1))

Os valores obtidos são:

r1 = 1.752376556396484375e-05
r2 = 4.76837158203125e-07
r3 = -1.3828277587890625e-05
r4 = 7.9870233093970455229282379150390625e-07

Com norma 1 (para efeitos de comparação posterior) igual a:

||r||_1 = 3.2627584e-05

Item c)

Resolve-se então o sistema linear \({Az=r}\), obtendo-se a solução "aprimorada", \({x+z}\), conforme segmento de código a seguir:

    #%% c)

    # encontra z tal que Az = r, realizando o cálculo com precisão simples:
    z = resolva(A, r, precisao)

    # calcula a solução "aprimorada", x + z:
    x_apr = np.add(x, z, dtype=precisao)

    print("\nitem c)\n")
    imprima("z", z)
    print("\nApós aprimoramento:")
    imprima("x", x_apr)

Em que \({z}\) obtido é:

z1 = 3.5762786865234375e-07
z2 = -4.76837158203125e-07
z3 = 4.76837158203125e-07
z4 = 0

E a solução aprimorada é:

x1 = -1
x2 = 2
x3 = -3
x4 = 4

Assim, com apenas um "aprimoramento", chega-se na solução exata por meio do procedimento descrito até então.

Para que seja possível explorar o item d da questão, ou seja, obter-se mais passos de refinamento, tentamos então reduzir a precisão do processo de Eliminação Gaussiana, "desligando" o passo de pivoteamento parcial, e assim mostrando que dessa forma realmente pioramos a precisão da resposta. Para isso, executa-se o seguinte segmento de código:

    #%% c) alternativa

    x = resolva(A, b, precisao, pivoteamento_parcial="off")

    # realiza a operação r = b - Ax:
    r = np.subtract(b, np.matmul(A, x, dtype="float64"), dtype="float64") \
        .astype(precisao)

    # encontra z tal que Az = r, realizando o cálculo com precisão simples:
    z = resolva(A, r, precisao, pivoteamento_parcial="off")

    # calcula a solução "aprimorada", x + z:
    x_apr = np.add(x, z, dtype=precisao)

    print("\nitem c) sem pivoteamento parcial:\n")
    imprima("x", x)
    imprima("r", r)
    print("||r||_1 =", np.linalg.norm((r), ord=1))
    imprima("z", z)
    print("\nApós aprimoramento:")
    imprima("x", x_apr)

Como resultado "pré aprimoramento", obtemos:

x1 = -0.99985682964324951171875
x2 = 2.0006465911865234375
x3 = -2.9977271556854248046875
x4 = 3.996625423431396484375

r1 = 5.60283660888671875e-06
r2 = -3.45706939697265625e-05
r3 = -1.239776611328125e-05
r4 = 6.04099886913900263607501983642578125e-06

Podemos notar que a norma 1 do resíduo é quase o dobro da solução que utiliza pivotamento parcial:

||r||_1 = 5.8612295e-05

E após resolvido \({Az=r}\), a solução "aprimorada" \({x+z}\) é:

x1 = -0.999999940395355224609375
x2 = 2.0000002384185791015625
x3 = -2.99999904632568359375
x4 = 3.9999988079071044921875

Ou seja, ainda há oportunidade de melhoria.

0 votos
respondida Out 26 por Carlos Piccioni (21 pontos)  
editado Nov 11 por Carlos Piccioni

Item d)

Por fim, implementamos a função abaixo que repete os passos realizados nos item c) e d), buscando aprimorar a solução até que a norma 1 de \(r\) seja zero ou pare de diminuir:

def resolva_apr(A, b, precisao, pivoteamento_parcial="on"):

    norma_r_ant = 10 ** 5
    count = 0

    x = resolva(A, b, precisao, pivoteamento_parcial="off")

    r = np.subtract(b, np.matmul(A, x, dtype="float64"), dtype="float64") \
        .astype(precisao)

    norma_r = np.linalg.norm((r), ord=1)

    print("Solução original:")
    imprima("x", x)
    imprima("r", r)
    print("||r||_1 =", norma_r)

    while(not ((norma_r_ant <= norma_r) or (norma_r == 0))):

        norma_r_ant = norma_r
        count = count + 1

        z = resolva(A, r, precisao, pivoteamento_parcial="off")

        x = np.add(x, z, dtype=precisao)

        r = np.subtract(b, np.matmul(A, x, dtype="float64"), dtype="float64") \
            .astype(precisao)

        norma_r = np.linalg.norm((r), ord=1)

        print("\nRefinamento", count, ":")
        imprima("z", z)
        imprima("x", x)
        imprima("r", r)
        print("||r||_1 =", norma_r)

Então a chamamos, sem utilizar pivoteamento parcial:

    x = resolva_apr(A, b, precisao, pivoteamento_parcial="off")

E como resultado obtemos:

item d)

Solução original:
x :
x1 = -0.99985682964324951171875
x2 = 2.0006465911865234375
x3 = -2.9977271556854248046875
x4 = 3.996625423431396484375
||r||_1 = 5.8612295e-05

Refinamento 1 :
z :
z1 = -0.00014311590348370373249053955078125
z2 = -0.0006463497993536293506622314453125
z3 = -0.00227199564687907695770263671875
z4 = 0.00337331625632941722869873046875
x :
x1 = -0.999999940395355224609375
x2 = 2.0000002384185791015625
x3 = -2.99999904632568359375
x4 = 3.9999988079071044921875
r :
r1 = -1.4126300811767578125e-05
r2 = -2.384185791015625e-06
r3 = -9.2983245849609375e-06
r4 = -5.1558026825659908354282379150390625e-06
||r||_1 = 3.0964613e-05

Refinamento 2 :
z :
z1 = -5.958870730182752595283091068267822265625e-08
z2 = -2.38346586911575286649167537689208984375e-07
z3 = -9.53421476879157125949859619140625e-07
z4 = 1.1917173878828180022537708282470703125e-06
x :
x1 = -1
x2 = 2
x3 = -3
x4 = 4
r :
r1 = 0
r2 = 0
r3 = 0
r4 = 0
||r||_1 = 0.0

Ou seja, são necessárias duas etapas de aprimoramento para que não haja mais redução do resíduo. Em nosso exemplo, com precisão simples, chegamos ao resultado exato, com resíduo zero.

Complemento:

Executando os passos a) e b) com precisão dupla e pivoteamento parcial, ou seja, o calculando o sistema linear e seu resíduo (para tal, basta alterar a variável precisao para float64 no código apresentado no item c), encontramos um \({r}\) com norma 1 na ordem \(1.1\times 10^{-13}\), o que indica precisão consideravelmente superior se comparado à norma 1 do resíduo na resolução por precisão simples, na ordem de \(3.3\times 10^{-05}\).

Referências:

www.lia.ufc.br/~valdisio/download/ieee.pdf
www.geeksforgeeks.org/floating-point-representation-basics
www.inf.pucrs.br/emoreno/undergraduate/EC/arqi/class_files/Aula09.pdf
docs.python.org/3/tutorial/floatingpoint.html
pt.wikipedia.org/wiki/V%C3%ADrgula_flutuante

...