domingo, 25 de outubro de 2015

Programação paralela em Python - parte 2

No último texto sobre programação, demos os primeiros passos na utilização do módulo multiprocessing, que permite a execução paralela de código em Python. Neste post iremos estudar mais a fundo as funcionalidades deste módulo usando como exemplo um pequeno programa de detecção de bordas em imagens.

Imagens digitais são normalmente representadas utilizando matrizes, onde a posição $(i, j)$ contém a "cor" presente no pixel correspondente da imagem. Neste texto trataremos somente de imagens em níveis de cinza. Portanto, nossas imagens serão matrizes de inteiros entre $0$ (preto) e $255$ (branco). O gradiente morfológico é uma operação muito simples que realça as bordas de uma imagem, seja ela binária ou em níveis de cinza. Para cada ponto da imagem de saída, analisamos os valores presentes em uma vizinhança de tamanho $3\times 3$ ao seu redor. O valor de saída é a diferença entre o pixel de maior valor e o de menor valor.

O código abaixo calcula o gradiente morfológico sequencialmente e mostra a imagem resultado.

In [1]:
%matplotlib inline
import numpy as np
import scipy as sp
import scipy.ndimage
import time
import matplotlib.pyplot as plt

def morpho_gradient(img, out, p):
    pi, pj = p
    minv = 255
    maxv = 0
    
    for i in range(-1, 2):
        for j in range(-1, 2):
            if img[pi+i, pj+j] < minv:
                minv = img[pi+i, pj+j]
            if img[pi+i, pj+j] > maxv:
                maxv = img[pi+i, pj+j]
    out[pi, pj] = maxv - minv

def versao_serial(img):
    out = np.zeros(img.shape, img.dtype)
    for i in range(1, img.shape[0]-1):
        for j in range(1, img.shape[1]-1):
            morpho_gradient(img, out, (i, j))
    return out

img = sp.misc.lena()
start = time.clock()
out = versao_serial(img)
print('Versao serial:', time.clock() - start)
plt.figure(figsize=(10, 10))
plt.subplot(121)
plt.gray()
plt.imshow(img)
plt.subplot(122)
plt.gray()
plt.imshow(out)
Versao serial: 2.75
Out[1]:
<matplotlib.image.AxesImage at 0x7f52d78bbda0>

Seguindo a nossa receita anterior, podemos paralelizar os dois for da versao_serial. Para facilitar a passagem de argumentos para nossa função iremos usar os módulos functools e itertools.

Como vimos anteriormente, a função chamada pelo Pool.map só recebe um argumento. Podemos usar a função partial do módulo functools para criar uma versão de morpho_gradient cujos parâmetros img e out são fixos. A função chama-se partial pois ela aplica parcialmente a função.

Para criar a lista de pontos a serem processados usamos o módulo itertools, que permite criar diversos tipos de iteradores. No nosso caso, queremos o produto cartesiano entre os pontos horizontais e verticais da imagem. A função itertools.product realiza exatamente este trabalho. Este módulo contém diversas funções úteis e vale a pena uma olhada na documentação para conhecê-lo melhor.

A nossa primeira tentativa de paralelização está abaixo.

In [2]:
import multiprocessing as mp
import functools
import itertools

def versao_paralela(img):
    out = np.zeros(img.shape, img.dtype)
    p = mp.Pool()
    pi = range(1, img.shape[0]-1)
    pj = range(1, img.shape[1]-1)
    point_list = itertools.product(pi, pj)
    par = functools.partial(morpho_gradient, img, out)
    p.map(par, point_list)
    p.close()
    p.join()
    
    return out

start = time.clock()
out2 = versao_paralela(img)
print('Versao paralela:', time.clock() - start)
plt.figure(figsize=(10, 10))
plt.subplot(121)
plt.gray()
plt.imshow(img)
plt.subplot(122)
plt.gray()
plt.imshow(out)
Versao paralela: 1.3299999999999992
Out[2]:
<matplotlib.image.AxesImage at 0x7f5328a5ce80>

Surpreendentemente (ou não), esta versão não é muito mais rápida que a serial e dependendo do ambiente de execução pode ser até mais lenta. A razão disto é que, da mesma maneira do exemplo anterior, fazemos uma cópia de img e out para cada processo do Pool. No exemplo anterior, cada elemento da lista era diferente e esta cópia realmente é necessária, mas neste caso todos os processos leem img e cada processo escreve em uma posição $(i,j)$ diferente de out.

Objetos "normais" de Python, como listas, dicionários e até arrays do Numpy, não são "compartilháveis" automaticamente. Todo objeto contém informações específicas do processo em que foram criadas e que podem se tornar inválidas se manipuladas paralelamente por outros processos. É preciso, portanto, utilizar versões destes objetos que suportem a utilização por mais de um processo.

O módulo multiprocessing.sharedctypes contém, basicamente, dois tipos de objetos com nomes bastante descritivos: Array e Value. Quando criamos um objeto destes tipos alocamos sua memória em uma área especial que pode ser acessada por mais de um processo ao mesmo tempo. Porém, é necessário especificar o tipo dos objetos armazenados e, para o Array, o número de elementos desejados. Também é necessário deletá-los manualmente. Ou seja, é praticamente uma chamada às funções malloc()/free() usadas em C.

A utilização de memória compartilhada deixa o código um pouco mais longo, mas as mudanças são mais mecânicas que complexas. Primeiramente, precisamos copiar nossas imagens para a memória compartilhada (variáveis img_shared e out_shared) e criar uma nova função morpho_grad_map que cria arrays do numpy utilizando a memória compartilhada que alocamos para então chamar morpho_gradient. Note que usamos variáveis globais chamadas img e out em morpho_grad_map. Essas variáveis globais são criadas na função init_pool, que é chamada ao inicializar cada processo do Pool e que recebe as nossas variáveis compartilhadas img_shared e out_shared. As variáveis globais só existem nos processo do Pool e não irão vazar para processo principal que estamos executando. Desta maneira, todos os processos contém uma referência à nossa área de memória compartilhada e evitamos copiar dados entre os processos.

Novamente, o código da versão paralela2 está abaixo.

In [3]:
import multiprocessing.sharedctypes
from ctypes import c_int32

def morpho_grad_map(p):
    img_np = np.frombuffer(img, dtype=np.int32).reshape(shape)
    out_np = np.frombuffer(out, dtype=np.int32).reshape(shape)
    morpho_gradient(img_np, out_np, p)

def init_pool(img_shared, out_shared, shape_):
    global img, out, shape
    img = img_shared
    out = out_shared
    shape = shape_

def versao_paralela2(img):
    out_shared = mp.sharedctypes.RawArray(c_int32, np.zeros(img.shape, np.int32).flat)
    img_shared = mp.sharedctypes.RawArray(c_int32, img.flat)

    p = mp.Pool(initializer=init_pool, initargs=(img_shared, out_shared, img.shape))

    pi = range(1, img.shape[0]-1)
    pj = range(1, img.shape[1]-1)
    point_list = itertools.product(pi, pj)
    p.map(morpho_grad_map, point_list)
    p.close()
    p.join()
    
    out_np = np.zeros(img.shape, np.int32)
    out_np.flat[:] = out_shared # usar um view aqui para copiar os dados
    
    del img_shared
    del out_shared

    return out_np


start = time.clock()
out2 = versao_paralela2(img)
print('Versao paralela 2:', time.clock() - start)
plt.figure(figsize=(10, 10))
plt.subplot(121)
plt.gray()
plt.imshow(img)
plt.subplot(122)
plt.gray()
plt.imshow(out)
Versao paralela 2: 0.40000000000000036
Out[3]:
<matplotlib.image.AxesImage at 0x7f5328548d30>

Os ganhos de desempenho desta versão são significativos e mostram a capacidade do módulo multiprocessing em acelerar nosso código. Os benefícios da computação paralela, porém, não são exatamente proporcionais ao número de processadores usados e muitas vezes um código paralelo que funciona bem em um sistema pode não ser tão eficiente em outros. Por exemplo, ao rodar este post em um computador com 4 processadores a diferença entre as duas versões paralelas foi bem menor, pois o tempo gasto para criar 4 cópias é menor que o tempo necessário para criar 24.

É importante notar, porém, que usar memória compartilhada é um trabalho mais delicado. No nosso caso cada processo usou uma matriz somente leitura e escreveu em uma posição diferente da matriz de saída. Esta condição não é sempre verdadeira. Nos próximos textos veremos como utilizar travas para compartilhar objetos em situações em que pode haver conflito entre os processos.

sábado, 17 de outubro de 2015

Processamento paralelo usando Python - parte 1

Atualização: Veja também a parte 2 da série sobre processamento paralelo.

Uma das coisas mais chatas (e angustiantes) de se fazer pesquisa é esperar aquele experimento de 3 horas acabar só para você ver que tinha um erro e você precisa executar tudo de novo (e perder mais 3 horas esperando). E o pior de tudo é que muitas vezes não estamos nem usando o máximo da capacidade de nossa máquina. Neste tutorial iremos ver uma técnica simples de programação em Python para paralelizar código que realiza tarefas independentes.

Frequentemente, nossos programas são muito similares ao trecho de código abaixo.

In [1]:
def funcao_demorada(el):
    pass
    # funcao demorada

def cria_lista(n):
    return []
    
long_list = cria_lista(1000000)
resultados = []
for el in long_list:
    resultados.append(funcao_demorada(el))
# processa resultados

Neste caso, processamos os elementos de uma lista um por vez e (idealmente) a função run_task não muda o valor de nenhuma variável global nem os atributos de nenhum dos objetos passados para a função. Como as execuções de run_task são independentes uma da outra, iremos paralelizar o loop da linha #2 e executar vários run_task ao mesmo tempo usando o módulo multiprocessing.

Nos próximos exemplos, funcao_demorada recebe uma matrizes A e B e faz uma série de operações de matrizes e cria_lista gera matrizes quadradas aleatoriamente.

Utilizaremos a classe Pool para criar um conjunto de processos que irão executar nosso código. Uma instância de Pool possui o método map, que faz o papel do for da versão serial. Ao chamar o map, o Pool seleciona um processo livre e o configura para chamar a nossa função passando como argumento um elemento de nossa lista (no nosso caso, duas matrizes). Existem alguns pontos importantes a serem considereados quando programamos usando multiprocessing:

  1. todo elemento da lista de argumentos deve ser picklable;
  2. a funcao_demorada só pode usar variáveis de sua lista de argumentos;
  3. modificações feitas nas variáveis (globais e locais) não serão propagadas para o processo original.

Além disto, é recomendado que o ponto de entrada do seu programa esteja em um if __name__ == "__main__": para evitar problemas ao rodar em Windows.

Como comparação, iremos rodar funcao_demorada 500 vezes para matrizes 250x250 e medir o tempo de uma versão serial e de uma versão paralelizada.

In [2]:
import numpy as np
import time
import multiprocessing as mp


def funcao_demorada(el):
    # isto não faz sentido algum...
    A, B = el
    A += B
    A *= B
    A = np.dot(A, B)
    A = np.linalg.inv(A)
    A = np.dot(B, A)
    B = np.dot(A, B)
    A = np.linalg.inv(A)
    A = np.dot(A, B)
    B = np.dot(B, A)
    A = np.dot(A, B)
    return A

def cria_lista(n, sz=(250, 250)):
    return [(np.random.rand(sz[0], sz[1]), np.random.rand(sz[0], sz[1])) for i in range(n)]

def versao_serial():
    lista = cria_lista(500)
    resultados = []
    for el in lista:
        resultados.append(funcao_demorada(el))

def versao_paralela():
    p = mp.Pool(mp.cpu_count())
    lista = cria_lista(500)
    resultados = p.map(funcao_demorada, lista)

def executa_varias_vezes(func, n):
    tempos = []
    for i in range(n):
        start = time.clock()
        func()
        tempos.append(time.clock() - start)
    return tempos

if __name__ == "__main__":
    tserial = executa_varias_vezes(versao_serial, 10)
    tparalela = executa_varias_vezes(versao_paralela, 10)
    
    print('Serial:', sum(tserial)/10)
    print('Paralela:', sum(tparalela)/10)
Serial: 18.259999999999998
Paralela: 4.306

O código acima foi executado em uma máquina de grande porte com 24 processadores (rodando Ubuntu) e os ganhos são significativos: a versão serial demora em média 18,25 segundos, enquanto a versão paralela demora apenas 4,3 segundos.

É importante notar, porém, que o programa não rodou 24x vezes mais rápido só por usar 24 processadores. O tempo de criar todos os processos auxiliares e de copiar os objetos necessários para cada processo pode ser significativo dependendo de cada caso. Existem também diversas diferenças entre Sistemas Operacionais que tornam o módulo multiprocessing mais ou menos efetivo. O que funciona bem em um Linux pode não funcionar igualmente bem em ambientes Windows. De qualquer maneira, mesmo máquinas com 2 ou 4 processadores podem obter aumentos de velocidade significativos usando multiprocessing.

Este exemplo simples já mostra o quão interessante pode ser a utilização de paralelismo em Python. Nos próximos posts mostrarei mais da API de processamento paralelo em Python e quanto podemos ganhar em eficiência utilizando-a.

Se este texto lhe for útil, não deixe de escrever um comentário contando a sua experiência com processamento paralelo e Python ;)

segunda-feira, 28 de setembro de 2015

Condições de otimalidade KKT

\[ \newcommand{\I}[0]{\mathcal{I}} \newcommand{\E}[0]{\mathcal{E}} \newcommand{\R}[0]{\mathbb{R}} \renewcommand{\L}[0]{\mathcal{L}} \newcommand{\eq}[2]{\begin{equation} #1 \label{#2} \end{equation}} \]
Técnicas de otimização são usadas em diversos contextos para encontrar a melhor solução de um problema, que pode ou não envolver diversas restrições nas soluções aceitas. A qualidade de uma solução é medida por uma função objetivo \(f(x)\) e toda solução deve obedecer uma série de restrições \(c_i(x)\). Neste texto apresentamos alguns resultados que levam às Condições de Otimalidade de primeira ordem, também conhecidas como condições KKT Formalmente, um problema de otimização pode ser expresso, de forma geral, como
\[ \begin{equation} \begin{cases} \textrm{minimize }f(x) \\ \textrm{ s.t.}\\ c_j(x) = 0, j \in \E \\ c_i(x) \geq 0, i \in \I \end{cases} \label{eq:formulation} \end{equation} \]
Primeiramente, estamos interessados somente nos pontos que satisfazem todas as restrições. O conjunto de pontos $\Omega$ que satisfaz todas as restrições é chamado conjunto viável.
\[\eq{ \omega = \{x \in \R^n : c_i(x) = 0, i \in \I ; c_j(x) \geq 0, j \in \E \} }{eq:viable-set} \]
Encontrar o mínimo global de \(f(x)\) é difícil mesmo quando não existem restrições. Portanto, estudamos o conjunto de soluções locais \(f(x)\). Sob certas condições (como no caso de Programação Linear ), mínimos locais também podem ser mínimos globais de \(f(x)\).
  • Uma solução \(x^* \in \omega\) é uma solução local se existe uma vizinhança \(\mathcal{N}\) de \(x^*\) tal que\(f(x^*)\geq f(x)\) para todo \(x \in \omega \cap \mathcal{N}\).
    Uma solução \(x^* \in \omega\) é uma solução local isolada se existe uma vizinhança \(\mathcal{N}\) de \(x^*\) tal que \(x^*\) é a única solução local. Isso implica \(f(x) > f(x^*), x \in \mathcal{N} \cap \omega\).
  • O conjunto ativo \(\mathcal{A}(x) = \{ i : c_i(x) = 0, i \in \I \cup \E \}\) contém todas as restrições satisfeitas com igualdade.
  • A suavidade de \(f(x)\) e das restrições \(c_i(x)\) fornece uma maneira de percorrer o conjunto viável \(\Omega\). Dada uma solução \(x \in \Omega\), suavidade garante que podemos buscar na vizinhança de \(x\) por uma solução \(x + \epsilon d\) que melhore a função objetivo e seja viável. A partir de agora tanto \(f(x)\) quanto as restrições \(c_i(x), i \in \I \cup \E\) são suaves e diferenciáveis uma vez.

    Aproximação de mínimos locais

    Dada uma solução \(x \in \R^n\) e um pequeno "passo" \(s \in \R^n\), iremos analisar a vizinhança de x para obter intuições de como navegar de solução em solução de modo a sempre diminuir o valor da função objetivo \(f(x)\). Esta análise serve como motivação para as condições KKT.
    Primeiramente, como queremos minimizar \(f(x)\), precisamos que \(f(x + s) - f(x) < 0\). Aproximando \(f(x+s)\) à primeira ordem temos
    \[ f(x + s) = f(x) + s^T\nabla f(X) < f(x) \Rightarrow s^T\nabla f(x) < 0. \]
    Como todas as restrições de igualdade devem ser mantidas, também queremos que \(c_j(x + s) = 0\). Fazendo novamente uma aproximação de primeira ordem, temos que
    \[ 0 = c_j(x + s) = c_j(x) + s^T\nabla c_j(x) \Leftarrow s^T \nabla c_j(x) = 0. \]
    Também é necessário manter viabilidade nas inequações (i.e. \(c_i(x) \geq 0, i \in \I\)), o que significa \(c_i(x+s) \geq 0\). Aproximando à primeira ordem,
    \[ c_i(x + s) = c_i(x) + s^T \nabla c_i(x) \geq 0. \]
    Caso I: A restrição não está ativa e \(c_i(x) > 0\). Então qualquer passo\(s\)pequeno o suficiente satisfaz a equação acima.
    Case II:A restrição está ativa e \(c_i(x) = 0\). Então as equações se tornam
    \[ s^T \nabla c_i(x) \geq 0 \textrm{ and } s^T\nabla f(x) < 0. \]
    Como buscamos por um mínimo local, queremos que a interseção de ambas as condições seja vazia. Isto é alcançado quando \(\nabla c_i(x)\) e \(\nabla f(x)\) possuem a mesma direção, (\(\nabla f(x) = \lambda_i \nabla c_i(x), \lambda_i > 0\)).
    Estas duas condições podem ser sumarizadas usando a função Lagrangiano definida abaixo:
    \[ \L(x, \lambda) = f(x) - \sum_j^\E \lambda_j c_j(x) - \sum_i^\I \lambda_i c_i(x) \]
    Em um mínimo local, é necessário que \[ \nabla \L_x(x^*, \lambda^*) = \nabla f(x) - \sum_j^\E \lambda_j c_j(x) - \sum_i^\I \lambda_i c_i(x) = 0 , \textrm{ for some } \lambda^* \geq 0 \] e \[ \lambda^*_i c_i(x) = 0, i \in \I. \]
    Esta última condição é uma condição de complementaridade: para todas inequações ativas precisamos que \(\lambda^*_i > 0\). Para as outras, ao definir seus \(\lambda^*_i = 0\) correspondentes os eliminamos do Lagrangiano e \(\nabla \L(x^*, \lambda^*) = f(x^*)\). Estas condições são compiladas nas Condições de primeira ordem KKT.

    Condições de otimalidade KKT

    As condições KKT sumarizam as condições necessárias de primeira ordem para que um ponto \(x^*\) seja um mínimo local. Estas condições estão associadas com os gradientes de \(f(x)\) e \(c_i(x)\). As condições são necessárias, mas não suficientes
    Dado um mínimo local \(x^*\) (tal que os gradientes de suas restrições ativas seja LI), existe \(\lambda^* \geq 0\) tal que:
    \[ \begin{equation} \label{eq:kkt1} \nabla_x \L(x^*, \lambda^*) = 0, \end{equation} \] \[ \begin{equation} c_i(x) = 0, i \in \E, \end{equation} \] \[ \begin{equation} c_i(x) \geq 0, i \in \I, \end{equation} \] \[ \begin{equation} \lambda_i^* \geq 0, i \in \I, \end{equation} \] \[ \begin{equation} \lambda^*_i c_i(x) = 0, i \in \E \cup \I. \end{equation} \]

    Referências:

    Jorge Nocedal, Stephen J Wright, "Numerical Optimization", Springer-Verlag New York, 2006.

    quarta-feira, 18 de março de 2015

    Lema de Farkas

    Lema de Farkas

    Author: Igor dos Santos Montagner
    Date: 25/02/2015

    O lema de Farkas é parte importante da Teoria de Otimização Convexa pois apresenta uma "alternativa": dados um cone convexo e um ponto \(g \in R^n\), ou \(g\) pertence ao cone ou existe um hiperplano que os separa. Construindo um cone adequado podemos usar este resultado para provar as condições de otimalidade de primeira ordem KKT (que serão alvo do próximo texto).

    Definição:

    Sejam \(B \in R^{n\times m}, C \in R^{n \times p}\) duas matrizes que contém vetores de \(R^n\) em suas colunas e \(y \in R^m, w \in R^p\). O conjunto \(K\) definido abaixo é um cone convexo.

    \[ K = \{ By + Cw | y \geq 0 \} \]

    Lema:

    Dado um ponto \(g \in R^n\) e um cone convexo \(K\), então somente uma das alternativas abaixo é verdadeira:

    1. \(g \in K\)
    2. Existe \(d \in R^n\) tal que:
      1. \(d^T g < 0\)
      2. \(B^T d \geq 0\)
      3. \(C^T d = 0\)

    O interessante da prova do lema é que ela mostra como construir a direção \(d\) que satisfaz estas condições.

    Prova:

    Primeiro iremos provar que as duas alternativas não podem ocorrer juntas.

    Suponha que as duas alternativas valham, ou seja, que \(g \in K\) e existe \(d \in R^n\) que satisfaz as condições acima.

    Primeiramente, \(g = By + Cw\) para algum \(y \in R^m, y \geq 0\) e \(w \in R^n\). Logo, \(0 > g^T d = (By + Cw)^T d = y^T B^T d + w^T C^T d\). Porém, como \(C^T d = 0\), isto é igual a

    \[ 0 > g^T d = (By + Cw)^T d = y^T B^T d \geq 0, \]

    pois \(B^T d \geq 0\) e \(y \geq 0\). Logo, chegamos em uma contradição e as duas condições são disjuntas.

    Vamos mostrar agora como construir a direção \(d\) no caso em que \(g \notin K\). Primeiramente, seja \(\hat{s} \in K\) o elemento do cone com menor distância euclidiana até \(g\). Logo, \(\hat{s}\) é a solução do seguinte problema de minimização:

    \[ min_{s \in K} ||s - g||_2^2 \]

    Definimos então \(d = \hat{s}-g\) e, consequentemente, \(g=\hat{s}-d\).

    A direção d em relação a g e s
    Figura: A direção d em relação a \(g\) e \(s\).

    Condição 1:

    Como o valor mínimo da função acima encontra-se em \(\hat{s}\), sabemos que a função \(|| t \hat{s} - g||_2^2\) possui mínimo em \(t=1\), pois \(t \hat{s} \in K, \forall t \geq 0\). Logo,

    \[ \frac{||t\hat{s} - g||_2^2}{dt} |_{t=1} = 0 \\ \frac{t^2 \hat{s}^T \hat{s} - 2t\hat{s}^Tg + ||g||_2^2}{dt} |_{t=1} = 0 \\ 2t\hat{s}^T \hat{s} - 2\hat{s}^Tg |_{t=1} = 0 \\ 2\hat{s}^T(t\hat{s} - g)|_{t=1} = 0 \\ \hat{s}^T(\hat{s} - g) = 0 \]

    Logo,

    \[ d^Tg = d^T(\hat{s}-d) = d^T\hat{s} - d^T d = \\ \hat{s}^T(s - d) - ||d||_2^2 = -||d||_2^2 < 0 \]

    Condições 2 e 3:

    Para mostrar que \(d\) satisfaz as condições 2 e 3 vamos primeiramente analisar o que acontece quando consideramos, em vez de \(\hat{s}\), \(\hat{s} + \theta(\hat{s} - s), s \in K\).

    \[ ||\hat{s} + \theta(\hat{s} - s) - g ||_2^2 \geq ||\hat{s} - g||_2^2 \\ ||\hat{s} - g||_2^2 +2\theta(\hat{s} - g)^T(\hat{s} - s) + \theta^2||\hat{s} - s||_2^2 \geq ||\hat{s} - g||_2^2 \\ 2\theta(\hat{s} - g)^T(\hat{s} - s) + \theta^2||\hat{s} - s||_2^2 \geq 0 \\ 2(\hat{s} - g)^T(\hat{s} - s) + \theta||\hat{s} - s||_2^2 \geq 0 \\ \theta \rightarrow 0 \\ 2(\hat{s} - g)^T(\hat{s} - s) \geq 0 \\ 2(\hat{s} - s)^T(\hat{s} - g) = \hat{s}^T(\hat{s} - g) - s^T(\hat{s} - g) \geq 0 \\ s^T(\hat{s} - g) = d^Ts \geq 0 \forall s \in K \\ \]

    Logo, \(d^Ts = d^T(By + Cw) \geq 0, y \geq 0\). Isto vale para todos os \(s \in K\) e, em especial para os casos abaixo:

    1. \(y = 0\), \(d^TCw \geq 0 \rightarrow d^TC = C^T d = 0\)
    2. \(w = 0\), \(d^TBy \geq 0 \rightarrow B^Td \geq 0\), pois \(y \geq 0\).

    Portanto, a direção \(d = \hat{s} - g\) satisfaz as propriedades do lema se \(g \notin K\) e tudo está provado :)


    quarta-feira, 11 de março de 2015

    Teoria de W-operadores - Operações entre imagens

    Teoria de W-operadores - Operações entre imagens

    \[ \newcommand{\E}[0]{\mathbb{E}} \newcommand{\Z}[0]{\mathbb{Z}} \newcommand{\L}[0]{\mathcal{L}} \newcommand{\P}[0]{\mathcal{P}} \]

    Imagens podem ser interpretadas como reticulados e operações entre imagens como operações entre reticulados. Neste texto procuro apresentar a relação os diversos aspectos de teoria de reticulados apresentados anteriormente e o processamento de imagens digitais. Em especial, mostrarei que o operador intervalo produz resultados visualmente intuitivos e como ele pode ser composto para produzir resultados mais complexos.

    Imagens como elementos de reticulados

    Dados um domínio \(\E \subset \Z^2\) e um intervalo \(K = [0, \dots, k]\), uma imagem em \(k+1\) níveis de cinza pode ser interpretada como um elemento \( I \in Fun[\E, K]\). Em outras palavras, \(I\) atribui para cada ponto do domínio um nível de cinza.

    O conjunto \(Fun[\E, K]\) é um reticulado com a operação \(\preccurlyeq\) definida abaixo.

    \[ f \preccurlyeq g \Leftrightarrow f(x) \leq g(x) \forall x \in \E \]

    De maneira equivalente, uma imagem binária \(B \in Fun[\E, [0,1]]\) também pode ser interpretada como um subconjunto de \(B' \subseteq \E\), onde \(p \in B' \Leftrightarrow B(p) = 1\). Seguindo a interpretação de conjuntos, o reticulado \((\P(\E), \subseteq)\) das imagens binárias é formado pelo conjunto de todos subconjuntos de \(\E\), \(\P(\E)\), e pela relação de inclusão usual em conjuntos \(\subseteq\).

    Note que a noção de ordem (parcial) entre imagens é fundamental para a interpretação de imagens como reticulados.

    Muitas vezes é interessante selecionar um pequeno recorte da imagem, selecionando os pixels que pertencem à vizinhança, chamada de janela, de um pixel. Estas imagens são chamadas de imagens-janela e são elementos em \(Fun[W, K]\), onde \(W \subset \Z^2\) é um conjunto de pontos que define a vizinhança considerada. A imagem janela \(I_z^{(W)}\) obtida ao centralizar a janela \(W\) no pixel \(z \in \Z^2\) na imagem \(I\) é dada por

    \[ I^{(W)}_z(p) = I(z + p) \forall p \in W.\]

    Processamento de imagens como uma operação de reticulados

    Um operador de imagens é uma função que transforma uma imagem em outra imagem diferente. Logo, um operador \(\Psi\) é um elemento de \(Fun[ Fun[\E, K], Fun[\E, K] ]\) e é tanto uma função entre reticulados quanto um elemento no reticulado dos operadores. Um classe de especial interesse entre os operadores de imagens é a classe dos \(W-\)operadores.

    Um operador \(\Psi\) é u \(W-\)operador se ele possui as duas seguintes propriedades:

    1. invariante à translação: \(\Psi(t(I, p)) = t(\Psi(I), p)\), onde \(t(I, p)\) representa a translação da imagem \(I\) por \(p \in \Z^2\);
    2. localmente definido: existe uma janela \(W \in \Z^2\) tal que \(\Psi(I)(p) = \Psi(I_p^{(W)}), \forall I \in Fun[\E, K], p \in \Z^2\).

    Todo \(W-\)operador \(\Psi\) pode ser unicamente caracterizado por uma função \(\psi \in Fun[ Fun[W, K], K]\), de maneira que

    \[ \Psi(I)(p) = \psi(I_z^{(W)}) \forall I \in Fun[\E, L], p \in \Z^2.\]

    Como \(\psi\) é um operador entre um reticulado (\(Fun[W, K]\)) e uma cadeia (\(K\)), \(\psi\) pode ser expresso utilizando a decomposição canônica:

    \[\psi(I)(p) = \sum_{y=0}^m \vee \{ \lambda[A,B](I_p^{(W)} : [A,B] \in \textbf{B}_\psi(k) \}.\]

    Lembrando que o operador \(\lambda[A,B] \in Fun[ Fun[W, K], \{0, 1\}]\) tem a seguinte forma:

    \[ \lambda_{A,B}(I) = \begin{cases} 1 \textrm{, if } I \in [A,B] \\ 0 \textrm{, otherwise. } \end{cases} \]

    Por um lado, o operador de imagens \(\Psi\), não importa tão complexo ele seja, pode ser decomposto em uma série de operadores de imagens fundamentais. Por outro lado, qualquer operação de imagens pode ser desenvolvida ou aprendida se for possível determinar quais intervalos devem ser usados para representar o operador.

    O operador intervalo em imagens

    Nesta seção mostrarei como utilizar operadores intervalo para construir um operador que detecta bordas em imagens binárias como a imagem abaixo.

    O operador intervalo é parametrizado por dois extremos \(A\) e \(B\). Todos os pontos cuja imagem janela contem o primeiro e está contida no segundo extremo terão como saída um pixel branco. Para identificar os cantos superiores direitos podemos utilizar os seguintes intervalos:

    De maneira similar, os seguintes intervalos podem ser utilizados para detectar os extremos verticais esquerdos das formas:

    O mesmo raciocínio pode ser aplicado para os outros cantos e extremos verticais e horizontais, resultando no seguinte conjunto de intervalos:

    Veja abaixo o resultado da união dos operadores intervalo acima na imagem de exemplo.

    Note que este conjunto de intervalos não é capaz de detectar todas as bordas diagonais, pois elas possuem uma variação maior que as bordas em 90 graus.


    segunda-feira, 2 de março de 2015

    Teoria de W-operadores - Reticulados e Operadores entre reticulados

    \[ \newcommand{\L}[0]{\mathcal{L}} \newcommand{\P}[0]{\mathcal{P}} \]
    A revisão de teoria continua, desta vez com aspectos fundamentais do treinamento de \(W-\)operadores. Neste documento irei apresentar algumas definições importantes de reticulados e operadores entre reticulados.

    Conjuntos parcialmente ordenados e Reticulados completos

    Seja \(\mathcal{L}\) um conjunto e \(\leq\) uma relação binária entre elementos de \(\mathcal{L}\). Se a relação \(\leq\) for:
    1. reflexiva (\(x\leq x \forall x \in \mathcal{L}\));
    2. anti-simétrica (\(x \leq y\) e \(y \leq x\) implica \(x = y\));
    3. transitiva (\(x \leq y\) e \(y \leq z\) implica \(x \leq x\)),
    então \((\mathcal{L}, \leq)\) é chamado de conjunto parcialmente ordenado ou poset (de partially ordered set). A relação \(\leq\) é chamada de relação de ordem parcial (pois não requer que todos elementos estejam relacionados).
    Sejam \(L, U \in \mathcal{L}\) e \(\Xi \subseteq \mathcal{L}\), então \(L\) e \(U\) são chamados de limite superior e limite inferior de \(\Xi\) se \(X \leq U \forall X \in \Xi\) e \(L \leq X, \forall X \in \Xi\), respectivamente. O menor limitante superior de \(\Xi\), se ele existir, é chamado de supremo. Da mesma maneira, o maior limitante inferior de \(\Xi\) é chamado de ínfimo. Tanto o ínfimo quanto o supremo são únicos, se existirem. O ínfimo e o supremo entre dois elementos \(X\) e \(Y\) são denotados \(X \wedge Y\) e \(X \vee Y\).
    Um subconjunto \(\Xi \subseteq \mathcal{L}\) é um intervalo se e somente se existem dois elementos \(A,B \in \mathcal{L}\) tal que
    \[ A \leq X \leq B \Leftrightarrow X \in \Xi, \forall X \in \mathcal{L}.\]
    Intervalos são denotados por \([A,B]\), sendo que \(A\) é a extremidade esquerda e \(B\) a extremidade direita do intervalo. O conjunto de intervalors \(\{[A,B] : A,B \in \L, A \leq B\}\) junto com a operação \(\preccurlyeq\) é um reticulado completo.
    \[ [A,B] \preccurlyeq [A',B'] \Leftrightarrow A \leq A' \textrm{ e } B \leq B' \]
    O conjunto \(Max(\{[A,B]\})\) contém todos os intervalos maximais.
    Um poset \(\mathcal{L}\) é um reticulado completo se todo subconjunto de \(\L\) possui um ínfimo e um supremo. Reticulados completos sempre possuem um máximo \(I\) e um mínimo \(O\). Quando todos os elementos de um reticulado completo são comparáveis ele é chamado de cadeia.

    Operadores entre reticulados

    Dados reticulados \(\L_1\) e \(\L_2\), o conjunto \(Fun[\L_1,\L_2]\) contém todas as funções de \(\L_1\) a \(\L_2\). Elementos deste conjunto são chamados de operadores ou mapeamentos. Elementos de \(\L_1\) serão denotados \(A, B\) e \(X\) e elementos de \(\L_2\) serão denotados \(V\) e \(Y\). Operadores serão denotados por letras gregas minúsculas \(\alpha, \beta, \dots\).
    Defina \(\lambda_{A,B}\) um operador (chamado de sup-gerador) de \(\L_1\) para \(\{0, 1\}\) da seguinte maneira:
    \[ \lambda_{A,B}(X) = \begin{cases} 1 \textrm{, if } X \in [A,B] \\ 0 \textrm{, otherwise } \end{cases} \]
    O kernel \(K_\psi \in Fun[L_2, \P(L_1)]\) de um operador \(\psi \in Fun[\L_1, \L_2]\) é dado por
    \[ K_\psi(Y) = \{X \in \L_1 : Y \leq \psi(X) \}. \]
    Essencialmente, \(K_\psi\) contém, para cada elemento \(Y\) todos os elementos de \(L_1\) que são maiores ou iguais a \(Y\) após a aplicação de \(\psi\).

    Teorema (Decomposição Canônica):

    Todo operador \(\psi \in Fun[\L_1, \L_2]\) pode ser decomposto em função de operadores sup-geradores:
    \[ \psi(X) = \vee \{ Y \in \L_2 : \vee\{\lambda_{A,B}(X) : [A,B] \in K_\psi(Y) \} = 1\}, \forall X \in \L_1 \]
    A representação de um operador utilizando a decomposição acima é redundante, pois se \(X \in [A,B]\) e \(X \in [A',B']\) tal que \(A\leq A'\) e \(B\leq B'\) então o intervalo \([A,B]\) é redundante. Desta forma, definimos
    \[ \textbf{B}_\psi(Y) = Max(K_\psi(Y)) = \{[A,B] : \nexists [A',B'] \textrm { s.t. } [A,B] \preccurlyeq [A', B']\} \]
    e o teorema acima se torna
    \[ \psi(X) = \vee \{ Y \in \L_2 : \vee\{\lambda_{A,B}(X) : [A,B] \in \textbf{B}_\psi(Y) \} = 1\}, \forall X \in \L_1 \]
    e a representação de \(\psi\) em termos de \(\textbf{B}_\psi\) é chamada de decomposição por um conjunto de operadores sup-geradores.

    Caso específico: Operadores entre reticulados e cadeias

    Suponha que \(\L_2 = M = [0, m]\). Então o conjunto \(K = \{K_\psi(y) : y \in M\}\) também é uma cadeia em \((\P(\L_1), \subseteq)\) para todo \(\psi \in Fun[\L_1, M]\). Logo, o teorema da decomposição canônica se torna:
    \[ \psi(X) = \sum_{y=1}^m \vee \{\lambda_{A,B}(X) : [A,B] \in \textbf{B}_\psi(y)\} \]
    Se \(m=1\), o teorema se torna:
    \[ \psi(X) = \vee \{\lambda_{A,B}(X) : [A,B] \in \textbf{B}_\psi(1)\} \]

    Caso específico II: Imagens

    Este outro texto explica o caso específico que interpreta imagens como reticulados e contextualiza a decomposição canônica como uma composição de operações de imagens.

    sábado, 7 de fevereiro de 2015

    Estudos de otimização: Programação linear

    Um dos tópicos de revisão de otimização que estou estudando é o de Programação Linear. Neste tipo de problema de otimização tanto a função objetivo quanto as restrições são funções lineares das variáveis.

    Apesar de problemas deste tipo poderem ser expressos de diversas maneiras, a maneira abaixo, chamada de forma canônica, é frequentemente utilizada por materiais da área (e também por pacotes computacionais).

    \[ \begin{cases} \textrm{minimize} f(x) = c^Tx \textrm{ s.a.}\\ Ax = b\\ x \geq 0 \end{cases} \]

    Conversão para a forma canônica

    Seguem abaixo alguns ajustes que podem ser feitos para converter um problema linear qualquer para sua forma canônica:

    1. Se problema for \(maximize c^Tx\), converta-o para \(minimize -c^Tx\);
    2. Se houverem restrições de \(\geq\) ou \(\leq\), adicione uma variável \(s_i\) para cada restrição e converta-as em restrições de igualdade da seguinte maneira: se \(A_ix \geq b\), então \(A_ix -s_i = b\); se \(A_ix \leq b\), então \(A_ix + s_i = b\).
    3. Se houverem variáveis que podem ser negativas (sem restrição \(x_i \geq 0\)), faça a substituição \(x_i = x_i^{(+)} - x_i^{(-)}, x_i^{(+)} \geq 0, x_i^{(-)} \geq 0\).

    As transformações (1) e (3) são relativamente fáceis de serem entendidas. No caso da transformação (2), é possível atribuir a \(s_i\) qualquer valor (pois s_i não faz parte da função objetivo). Desta maneira, se durante o proceso de otimização \(A_ix = b\), então basta que \(s_i = 0\).

    O problema Dual

    Todo problema de minimização (primal) em PL possui um problema de maximização (dual) relacionado de modo que ambos possuem o mesmo valor de função objetivo e é trivial encontrar a resposta de um a partir do outro.

    Dado o problema linear na forma canônica acima, o problema dual consiste em:

    \[ \begin{cases} \textrm{maximize} f(y) = b^Ty \textrm{ s.a.}\\ A^Ty \leq c\\ \end{cases} \]

    A partir da solução ótima \(y^*\) do problema dual, é possível obter a solução ótima \(x^*\) do problema primal.

    Seja \(u = c - A^Ty\). O vetor \(u\) contem a "distância" que cada restrição dual está da igualdade. Uma restrição somente restringe a função objetivo se está "saturada", ou seja, \(u_i = 0\). Desta maneira, podemos ignorar toda coluna de \(A\) tal que \(u_j \neq 0\), pois esta restrição não está ativa na solução do dual.Se o problem dual tiver solução única o número de componentes não zero de \(u\) será igual ao número de linhas de \(A\) (número de restrições primais). Logo, podemos obter \(x\) resolvendo o seguinte sistema linear:

    \[ A[:,u_i = 0]x^* = b \].

    As posições de x^* tal que \(u_i \neq 0\) são completadas com \(x^*_i = 0\), pois todas as restrições de \(A\) já são satisfeitas com igualdade pelas outras componentes de \(x^*_i\).

    Exercício exemplo

    Como exemplo, resolveremos um exercício simples de programação linear que ilustra os conceitos revistos.

    Enunciado:

    Uma startup busca fabricar máquinas de lavar falantes ao menor custo possível. Existem três maneiras de produzí-las: manualmente, semi automaticamente e automaticamente. A produção manual de uma máquina de lavar requer 1 minuto de trabalho qualificado, 40 minutos de trabalho não qualificado e 3 minutos de montagem. Para a produção semi automáticas os valores são, respectivamente, 4, 30 e 2 minutos e para a produção automática os valores são 8, 20 e 4 minutos. A startup dispões de 4500 minutos de trabalho qualificado, 36000 de trabalho não qualificado e 2700 minutos de montagem. Os custos de produção manual, semi automático e automático são, respectivamente, 70, 80 e 85 euros.

    1. Formule o problema como um PL e dê uma solução.

    No problema acima, as variáveis são o número \(x_1, x_2 e x_3\) de máquinas feitas por cada tipo de produção (manual, semi-automática e automática). O custo a ser minimizado é, portanto \(70x_1 + 80x_2 +85x_3\). Existem restrições quanto ao número de máquinas a serem produzidas e à capacidade de produção da fábrica (número de minutos para cada tipo de trabalho).

    Abaixo o problema completo:

    \[ \begin{cases} \textrm{minimize} f(x) = 70x_1 + 80x_2 +85x_3 \textrm{ s.a.}\\ x_1 + x_2 + x_3 = 999\\ x_1 + 4x_2 + 8x_3 \leq 4500\\ 40x_1 + 30x_2 + 20x_3 \leq 36000\\ 3x_1 + 2x_2 + 4x_3 \leq 2700\\ x \geq 0 \end{cases} \]

    Solução usando scipy:

    import numpy as np
    from scipy.optimize import linprog
    from numpy.linalg import solve
    
    A_eq = np.array([[1,1,1]])
    b_eq = np.array([999])
    
    A_ub = np.array([
    [1, 4, 8],
    [40,30,20],
    [3,2,4]])
    
    b_ub = np.array([4500, 36000,2700])
    
    c = np.array([70, 80, 85])
    
    res = linprog(c, A_eq=A_eq, b_eq=b_eq, A_ub=A_ub, b_ub=b_ub,
    bounds=(0, None))
    print('Valor otimo:', res.values()[3], '\nX:', res.values()[4])
    
    Valor otimo: 73725.0
    X: [ 636.  330.   33.]
    

    2. Converta o problema para sua forma canônica

    Para converter o problema acima para a forma canônica é necessário adicionar variáveis de folga para as restrições de desigualdade. Desta forma, o problema se torna:

    \[ \begin{cases} \textrm{minimize} f(x) = 70x_1 + 80x_2 +85x_3 \textrm{ s.a.}\\ x_1 + x_2 + x_3 = 999\\ x_1 + 4x_2 + 8x_3 + s_1 = 4500\\ 40x_1 + 30x_2 + 20x_3 + s_2 = 36000\\ 3x_1 + 2x_2 + 4x_3 + s_3 = 2700\\ x \geq 0, s \geq 0 \end{cases} \]

    Solução usando scipy:

    A = np.array([
    [1, 1, 1, 0, 0, 0],
    [1, 4, 8, 1, 0, 0],
    [40, 30, 20, 0, 1, 0],
    [3, 2, 4, 0, 0, 1]])
    
    b = np.array([999, 4500, 36000, 2700])
    c = np.array([70, 80, 85, 0, 0, 0])
    
    res = linprog(c, A_eq=A, b_eq=b, bounds=(0, None))
    print('Valor otimo:', res.values()[3], '\nX:', res.values()[4])
    
    Valor otimo: 73725.0
    X: [  636.   330.    33.  2280.     0.     0.]
    

    Note que \(s_1 = 2280\) implica que \(x_1 + 4x_2 + 8x_3 \leq 4500\).

    3. Apresente o dual do problema e mostre que ambas formulações obtem os mesmos resultados.

    O problema dual correspondente é dado por:

    \[ \begin{cases} \textrm{maximize} f(y) = 999y_1 + 4500y_2 + 360000y_3 + 2700y_4 \textrm{ s.a.} \\ y_1 + y_2 + 40y_3 + 3y_4 \leq 70 \\ y_1 + 4y_2 + 30y_3 + 2y_4 \leq 80 \\ y_1 + 8y_2 + 20y_3 + 4y_4 \leq 85 \\ y_2 \leq 0 \\ y_3 \leq 0 \\ y_4 \leq 0 \end{cases} \]

    Solução usando scipy: Para a solução usando \(scipy\) podemos utilizar as variáveis declaradas anteriormente. Note que como linprog resolve problemas de minimização, o sinal da função objetivo está trocado.

    res = linprog(-b, A_ub=A.T, b_ub=c, bounds=[(None,None), (None,None),
    (None,None), (None,None)])
    y = res.values()[4]
    print('Valor otimo:', -res.values()[3], '\nY:', y)
    
    u = c - A.T.dot(y)
    Ar = A[:, np.abs(u)< 1e-10]
    x_1 = solve(Ar, b)
    print(x_1)
    x = np.array([0.0] * len(c))
    x[np.abs(u)< 1e-10] = x_1
    x[np.abs(u)> 1e-10] = 0
    print('Solucao X a partir de Y:', x)
    
    Valor otimo: 73725.0
    Y: [ 108.33333333    0.           -0.83333333   -1.66666667]
    [  636.   330.    33.  2280.]
    Solucao X a partir de Y: [  636.   330.    33.  2280.     0.     0.]