Aula 4 - SciPy

Observação: As notas a seguir estão em constante desenvolvimento, aperfeiçoamento e revisões.

Introdução


Essa aula é baseada principalmente no material disponibilizado gratuitamente em http://scipy-lectures.org/intro/scipy.html sob a licença CC BY 4.0.

  • O que é?

O pacote SciPy contém um conjunto de toolboxes para solução dos problemas mais frequentes da Computação Científica. Cada um de seus sub-módulos é dedicado a um problema específico. A utilização da SciPy foi feita para ser aplicada eficientemente com a NumPy.

  • Por que eu deveria usar?

Antes de se implementar algo objetivando solucionar um problema, é interessante verificar se já não existe uma solução disponibilizada, ou mesmo parte da solução do problema, objetivando ganhar tempo ou se focar apenas no que se deseja adicionar (a sua contribuição original). No meio acadêmico, composto em vasta maioria por programadores não-profissionais, existe uma tendência a reinventar a roda, e isso leva à criação de códigos que apresentam bugs, não-otimizados, difíceis de compartilhar para outros usuários, com desenvolvimento sem elaboração de testes apropriados e com baixa manutenibilidade, todos esses atributos indesejados em termos de Arquitetura e Projeto de Software. Por outro lado, a SciPy oferece rotinas otimizadas e testadas, o que torna seu uso como preferencial sempre que possível.

Álgebra Linear computacional (solução de sistemas lineares)

Esse módulo da scipy oferece recursos para solucionamento de sistemas lineares e operações em matrizes. Ele basicamente herda todas as funcionalidades disponibilizadas pela numpy, com ligeiras mudanças, e algumas rotinas adicionais.

In [1]:
import scipy.linalg as linalg

# Alguns pré-requisitos
import numpy as np
import numpy.random as rnd

Solução de sistemas lineares (abordagem densa)

A rotina é bastante similar ao equivalente em numpy, inclusive é um wrapper da LAPACK.

In [2]:
n = 100
shape = (n, n)
A = rnd.randint(2, size=shape)
b = rnd.randint(5, size=(n,))

A, b
Out[2]:
(array([[1, 0, 1, ..., 1, 1, 0],
        [1, 0, 1, ..., 0, 0, 1],
        [0, 1, 0, ..., 1, 1, 0],
        ...,
        [0, 0, 0, ..., 0, 1, 0],
        [1, 0, 0, ..., 1, 0, 0],
        [0, 0, 1, ..., 1, 1, 0]]),
 array([2, 3, 3, 2, 1, 3, 3, 3, 1, 4, 2, 3, 4, 4, 0, 0, 0, 2, 0, 2, 1, 2,
        4, 4, 0, 3, 0, 2, 2, 0, 3, 4, 1, 2, 0, 2, 4, 4, 4, 3, 1, 0, 2, 0,
        4, 4, 0, 1, 4, 4, 1, 4, 0, 4, 1, 3, 0, 2, 0, 2, 3, 2, 3, 4, 4, 1,
        4, 4, 4, 1, 2, 2, 3, 1, 2, 2, 2, 4, 4, 4, 2, 1, 4, 1, 3, 2, 1, 0,
        4, 0, 3, 4, 1, 1, 4, 3, 1, 1, 2, 3]))
In [3]:
linalg.det(A)
Out[3]:
-1.9738270951972973e+47
In [4]:
import matplotlib.pyplot as plt

plt.spy(A)
plt.show()
<Figure size 640x480 with 1 Axes>

Solucionando:

In [5]:
x = linalg.solve(A, b)

x
Out[5]:
array([ -82.7164806 ,   -2.96270617, -151.06487509,  -19.08967755,
          8.83862836,   11.34653803,   -5.44776064,   26.92565362,
        158.83444399,  -51.12905166, -134.08341032,  -17.66280145,
          0.80419445,  -11.68831808,  -29.90892537,  -12.23240838,
         29.86457556,  -46.29789116,   72.93453583,  -42.22197583,
         13.79321434,   56.09927433,   59.40283784,   16.94974588,
        -65.30227781,   -3.25734712,   62.61785503,   51.00830103,
        102.76151776,  -17.12401128,  150.30331322,   83.10413748,
        -15.60689603,  -15.07398606, -137.37444284,   24.86303914,
         48.89025745,  -41.41852994,  -47.71515084,   73.19304835,
        -37.72700976,   26.9577809 ,  -83.98343844, -189.55226606,
         35.29619078,   15.41808593,  -41.59848702,    0.78453521,
        -64.84137661,  -31.87531852,   13.0402689 ,  -24.99819604,
        101.63022238,   45.9238573 , -109.83604935,   61.88247869,
         72.5279913 ,   27.88962286,  -38.20705127,  -19.0390928 ,
       -108.84023019,   24.44441981,    5.58167811,  -55.41652727,
         12.70834201,  -47.58408435,   -5.67149441,   30.72536621,
        -63.37401229,   30.03328435,   71.65085963,   73.60675895,
         -5.42522681,   32.80862951,   -7.87702417,   85.98257389,
         23.84319741,   16.62380602,  -11.1597425 ,  -91.69084818,
          7.66973906,   33.8967038 ,  -42.05255591,   31.82730897,
         14.11558382,  -67.7580847 ,  -25.11368995,   48.8321148 ,
        -49.46817067,   13.35482847,    0.29641158,   57.19453605,
        -70.51908646,  -24.5865773 ,  -18.50695833,   25.92815729,
         43.63689662,   47.7744739 ,   81.81421388,   12.51963215])
In [6]:
np.allclose(np.dot(A, x), b)
Out[6]:
True

Solucionando com Fatoração LU

A fatoração é feita da seguinte forma: $A \equiv P L U$. Para resolver um sistema linear com fatoração $LU$, os seguintes passos são aplicados:

  1. Resolver o sistema $L y = P b$ para obter $y$;
  2. Por fim, obtemos $x$ resolvendo $U x = y$.

Vamos obter os fatores $P L U$ (permutador de linhas, diag. inferior, diag. superior):

In [7]:
P, L, U = linalg.lu(A)

P, L, U
Out[7]:
(array([[1., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 1., 0., ..., 0., 0., 0.],
        ...,
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.],
        [0., 0., 0., ..., 0., 0., 0.]]),
 array([[ 1.        ,  0.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 0.        ,  1.        ,  0.        , ...,  0.        ,
          0.        ,  0.        ],
        [ 1.        ,  1.        ,  1.        , ...,  0.        ,
          0.        ,  0.        ],
        ...,
        [ 1.        ,  1.        , -0.        , ...,  1.        ,
          0.        ,  0.        ],
        [ 0.        ,  1.        , -1.        , ..., -0.32079621,
          1.        ,  0.        ],
        [ 1.        ,  0.        , -0.        , ...,  0.9338438 ,
          0.10281833,  1.        ]]),
 array([[ 1.        ,  0.        ,  1.        , ...,  1.        ,
          1.        ,  0.        ],
        [ 0.        ,  1.        ,  0.        , ...,  1.        ,
          1.        ,  0.        ],
        [ 0.        ,  0.        , -1.        , ..., -1.        ,
         -1.        ,  0.        ],
        ...,
        [ 0.        ,  0.        ,  0.        , ..., -3.96413061,
          2.17530996,  0.89069924],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         -0.22062118,  1.37243093],
        [ 0.        ,  0.        ,  0.        , ...,  0.        ,
          0.        , -0.65385452]]))
In [8]:
np.allclose(A - P @ L @ U, np.zeros(shape))
Out[8]:
True
In [9]:
plt.spy(L)
Out[9]:
<matplotlib.image.AxesImage at 0x7f6b7629b7b8>
In [10]:
plt.spy(U)
Out[10]:
<matplotlib.image.AxesImage at 0x7f6b76210908>
In [11]:
plt.spy(P)
Out[11]:
<matplotlib.image.AxesImage at 0x7f6b76171860>
  • Passo 1:
In [12]:
y = np.dot(linalg.inv(L) @ linalg.inv(P), b)

y
Out[12]:
array([ 2.00000000e+00,  3.00000000e+00, -3.00000000e+00,  1.00000000e+00,
        0.00000000e+00,  3.00000000e+00,  5.00000000e+00, -7.00000000e+00,
        1.25000000e+00,  9.09090909e-02,  7.46031746e-01,  1.91812865e+00,
       -5.09248555e+00, -3.34549878e+00, -1.18552036e+00, -1.64918970e+00,
       -4.37334485e+00,  5.68544385e+00,  2.23095651e+00, -6.88061992e-01,
       -6.63810108e-01,  1.37091003e+00, -3.29277133e+00,  1.46491231e+00,
       -4.15189629e-01,  8.27074413e-02,  5.12804057e+00, -5.02379511e-01,
       -3.38545905e+00,  1.99528521e+00, -5.21577765e+00, -9.96031400e+00,
       -5.02156203e+00,  8.73090859e+00, -1.20800766e+01, -8.33179578e+00,
        1.60915685e+00,  2.96917710e+00, -8.29789524e-01,  5.79510066e+00,
       -3.43599097e+00,  1.56036856e+00,  4.73979046e+00,  2.47906052e+00,
       -6.43738781e+00,  3.09555037e+00, -2.33941211e+00, -1.23910986e+00,
        7.41871640e-01, -8.25333255e+00,  1.00820310e+01,  4.39949192e+00,
       -2.00962641e+00,  1.18822812e-01,  6.05320961e+00, -2.53885729e+00,
       -6.01627560e-01,  5.12565866e-03, -3.99906934e+00, -2.69806990e+00,
       -2.08615118e+00,  2.05599280e+00,  3.03696550e+00,  1.19372849e+00,
       -3.16060549e+00, -2.08547544e+00, -2.28442021e+00, -3.16178591e+00,
       -4.30598592e+00,  3.27770729e+00, -8.80404019e-01, -9.24769571e-01,
       -6.72847811e+00, -6.69723268e+00, -5.40733747e+00,  6.91752292e+00,
        7.11884419e+00, -1.18182432e+00,  6.98971354e+00, -9.46066366e-01,
       -5.00368343e+00, -5.02265003e+00,  6.85190500e+00,  7.98481987e+00,
        6.16031424e+00, -1.39713794e+01, -7.10369209e+00, -6.68454489e+00,
       -9.12625079e+00,  4.41008589e+00, -2.53760305e-01, -6.24515354e+00,
       -2.99698779e+00, -2.40608179e+00, -3.38808780e+00, -3.04715189e+00,
       -3.98671724e-01, -2.61753197e-01, -8.67617993e-01, -8.18601803e+00])
  • Passo 2:
In [13]:
x = np.dot(linalg.inv(U), y)

x
Out[13]:
array([ -82.7164806 ,   -2.96270617, -151.06487509,  -19.08967755,
          8.83862836,   11.34653803,   -5.44776064,   26.92565362,
        158.83444399,  -51.12905166, -134.08341032,  -17.66280145,
          0.80419445,  -11.68831808,  -29.90892537,  -12.23240838,
         29.86457556,  -46.29789116,   72.93453583,  -42.22197583,
         13.79321434,   56.09927433,   59.40283784,   16.94974588,
        -65.30227781,   -3.25734712,   62.61785503,   51.00830103,
        102.76151776,  -17.12401128,  150.30331322,   83.10413748,
        -15.60689603,  -15.07398606, -137.37444284,   24.86303914,
         48.89025745,  -41.41852994,  -47.71515084,   73.19304835,
        -37.72700976,   26.9577809 ,  -83.98343844, -189.55226606,
         35.29619078,   15.41808593,  -41.59848702,    0.78453521,
        -64.84137661,  -31.87531852,   13.0402689 ,  -24.99819604,
        101.63022238,   45.9238573 , -109.83604935,   61.88247869,
         72.5279913 ,   27.88962286,  -38.20705127,  -19.0390928 ,
       -108.84023019,   24.44441981,    5.58167811,  -55.41652727,
         12.70834201,  -47.58408435,   -5.67149441,   30.72536621,
        -63.37401229,   30.03328435,   71.65085963,   73.60675895,
         -5.42522681,   32.80862951,   -7.87702417,   85.98257389,
         23.84319741,   16.62380602,  -11.1597425 ,  -91.69084818,
          7.66973906,   33.8967038 ,  -42.05255591,   31.82730897,
         14.11558382,  -67.7580847 ,  -25.11368995,   48.8321148 ,
        -49.46817067,   13.35482847,    0.29641158,   57.19453605,
        -70.51908646,  -24.5865773 ,  -18.50695833,   25.92815729,
         43.63689662,   47.7744739 ,   81.81421388,   12.51963215])
In [14]:
np.allclose(np.dot(A, x), b)
Out[14]:
True

Observação: Não precisamos fazer isso na mão, a linalg.solve() nada mais é do que uma fatoração LU por baixo dos panos, toda otimizada e testada, wrapper do LAPACK.

Solução de sistemas esparsos

Em aplicações científicas, frequentemente nos deparamos com sistemas esparsos (com grande número de entradas nulas). Um exemplo típico de uma matriz advinda da discretização por elementos finitos seria:

In [15]:
import scipy.io as io

A_lshp = io.mmread('data/lshp3466.mtx')
plt.spy(A_lshp, markersize=1)
plt.show()

Esta matriz foi retirada do Matrix Market, um conhecido site de matrizes para benchmarks de solvers.

Podemos notar que apenas algumas entradas (marcadas em azul) são não-nulas. As outras entradas ocupam igual espaço na memória, o que é bastante desnecessário, já que são entradas nulas. Inclusive estas podem adicionar erros numéricos ao resultado. O scipy já realiza a leitura da matriz do formato Matrix Market (.mtx) como uma estrutura adequada à representação de matrizes esparsas.

In [16]:
import scipy.sparse as sp

sp.isspmatrix(A_lshp)
Out[16]:
True

A quantidade de entradas não-nulas é:

In [17]:
nnz = A_lshp.getnnz()

nnz
Out[17]:
23896

E o número total de entradas (incluindo as nulas) seria:

In [18]:
n = A_lshp.shape[0]
n_entries = n * n

n_entries
Out[18]:
12013156

Ou seja, a porcentagem de entradas não-nulas é:

In [19]:
nnz / n_entries * 100
Out[19]:
0.19891525590777312

Menos de 1%!!!

Agora, para solucionar esse sistema, vamos criar um vetor $b$ unitário apenas para testar.

In [20]:
b = np.ones(n)

E, para comparação, vamos pegar a matriz cheia:

In [21]:
A_full = A_lshp.toarray()

E vamos tranformar de COO (formato padrão de leitura do mmread) para o CRS (Compressed Row Storage), que é necessário para o solver da SciPy:

In [22]:
A_lshp = A_lshp.tocsr()
  • Solucionando o sistema denso:
In [23]:
x_dense = linalg.solve(A_full, b)
In [24]:
x_dense
Out[24]:
array([ 0.19755386, -0.0824912 ,  0.14322   , ...,  0.60159886,
        0.36378233, -0.42355086])

Agora, para avaliar a performance em termos de custo computacional temporal:

In [25]:
%%timeit
x_dense = linalg.solve(A_full, b)
880 ms ± 41 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
  • Solucionando o sistema esparso:

Vamos utilizar o solver esparso linalg.spsolve() equivalente ao solver denso:

In [26]:
import scipy.sparse.linalg as slinalg

Solucionando:

In [27]:
x_sparse = slinalg.spsolve(A_lshp, b)
In [28]:
x_sparse
Out[28]:
array([ 0.19755386, -0.0824912 ,  0.14322   , ...,  0.60159886,
        0.36378233, -0.42355086])
In [29]:
np.allclose(x_dense, x_sparse)
Out[29]:
True

E avaliando a performance em termos de tempo de execução:

In [30]:
%%timeit
slinalg.spsolve(A_lshp, b)
16.8 ms ± 915 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Disso, vemos que a diferença em tempo é quase 100 vezes menor! E esse não é o maior dos sistemas, obviamente.

Caso deseje obter maior profundidade do que é possível fazer em termos de solucionamento de sistemas esparsos e manipulações, a documentação oficial é excelente. Maiores informações sobre a classe de matrizes esparsas da scipy pode ser encontrada lá também. Aqui, mostramos apenas o básico para efeitos de conhecimento.

Encontrando raízes de equações com scipy

As rotinas de solucionar equações fazem parte do módulo de otimização da SciPy scipy.optimize. Isso se deve principal pela relação entre a relação entre problemas de otimização não-linear com problemas de sistemas de equações não-lineares.

In [31]:
import scipy.optimize as opt

Vamos usar como problema base a equação de Colebrook-White para o fator de fricção de escoamento turbulento em tubos:

$$ \frac{1}{\sqrt{f}} = 1.14 - 2 \log_{10}\left(\frac{\varepsilon}{D} + \frac{9.35}{\text{Re}\sqrt{f}}\right) $$

onde $f$ é o fator de fricção, $\varepsilon$ é a rugosidade absoluta, $D$ é o diâmetro interno e Re é o número de Reynolds associado ao escoamento do fluido. Essa correlação - uma equação transcendental - é válida para regimes onde $\text{Re} > 4000$.

Por ser uma equação transcendental, é impossível isolar nossa quantidade de interesse, que é o fator de fricção, por meio de manipulações matemáticas. A única forma de solucionar esse problema é numericamente. Então, vamos formular esse problema como um problema de encontrar raízes:

In [32]:
def colebrook_white(f, eps, D, Re):
    return 1.14 - 2 * np.log10(eps / D + 9.35 / (Re * np.sqrt(f))) - 1 / np.sqrt(f)

Os parâmetros serão dados para um duto de ferro fundido:

In [33]:
eps = 0.26  # Em mm
D = 200.0  # Em mm

E o número de Reynolds:

In [34]:
Re = 10000.

Resolvendo com métodos livres de informação de derivadas

De início, vamos começar com métodos sem informações de derivadas, que tem garantia de convergência dado um intervalo de busca. Para exemplificar, veremos uma simples bisseção:

In [35]:
x_min, x_max = 1e-3, 1.0  # Intervalo de busca
xsol = opt.bisect(colebrook_white, x_min, x_max, args=(eps, D, Re,))

xsol
Out[35]:
0.032825352944912074

Caso interesse verificar informações do processo de solucionamento, basta fazer:

In [36]:
xsol, log = opt.bisect(colebrook_white, x_min, x_max, args=(eps, D, Re,), full_output=True)

log
Out[36]:
      converged: True
           flag: 'converged'
 function_calls: 41
     iterations: 39
           root: 0.032825352944912074

Dentre os métodos livres de derivada disponíveis pela scipy, o melhor é o TOMS 748, com taxa de convergência chegando até 2.7, superior ao método de Newton, que utiliza derivadas. Veremos o TOMS 748 em ação:

In [37]:
xsol, log = opt.toms748(colebrook_white, x_min, x_max, args=(eps, D, Re,), full_output=True)

log
Out[37]:
      converged: True
           flag: 'converged'
 function_calls: 18
     iterations: 7
           root: 0.03282535294547174

Bem melhor, hein? A principal vantagem desse método é a necessidade de definir apenas o intervalo de busca (sem necessidade de estimativa inicial), além de que a convergência será atingida, mesmo que seja em um dos limites.

Resolvendo com métodos que utilizam derivadas

Para resolver com o método de Newton, em vez de utilizar um intervalo de busca, fornecemos uma estimativa inicial para a solução.

In [38]:
x0 = 0.05
xsol, log = opt.newton(colebrook_white, x0, args=(eps, D, Re,), full_output=True)

log
Out[38]:
      converged: True
           flag: 'converged'
 function_calls: 8
     iterations: 7
           root: 0.03282535294547166

Em geral, a taxa de convergência de métodos que utilizam derivadas, em geral, é maior. Porém, a convergência é condicionada à qualidade do chute inicial, podendo não convergir.

Resolvendo sistemas não-lineares

Para sistemas não-lineares, a interface muda ligeiramente em relação ao caso de uma única variável independente. Neste caso, se usa a interface geral optimize.root(). O método padrão dessa chamada é uma modificação do método de Powell híbrido da MINPACK.

$$ \left\{ \begin{align} &x + 0.5 (x - y)^3 = 1\\ &0.5 (y - x)^3 + y = 0 \end{align} \right. $$

In [39]:
def func(x):
    f = (
        x[0]  + 0.5 * (x[0] - x[1])**3 - 1.0,  # Eq 1
        0.5 * (x[1] - x[0])**3 + x[1]  # Eq 2
    )  # retorna uma lista onde cada entrada é uma linha do sistema
    return f
In [40]:
x0 = [0.0, 0.0]
solution = opt.root(func, x0)

solution.x
Out[40]:
array([0.8411639, 0.1588361])

Problemas de Otimização Não-linear

Problemas irrestritos

Aproveitando o embalo, já que é no mesmo sub-módulo, vamos logo ver como solucionar problema de otimização (local e global). Como primeiro exemplo, vamos tomar a função de Rosenbrock como função-objetivo, que é um conhecido benchmark para problemas de otimização. Para efeitos práticos, vamos tomar a seguinte forma em 2D:

$$ f(x, y) := (a - x)^2 + b \left(y - x^2\right)^2 $$

O mínimo local (que é único e global neste caso) ocorre quando $f(x, y) = 0$, que é obtido com $(x, y) = (a, a^2)$. Aqui usaremos $a = 1$ e $b = 100$.

In [41]:
rosenbrock = lambda x: (1 - x[0])**2.0 + 100 * (x[1] - x[0]**2.0)**2.0

Minimização local

Para resolver o problema de otimização local, usamos a interface genérica optimize.minimize(). Com ela, podemos acessar todos os métodos baseados em derivadas e os livres de derivada, mas para otimização local.

In [42]:
x0 = [-10, 10]
result = opt.minimize(rosenbrock, x0)

result
Out[42]:
      fun: 2.6957239738060108e-11
 hess_inv: array([[0.49170495, 0.98335213],
       [0.98335213, 1.97158867]])
      jac: array([-4.63905237e-07, -4.81389906e-07])
  message: 'Optimization terminated successfully.'
     nfev: 504
      nit: 98
     njev: 126
   status: 0
  success: True
        x: array([0.99999481, 0.99998961])

Lembrando que o mínimo local, nesse caso, é global e ocorre em $(1, 1)$. Para acessar apenas o resultado, acesse o atributo:

In [43]:
result.x
Out[43]:
array([0.99999481, 0.99998961])

O método padrão é o BFGS, um método de quasi-Newton. Digamos que desejamos resolver com Gradiente Conjugado, então devemos passar no argumento opcional method:

In [44]:
x0 = [-10, 10]
result = opt.minimize(rosenbrock, x0, method='CG')

result
Out[44]:
     fun: 3.2706500031004414e-11
     jac: array([-7.33855221e-06,  2.42827904e-06])
 message: 'Optimization terminated successfully.'
    nfev: 292
     nit: 30
    njev: 73
  status: 0
 success: True
       x: array([0.99999428, 0.99998857])

Uma lista completa de métodos e argumentos possíveis está disponibilizado aqui.

Minimização global

A SciPy oferece algumas opções de métodos de Otimização Global. Um dos métodos que estão disponibilizados é bastante utilizado (e desenvolvido em pesquisas!) no LNCC, que é o Evolução Diferencial. Por isso, em homenagem aos colegas, vamos esse método. Trata-se de um método meta-heurístico estocástico que não é bioinspirado, porém bastante eficiente, paralelizável e (inesperadamente) simples. Vamos ver em ação!

Para utilizá-lo, devemos fornecer o intervalo de busca em formato de hipercubo.

In [45]:
bounds = [(0, 2), (0, 2)]  # Informando o intervalo de busca em cada componente
In [46]:
result = opt.differential_evolution(rosenbrock, bounds)

result
Out[46]:
     fun: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 3783
     nit: 125
 success: True
       x: array([1., 1.])

A estratégia padrão de evolução utilizada é a best1bin. Digamos que o usuário deseje utilizar a estratégia clássica e originalmente proposta no primeiro paper de Storn e Price, que seria a rand1bin. Além disso, desejamos utilizar uma população de 50 indivíduos, sendo que o padrão é 15. Para termos essa configuração, basta utilizar os parâmetros opcionais como abaixo:

In [47]:
result = opt.differential_evolution(rosenbrock, bounds, popsize=50, strategy='best1bin', mutation=0.2)

result
Out[47]:
     fun: 0.0
 message: 'Optimization terminated successfully.'
    nfev: 4003
     nit: 39
 success: True
       x: array([1., 1.])

Problemas com restrições

Vamos tomar como exemplo o problema:

Achar:

$$ \min_{x_1, x_2 \in \mathbb{R}}\, (x_1 + 1)^2 + (x_2 - 1)^2 $$

sujeito a

$$ \begin{align} 2 x_2 - 1 &= 0\\ (1- x_1)\left(4 - x_1^2 - x_2^2\right) &\leq 0\\ 2x_1^2 + x_2^2 - 100 &\leq 0 \end{align} $$

que possui o mínimo global marcado em x (sem restrições) na figura abaixo. As regiões preenchidas e ao longo da linha roxa são o conjunto viável.

ineqs2.png

Sob esse conjunto de restrições, a função possui dois mínimos locais em

$$ x_1^* = \left(1, \frac{1}{2} \right); \quad x_2^* = \left(-\frac{\sqrt{15}}{2}, \frac{1}{2}\right) $$

Porém, $f_{objetivo}(x_1^*) > f_{objetivo}(x_2^*)$.

A função-objetivo e sua Jacobiana é definida como:

In [48]:
def f_obj(x):
    return (x[0] + 1)**2 + (x[1] - 1)**2


def f_deriv(x):
    dfdx0 = 2 * x[0] + 2
    dfdx1 = 2 * x[1] - 2
    return np.array([dfdx0, dfdx1])

As restrições são passadas para a rotina com a seguinte forma:

In [49]:
constraints = (
    {'type': 'eq', 'fun': lambda x: np.array([2 * x[1] - 1])},
    {'type': 'ineq', 'fun': lambda x: np.array([-(1 - x[0]) * (4 - x[0]**2 - x[1]**2)])},
    {'type': 'ineq', 'fun': lambda x: np.array([-1 * (2 * x[0]**2 + x[1]**2 - 100)])}
)

Por fim, obtemos o mínimo como:

In [50]:
result = opt.minimize(f_obj, [-3.0, 1.0], jac=f_deriv,
                  constraints=constraints, method='SLSQP', options={'disp': True})
result.x
Optimization terminated successfully.    (Exit mode 0)
            Current function value: 1.127016654329262
            Iterations: 6
            Function evaluations: 6
            Gradient evaluations: 6
Out[50]:
array([-1.93649167,  0.5       ])

Onde foi utilizado o método Sequential Least SQuares Programming. Neste exemplo, podemos ver como passar restrições de igualdade e desigualdade. É importante notar que as desigualdades são fornecidas na forma $g(x_1, x_2) \ge 0$, ao contrário de como está na formulação do problema.

Interpolação

Interpolação é uma técnica de aproximação que pode ser tanto usada para obter funções mais "simples" (em algum sentido) a partir de funções mais complicadas ou obter uma função a partir de um conjunto de dados experimentais (que também pode ser feito com técnicas de ajuste. Aqui, vamos ver apenas um simples caso para conhecermos a ferramenta pela scipy.

Vamos simular dados experimentais que tem comportamento similar à função seno ao longo do tempo:

In [51]:
measured_points_number = 10
measured_time = np.linspace(0, 1, measured_points_number)
noise = np.random.random(measured_points_number) * 2 - 1
measures = np.sin(2 * np.pi * measured_time) + noise

Vamos agora importar a função de interpolação linear da scipy:

In [52]:
from scipy.interpolate import interp1d

E realizar a interpolação:

In [53]:
linear_interp = interp1d(measured_time, measures)
In [54]:
linear_interp(0.5)
Out[54]:
array(-0.30457695)

Agora temos o interpolador montado a partir dos dados medidos (measures). Com isso, podemos entrar com um novo conjunto de pontos (dentro do intervalo dos pontos, já que é uma interpolação) para ver os resultados interpolados desse novo conjunto:

In [55]:
interpolation_time = np.linspace(0, 1, 50)
linear_results = linear_interp(interpolation_time)

Podemos fazer uma rápida visualização dos dados experimentais virtuais e da interpolação da seguinte maneira:

In [56]:
import matplotlib.pyplot as plt

plt.figure(figsize=(5, 3))
plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(interpolation_time, linear_results, '-', label='Linear interp')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.grid(True)
plt.savefig('figura.png')
plt.show()

Lembrando que a interpolação é um procedimento "exato" nos pontos evaluados, o que pode gerar uma característica indesejada chamada overfitting, que prejudica fortemente a capacidade preditiva da interpolação.

Para complementar, vamos ver como seria uma interpolação cúbica. Para isso, usamos a rotina genérica de interpolação 1D interp1d.

In [57]:
cubic_interp = interp1d(measured_time, measures, kind='cubic')
cubic_results = cubic_interp(interpolation_time)

E agora, vamos visualizar novamente:

In [58]:
import matplotlib.pyplot as plt

plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(interpolation_time, linear_results, '-', label='Linear interp')
plt.plot(interpolation_time, cubic_results, '-', label='Cubic interp')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.show()

Caso seja do interesse de alguém da audiência, uma série de outras possibilidades estão disponíveis na SciPy para interpolação, inclusive com multivariáveis. Para maiores detalhes, verifique a documentação.

Ajuste de curvas

Outra forma de obtermos funções a partir de dados experimentais é por meio de ajuste de curvas. Pela SciPy, podemos fazer isso a partir de algumas maneiras. As mais notáveis são ajustes de parâmetros de uma função por minimização (global) de alguma métrica de erro, ou fazendo o ajuste de curvas já pronto pela SciPy.

Ajuste por minimização de quadrado do resíduo

Essa é uma forma praticamente black box de solucionar o problema com otimização global. Dizemos qual o modelo e os parâmetros que desejamos ajustar nesse modelo. Aqui, por conveniência, escolhemos uma função seno para efeitos de simplificação. O problema que desejamos resolver é o seguinte:

Achar os parâmetros $a, b, c \in \mathbb{R}$ tal que $$ \min \sum_{i=1}^{N_{exp}} (y_i - f(x, a, b, c))^2 $$

sendo $y_i$ os valores medidos experimentalmente e

$$ f(x, a, b) := a \cdot \sin(b x) $$

a função que desejamos ajustar. Note que essa função poderia ser qualquer coisa que dependa de parâmetros que desejamos encontrar para minimizar o quadrado do resíduo. Vamos definir aqui essas funções como:

In [59]:
def square_error(par, x, f_exp):
    a, b, c = par
    residual = f_exp - f(x, a, b, c)
    return np.sum(residual**2.0)


def f(x, a, b, c):
    return a * np.sin(b * x) + c

E resolver o problema de minimização global da forma:

In [60]:
bounds = [(-10, 10), (-10, 10), (-10, 10)]
result = opt.differential_evolution(
    square_error, 
    bounds=bounds, 
    args=(measured_time, measures,), 
    popsize=50,
    strategy='rand1bin',
    tol=1e-3,
    recombination=0.9,
    mutation=0.2
)

result
Out[60]:
     fun: 2.3595463311342963
     jac: array([ 1.15463195e-06, -3.10862447e-07, -1.19904087e-06])
 message: 'Optimization terminated successfully.'
    nfev: 5728
     nit: 37
 success: True
       x: array([-1.09337942, -6.76241452, -0.31493112])

Pegando os parâmetros ajustados:

In [61]:
a, b, c = result.x

Por fim, visualizando o resultado dos ajustes em relação aos dados experimentais:

In [62]:
import matplotlib.pyplot as plt

f_fitted_DE = f(measured_time, a, b, c)
plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(measured_time, f_fitted_DE, '-', label='DE fitting')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.show()

Note que, neste tipo de abordagem, ao contrário da interpolação, temos um problema de overfitting bem menor.

Ajuste de curva com a interface do scipy

Agora vamos realizar o ajuste de uma dada curva em dados experimentais a partir da interface da scipy que é a optimize.curve_fit. Nela, também se ajusta utilizando mínimos quadrados. No entanto, o método que utilizaremos aqui no coração da rotina é o Levenberg-Marquardt.

In [63]:
p0 = (3., 3., 2.)
params, params_covariance = opt.curve_fit(f, measured_time, measures, p0=p0, method='lm')

params, params_covariance
Out[63]:
(array([-0.95009848,  2.24907996,  0.37840766]),
 array([[ 0.94577318, -0.22953039, -0.64684902],
        [-0.22953039,  1.46367239,  0.14041387],
        [-0.64684902,  0.14041387,  0.5375403 ]]))
In [64]:
a, b, c = params
In [65]:
import matplotlib.pyplot as plt

f_fitted_LM = f(measured_time, a, b, c)
plt.plot(measured_time, measures, 'o', label='Experimental')
plt.plot(measured_time, f_fitted_LM, '-', label='DE fitting')
plt.plot(measured_time, f_fitted_DE, '-', label='LM fitting')
plt.legend()
plt.xlabel('time')
plt.ylabel('y')
plt.show()

O que mostra que um solver local pode não ser adequado para todos os problemas de ajuste. Portanto, muita atenção! Com frequência, o adequado e mais robusto pode ser um otimizador global. Que tal testar a estimativa inicial para o Levenberger-Marquardt com $(3, 3)$. O que foi obtido?

Quadraturas Numéricas

Em computação científica, frequentemente temos que calcular integrais. Por meio da SciPy, podemos fazer isso de maneira numérica muito facilmente.

In [66]:
from scipy import integrate

Integral simples

Vamos utilizar aqui a rotina integrate.quad(), que vem da implementação do QUADPACK em Fortran.

Um simples exemplo:

$$ F = \int_0^{\infty} e^{-x} dx $$

In [67]:
f = lambda x: np.exp(-x)
F, err = integrate.quad(f, 0, np.inf)

F, err
Out[67]:
(1.0000000000000002, 5.842606742906004e-11)

Com erro absoluto

In [68]:
err
Out[68]:
5.842606742906004e-11

Ou ainda, podemos fazer uma integração de uma função que tenha vários parâmetros. Por exemplo,

$$ F = \int_0^1 (a x + b) dx $$

In [69]:
f = lambda x, a, b: a * x + b
F, err = integrate.quad(f, 0, 1, args=(1, 1,))

F, err
Out[69]:
(1.5, 1.6653345369377348e-14)

Integrais duplas

A função integrate.dblquad() trata da integração dupla. Agora, temos que fornecer, em sequência nos argumentos, os limites de integração em x e depois os limites em y, porém os de y são em forma de generalizada de funções $gf(x)$ até $hf(x)$.

Como exemplo simples, vamos integrar:

$$ F = \int_{gf(x) = 0}^{hf(x) = 1}\int_0^1 x y^2 dx dy $$

In [70]:
f = lambda y, x: x * y**2
F, err = integrate.dblquad(
    f, 
    0,  # x_0
    1,  # x_1
    lambda x: 0,  # gf(x)
    lambda x: 1   # hf(x)
)

F, err
Out[70]:
(0.16666666666666669, 3.6927075527489535e-15)

Integrais triplas

A lógica da integral dupla se mantém e agora temos os limites em z que dependerá das variáveis integradas anteriormente. Portanto, os limites em z são definidos genericamente por funções $qf(x, y)$ e $rf(x, y)$.

Como exemplo, tomemos o caso:

$$ F = \int_{qf(x, y) = 0}^{rf(x, y) = 1}\int_{gf(x) = 2}^{hf(x) = 3}\int_1^2 x y z dx dy dz $$

In [71]:
f = lambda z, y, x: x * y * z
F, err = integrate.tplquad(f, 1, 2, lambda x: 2, lambda x: 3,
                           lambda x, y: 0, lambda x, y: 1)

F, err
Out[71]:
(1.8750000000000002, 3.324644794257407e-14)

Equações Diferenciais Ordinárias

No mesmo módulo que estão as quadraturas, a scipy também disponibiliza uma série de métodos para resolução de EDOs tanto para problemas rígidos como não-rígidos. A interface é bem simples. Considere o problema genérico de valor inicial:

$$ \dot{y} = f(t, y) $$

que pode ser inclusive um sistema de EDOs. Para a scipy, o usuário basicamente tem que fornecer o lado direito. Ou seja, deve fornecer $f(t, y)$.

Inicialmente, vamos considerar um dos problemas mais simples possível:

$$ \dot{y} = -2 y $$ com $y(t = 0) = 1$. Este é um problema dito autônomo, pois o lado direito não contém dependência temporal. Ou seja, $f(t, y) \equiv f(y)$.

Portanto, vamos definir o lado direito inicialmente:

In [72]:
def rhs(ypos, time):
    return -2 * ypos

Note que mesmo para problemas autônomos, devemos colocar o tempo como argumento.

Aqui, vamos usar a integrate.odeint():

In [73]:
from scipy.integrate import odeint

Por fim, integramos a EDO no tempo de 0 até 4, com 40 passos de tempo:

In [74]:
time = np.linspace(0, 4, 40)
y0 = 1
y = odeint(rhs, y0=y0, t=time)

Agora uma rápida visualização:

In [75]:
import matplotlib.pyplot as plt

plt.plot(time, y)
plt.xlabel('time')
plt.ylabel('y')
plt.show()

Um problema mais interessante: Sistema de Lorenz

Vamos considerar agora o problema do atrator de Lorenz, um sistema caótico que nasceu da modelagem de fluidos, mais especificamente em convecção atmosférica em torno da década de 1960. O sistema contém 3 equação diferenciais ordinárias não-lineares acopladas e é escrito da forma:

$$ \left\{ \begin{align} &\dot{u} = \sigma (v - u)\\ &\dot{v} = \rho u - v - uw\\ &\dot{w} = u v - \beta w\\ \end{align} \right. $$

Como todo bom sistema dinâmico não-linear (e caótico), é um problema rígido (stiff) e requer muita atenção na hora de resolvê-lo numericamente. Se o método não for adequado e capaz de resolver essa classe de problemas, a solução simplesmente pode não condizer com a "realidade" por motivos numéricos. Ou seja, é um problema de verificação numérica. Podemos estar resolvendo o modelo certo, porém o método o corrompe.

Para essa classe de problema, em geral se escolhe métodos implícitos do tipo Runge-Kutta ou BDF (Backward Differentiation Formulas), dentre outros. Aqui, vamos utilizar a implementação de Runge-Kutta implícito de quinta ordem conhecida como RADAU-IIA e fornecida pela SciPy.

Aqui vamos utilizar $\sigma = 10$, $\beta = 8 / 3$ e $\rho = 28$.

  • Condições iniciais:
In [76]:
u0, v0, w0 = 0.5, 0, 0.1
  • O lado direito do sistema:
In [77]:
def lorenz(t, X):
    sigma, beta, rho = 10, 8 / 3, 28  # parâmetros do modelo
    u, v, w = X
    u_prime = -sigma * (u - v)
    v_prime = rho * u - v - u * w
    w_prime = -beta * w + u * v
    return u_prime, v_prime, w_prime
  • Parâmetros da integração no tempo:
In [78]:
t0 = 0  # tempo inicial
tmax = 100  # tempo máximo
t_points = 500  # número de pontos no tempo
  • Solucionando o problema com RAUDAU-IIA e com RK45 (Runge-Kutta explícito de ordem 5) para efeitos comparativos:
In [79]:
solution_RADAU = integrate.solve_ivp(lorenz, t_span=(t0, tmax), y0=(u0, v0, w0,), method='Radau')
solution_RK = integrate.solve_ivp(lorenz, t_span=(t0, tmax), y0=(u0, v0, w0,), method='RK45')

t_computed_RADAU, y_computed_RADAU = solution_RADAU.t, solution_RADAU.y
t_computed_RK, y_computed_RK = solution_RK.t, solution_RK.y
# Separando as soluções
u_RADAU, v_RADAU, w_RADAU = y_computed_RADAU
u_RK, v_RK, w_RK = y_computed_RK
  • Visualizando a famosa "borbuleta de Lorenz":

RADAU:

In [80]:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(12, 12))
ax = fig.gca(projection='3d')
ax.plot(u_RADAU, v_RADAU, w_RADAU)
# ax.set_axis_off()
plt.show()

RK45:

In [81]:
fig = plt.figure(figsize=(12, 12))
ax = fig.gca(projection='3d')
ax.plot(u_RK, v_RK, w_RK)
# ax.set_axis_off()
plt.show()
  • Analisando o comportamento separado das componentes da solução em função do tempo: