49
TÓPICOS DE INVERSÃO EM GEOFÍSICA Vanderlei Coelho de Oliveira Junior Leonardo Uieda Versão 1.1 2011

Tópicos de inversão em geofísica - Cloud Object Storage | Store … · uma hipótese) e dados preditos é o que se conhece como problema direto. Nessesentido,dadoumconjuntodeparâmetros,épossívelcalcularosdadospredi-tos

Embed Size (px)

Citation preview

TÓPICOS DE

INVERSÃO EM GEOFÍSICA

Vanderlei Coelho de Oliveira Junior

Leonardo Uieda

Versão 1.1

2011

Tópicos inversão em geofísica

Copyright c© 2011-2014 Vanderlei Coelho de Olivera Junior & Leonardo Uieda

This work is licensed under theCreative Commons Attribution 4.0 International License.

To view a copy of this license, visithttp://creativecommons.org/licenses/by/4.0/

ii

Sumário

1 Introdução 1

2 Formulação matemática do problema inverso 8

2.1 Problemas lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2 Problemas não-lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Regularização 19

3.1 Norma mínima (Tikhonov de ordem 0) . . . . . . . . . . . . . . . . . . 23

3.2 Igualdade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.3 Suavidade (Tikhonov de ordem 1) . . . . . . . . . . . . . . . . . . . . . 26

3.4 Variação total . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

4 Procedimentos práticos 31

4.1 Ajustar os dados observados . . . . . . . . . . . . . . . . . . . . . . . . 31

4.2 Determinar um valor para o parâmetro de regularização . . . . . . . . . 32

5 Leitura recomendada 35

Referências Bibliográficas 36

A Operações com matrizes 37

A.1 Definições . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

A.2 Derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

A.3 Gradientes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

SUMÁRIO iii

A.4 Hessianas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

A.5 Expansão em série de Taylor . . . . . . . . . . . . . . . . . . . . . . . . 42

A.5.1 Até primeira ordem . . . . . . . . . . . . . . . . . . . . . . . . . 43

A.5.2 Até segunda ordem . . . . . . . . . . . . . . . . . . . . . . . . . 43

1

Capítulo 1

Introdução

Em ciências aplicadas, ao nos depararmos com um determinado fenômeno a ser estu-

dado, a pergunta imediata é: como estudar tal fenômeno? Em problemas geofísicos,

esse fenômeno é, em geral, um sistema físico. O procedimento científico utilizado para

estudá-lo começa com a realização de observações de uma determinada grandeza física

que, supostamente, seja causada pelo sistema físico em questão.

Exemplos de sistemas físicos encontrados na geofísica são:

• Propagação de ondas elásticas: sísmica (Figura 1)

• Teoria do potencial: gravimetria (Figura 2) e magnetometria

• Difusão de correntes elétricas: SEV (Figura 3)

Uma vez que não conhecemos o sistema físico, a única informação que temos são as

observações (ou dados observados). Embora o sistema físico seja desconhecido, alguma

informação adicional (em geral proporcionada pela geologia) nos fornece um conjunto

de hipóteses (Figuras 4, 5 e 6).

Isso nos leva a seguinte pergunta: como testar se as observações podem ser expli-

cadas por uma determinada hipótese? Para responder esta pergunta, precisamos ser

capazes de calcular os dados preditos pela hipótese em questão. Matematicamente,

1. Introdução 2

isso significa que precisamos encontrar uma função que relaciona a hipótese aos dados

preditos.

hipótesedadospreditos

Uma função é descrita em termos de parâmetros. Portanto, calcular os dados

preditos por uma determinada hipótese significa encontrar uma função f tal que:

dados preditos = f(hipótese).

Para tanto, é preciso descrever a hipótese em termos de parâmetros (Figuras 7, 8, 9),

isto é:

dados preditos = f(parâmetros).

O raciocínio desenvolvido acima nos permite definir um importante conceito:

Estabelecer a relação funcional entre parâmetros (que descrevem

uma hipótese) e dados preditos é o que se conhece como problema

direto.

Nesse sentido, dado um conjunto de parâmetros, é possível calcular os dados predi-

tos. Isso nos permite testar se uma hipótese é capaz de explicar os dados observados.

Esse teste é feito por meio da comparação entre os dados observados e os dados preditos

por uma determinada hipótese. Se os dados preditos reproduzem os dados observados,

pode-se dizer que a hipótese explica os dados observados e, portanto, é válida. Nesse

contexto, é possível definir os seguintes conceitos importantes:

Fornecer os parâmetros e comparar os respectivos dados preditos

com os dados observados é o que se conhece como modelagem di-

reta.

1. Introdução 3

e

Estimar os parâmetros automaticamente, de forma que os dados pre-

ditos sejam os mais próximos possíveis aos dados observados, é

conhecido como problema inverso.

ParâmetrosDados

preditos

Modelagem direta

Dadosobservados

calcula compara

ParâmetrosDados

observados

Problema inverso

estima

1. Introdução 4

x

z v2

v1

v3

v4

t

u

Sismograma

fonte receptores

raios

Figura 1: Exemplo de sistema físico que caracteriza uma aquisição sísmica. Nestecaso o sistema envolve ondas sísmicas (representadas pelos raios) que se propagam emmeios com velocidades V1, V2, V3 e V4. A grandeza física que podemos observar é odeslocamento dos sismômetros u (neste caso somente na vertical). Estas observaçõessão tipicamente dispostas em um sismograma.

x

z

gz

Δρ

Figura 2: Exemplo de sistema físico que caracteriza um levantamento gravimétrico.Neste caso o sistema envolve o efeito gravitacional causado por um contraste de densi-dade anômalo ∆ρ. A grandeza física que podemos observar é a anomalia da gravidadegz.

1. Introdução 5

x

z

ρa

AB/2

A BM N

ρ(z)

Figura 3: Exemplo de sistema físico que caracteriza uma sondagem elétrica vertical(SEV). Neste caso o sistema envolve a condução de correntes elétricas em um meiocom resistividades elétrica ρ que varia com a profundidade. A grandeza física quepodemos observar é diferência de potencial elétrico entre os eletrodos M e N. Estagrandeza é geralmente convertida para resistividade aparente ρa.

x

z

Figura 4: Exemplo de hipótese que descreve o sistema físico referente a uma aquisiçãosísmica. Modelo de duas camadas plano-paralelas, homogêneas e isotrópicas.

x

z Δρ

Figura 5: Exemplo de hipótese que descreve o sistema físico referente a uma levanta-mento gravimétrico. Modelo de um corpo elipsoidal homogêno.

1. Introdução 6

x

z

A BM N

Figura 6: Exemplo de hipótese que descreve o sistema físico referente a uma sondagemelétrica vertical. Modelo de três camadas plano-paralelas, homogêneas e isotrópicas.

x

z

v1

v2

h1

t

u

Sismograma

Figura 7: Exemplo de parametrização da hipótese na Figura 4. Neste caso, a hipóteseé descrita em termos dos parâmetros h1, V1 e V2. Esta hipótese produz o dado preditorepresentado pela linha pontilhada vermelha.

1. Introdução 7

x

z Δρ

gz

Figura 8: Exemplo de parametrização da hipótese na Figura 5. Neste caso, a hipótese édescrita em termos dos parâmetros ∆ρ1, ∆ρ2, . . . , ∆ρM que representam os constrastesde densidade de cada célula. Esta hipótese produz o dado predito representado pelalinha pontilhada vermelha.

x

z

A BM N

ρ1

ρ2

ρ3

h1

h2

ρa

AB/2

Figura 9: Exemplo de parametrização da hipótese na Figura 6. Neste caso, a hipóteseé descrita em termos dos parâmetros h1, ρ1, h2, ρ2 e ρ3. Esta hipótese produz o dadopredito representado pela linha pontilhada vermelha.

8

Capítulo 2

Formulação matemática do problema

inverso

Dado um conjunto de N observações, feitas em diferentes posições, tempos, etc., defi-

nimos um vetor de dados observados

d o =

d o1

d o2...

d oN

, (2.1)

em que d oi , i = 1, 2, 3, . . . , N , é o dado observado na i-ésima posição, tempo, etc. De

forma análoga, definimos um vetor de dados preditos

d p =

d p1

d p2...

d pN

, (2.2)

em que d pi , i = 1, 2, 3, . . . , N , é o dado predito calculado na mesma posição, tempo,

etc., que o i-ésimo dado observado. Continuando no espírito de definição de vetores,

definimos também um vetor que agrupa todos os M parâmetros, denominado vetor de

2. Formulação matemática do problema inverso 9

parâmetros

p =

p1

p2

...

pM

, (2.3)

em que pj, j = 1, 2, 3, . . . ,M , é o j-ésimo parâmetro.

Como vimos no Capítulo 1, os dados preditos são descritos por uma função dos

parâmetros, ou seja,

d pi = fi(p) . (2.4)

Desta forma, podemos dizer que o vetor de dados preditos é uma função dos parâmetros

d p = f(p) =

f1(p)

f2(p)

...

fN(p)

. (2.5)

O problema inverso consiste em encontrar um vetor de parâmetros p que produza

os dados preditos mais próximos possívies dos dados observados. Para determinar a

“proximidade” entre os dados observados e os dados preditos, é necessário quantificar

a distância entre eles. Isto é feito em termos da norma do vetor de resíduos

r = d o − f(p), (2.6)

em que d o é o vetor de dados observados e f(p) é o vetor de dados preditos.

Para quantificar a distância entre os dados observados e os dados preditos utiliza-se,

usualmente, o quadrado da norma quadrática (também conhecida como norma `2 ou

norma Euclidiana) do vetor de resíduos

2. Formulação matemática do problema inverso 10

‖r‖22 =

N∑i=1

[d oi − fi(p)]2 . (2.7)

Esta equação é também uma função escalar dos parâmetros. Assim sendo, definimos

uma função φ(p), chamada de função do ajuste, como

φ(p) = ‖r‖22 =

N∑i=1

[d oi − fi(p)]2 . (2.8)

Lembramos que o quadrado da norma Euclidiana de um vetor é igual ao produto

escalar do vetor com ele mesmo, ou seja,

‖r‖22 = r · r = r1r1 + r2r2 + · · ·+ rNrN . (2.9)

Como o vetor r é um vetor coluna (matriz com uma coluna), podemos escrever o

produto escalar da equação 2.9 como

r · r = rT r =

[r1 r2 . . . rN

]

r1

r2

...

rN

. (2.10)

Dessa forma podemos reescrever a função do ajuste (equação 2.8) como

φ(p) = rT r =[d o − f(p)

]T [d o − f(p)

]. (2.11)

Neste ponto é importante ressaltar o seguinte conceito:

A função do ajuste é uma função escalar que quantifica a distân-

cia entre os dados observados e os dados preditos para um deter-

minado vetor de parâmetros p.

Perante este conceito, o problema inverso consiste em determinar um vetor p que

minimiza a função do ajuste φ(p). Matematicamente, isso equivale a encontrar o vetor

2. Formulação matemática do problema inverso 11

p tal que o gradiente da função φ(p) avaliado em p seja igual ao vetor nulo. Este vetor

p é um ponto extremo da função φ(p).

O gradiente da função φ(p), avaliado em um p qualquer, é um vetorM -dimensional

definido como (ver Apêndice A)

∇φ(p) =

∂φ(p)

∂p1

∂φ(p)

∂p2...

∂φ(p)

∂pM

, (2.12)

em que ∇ é o operador gradiente (Apêndice A). A partir da equação 2.11, a expressão

para o i-ésimo elemento do vetor gradiente é

∂φ(p)

∂pi=

∂pi

[d o − f(p)

]T [d o − f(p)

]

=

−∂f(p)

∂pi

T [d o − f(p)

]︸ ︷︷ ︸

escalar

+

−[d o − f(p)

]T ∂f(p)

∂pi

︸ ︷︷ ︸

escalar

. (2.13)

Lembrando que o transposto de um escalar é igual a ele mesmo, podemos tomar o

transposto do segundo termo do lado direito da equação 2.13, obtendo

∂φ(p)

∂pi= −2

∂f(p)

∂pi

T [d o − f(p)

], (2.14)

em que∂f(p)

∂pié um vetor N -dimensional (ver Apêndice A) dado por

2. Formulação matemática do problema inverso 12

∂f(p)

∂pi=

∂f1(p)

∂pi

∂f2(p)

∂pi...

∂fN(p)

∂pi

. (2.15)

Substituindo as equações 2.14 e 2.15 na expressão do gradiente (equação 2.12)

obtemos

∇φ(p) =

−2∂f(p)

∂p1

T [d o − f(p)

]−2

∂f(p)

∂p2

T [d o − f(p)

]...

−2∂f(p)

∂pM

T [d o − f(p)

]

= −2

∂f(p)

∂p1

T

∂f(p)

∂p2

T

...

∂f(p)

∂pM

T

[d o − f(p)

].

(2.16)

Por fim, o gradiente da função do ajuste φ(p) pode ser escrito como

∇φ(p) = −2 ¯G(p)T[d o − f(p)

]. (2.17)

em que

2. Formulação matemática do problema inverso 13

¯G(p) =

[∂f(p)

∂p1

∂f(p)

∂p2

. . .∂f(p)

∂pM

]=

∂f1(p)

∂p1

∂f1(p)

∂p2

. . .∂f1(p)

∂pM

∂f2(p)

∂p1

∂f2(p)

∂p2

. . .∂f2(p)

∂pM...

... . . . ...∂fN(p)

∂p1

∂fN(p)

∂p2

. . .∂fN(p)

∂pM

. (2.18)

A matriz ¯G(p) de dimensão N ×M é a matriz Jacobiana de f(p).

Em problemas inversos, essa matriz ¯G é comumente denominada ma-

triz de sensibilidade, uma vez que o i-ésimo elemento de sua j-ésima

coluna expressa a sensibilidade do i-ésimo dado predito em rela-

ção à variações do j-ésimo parâmetro.

As equações 2.17 e 2.18 mostram que o gradiente da função do ajuste φ(p) depende

do vetor de dados preditos f(p) (equação 2.5) e sua derivada em relação aos parâmetros

p. Sendo assim, a função f que relaciona os parâmetros aos dados preditos determina

o comportamento do gradiente de φ(p). Nas próximas seções analisaremos os casos

em que a função f é linear ou não-linear em relação aos parâmetros. Essa análise nos

permitirá compreender a influência da função f na busca pelo vetor de parâmetros p

que minimiza a função φ(p).

2.1 Problemas lineares

Se a função fi(p) (equação 2.4), que relaciona o vetor de parâmetros p ao i-ésimo dado

predito, for uma combinação linear dos parâmetros, ela terá a seguinte forma

fi(p) = gi1p1 + gi2p2 + · · ·+ giMpM + bi , (2.19)

em que gij, j = 1, 2, · · · ,M , e bi são constantes. Portanto, a derivada de fi(p) em

2. Formulação matemática do problema inverso 14

relação ao j-ésimo parâmetro pj é dada por

∂fi(p)

∂pj= gij , (2.20)

que não depende dos parâmetros. Para um conjunto de N dados, o vetor de dados

preditos tem a seguinte forma

f(p) = ¯Gp+ b , (2.21)

sendo ¯G a matriz de sensibilidade (equação 2.18) e b um vetor de constantes. Substi-

tuindo a expressão acima no gradiente da função do ajuste (equação 2.17), chegamos

a

∇φ(p) = −2 ¯GT(d o − ¯Gp− b

). (2.22)

Seja p o vetor que minimiza a função φ(p), o gradiente desta função avaliado em p

é igual ao vetor nulo. Isto nos leva ao sistema de equações

¯GT ¯Gp = ¯GT(d o − b

). (2.23)

Este sistema é conhecido como sistema de equações normais. A solução para este

sistema é

p = ( ¯GT ¯G)−1 ¯GT(d o − b

), (2.24)

que é conhecido como estimador de mínimos quadrados.

As equações 2.23 e 2.24 mostram que o vetor p que minimiza a função φ(p) pode

ser obtido diretamente a partir da matriz ¯G e dos vetores b e d o. Isto caracteriza um

problema inverso linear.

2. Formulação matemática do problema inverso 15

2.2 Problemas não-lineares

Nesta seção analisaremos o caso em que a função fi(p) (equação 2.4) não é uma combi-

nação linear dos parâmetros. Neste caso, a derivada de fi(p) em relação aos parâmetros

também será uma função dos parâmetros. Logo, dependendo das características da fun-

ção fi(p), o gradiente da função do ajuste φ(p) pode ser nulo para mais de um valor

de p. Em outras palavras, a função φ(p) pode possuir vários pontos extremos além

daquele em que esta função seja mínima. Por esta razão, o cálculo do vetor p que

minimiza a função do ajuste deve ser feito de forma iterativa. Isto difere do problema

inverso linear e caracteriza um problema inverso não-linear.

A busca pelo vetor p que minimiza uma determinada função faz parte de uma área

da matemática conhecida como otimização, dentro da qual existem diversos métodos

(Kelley, 1999). O procedimento padrão para realizar esta busca de forma iterativa é

começar com uma determinada aproximação inicial p0 e calcular uma correção ∆p.

Esta correção é então aplicada à aproximação inicial dando origem a um novo vetor p1.

Este novo vetor serve como aproximação inicial para o cálculo de um segundo vetor p2,

e assim sucessivamente. O processo termina quando é encontrado um vetor p que seja

próximo ao vetor p que minimiza a função em questão. Em geral, não se tem garantia

de que o vetor p seja igual ao vetor p.

Dentre os vários métodos existentes, apresentaremos a seguir aquele conhecido como

método de Gauss-Newton. O procedimento para calcular a correção ∆p começa com a

expansão em série de Taylor da função a ser minimizada, neste caso φ(p). A expansão

é feita até segunda ordem e em torno da aproximação inicial p0 (ver Apêndice A)

p = p0 + ∆p , (2.25)

φ(p0 + ∆p) ≈ φ(p0) + ∇φ(p0)T∆p+1

2∆pT ¯∇φ(p0)∆p = ψ(p) , (2.26)

sendo ∇φ(p0) o vetor gradiente e ¯∇φ(p0) a matriz Hessiana da função φ(p), calculados

2. Formulação matemática do problema inverso 16

em p0.

A função ψ(p) (equação 2.26) é uma aproximação de segunda ordem para a função

φ(p) em torno do ponto p0. Tal como no problema inverso linear, desejamos encontrar

um ponto p que minimiza a função ψ(p), ou seja, onde seu gradiente é nulo. O gradiente

de ψ(p) é dado por (ver Apêndice A)

∇ψ(p) = ∇[φ(p0) + ∇φ(p0)T∆p+

1

2∆pT ¯∇φ(p0)∆p

]= ∇φ(p0) + ¯∇φ(p0)∆p . (2.27)

Logo, o incremento ∆p = p − p0, que leva da aproximação inicial ao ponto p onde o

gradiente de ψ(p) é nulo, é a solução do sistema de equações

¯∇φ(p0)∆p = −∇φ(p0) . (2.28)

O cálculo do vetor gradiente de φ(p) foi feito anteriormente (equação 2.17), restando

então o cálculo da matriz Hessiana ¯∇φ(p). Para tanto, deriva-se o i-ésimo elemento

do vetor gradiente (equações 2.12 e 2.14) em relação ao j-ésimo elemento do vetor de

parâmetros pj (equação 2.3)

∂pj

(∂φ(p)

∂pi

)=

∂pj

(−2

∂f(p)

∂pi

T [d o − f(p)

])

=

(−2

∂2f(p)

∂pj∂pi

T [d o − f(p)

])︸ ︷︷ ︸

escalar

+

(2∂f(p)

∂pi

T∂f(p)

∂pj

)︸ ︷︷ ︸

escalar

.

(2.29)

Esta equação 2.29 representa o j-ésimo elemento da i-ésima linha da matriz Hessi-

ana (ver Apêndice A). No método de Gauss-Newton, o termo envolvendo as segundas

derivadas de f(p) é negligenciado. Desta forma a equação 2.29 pode ser aproximada

por

2. Formulação matemática do problema inverso 17

∂pj

(∂φ(p)

∂pi

)≈ 2

∂f(p)

∂pi

T∂f(p)

∂pj. (2.30)

Sendo assim, matriz Hessiana avaliada em p0 é dada por

¯∇φ(p0) ≈ 2

∂f(p0)

∂p1

T∂f(p0)

∂p1

∂f(p0)

∂p1

T∂f(p0)

∂p2

. . .∂f(p0)

∂p1

T∂f(p0)

∂pM

∂f(p0)

∂p2

T∂f(p0)

∂p1

∂f(p0)

∂p2

T∂f(p0)

∂p2

. . .∂f(p0)

∂p2

T∂f(p0)

∂pM...

... . . . ...

∂f(p0

∂pM

T∂f(p0)

∂p1

∂f(p0)

∂pM

T∂f(p0)

∂p2

. . .∂f(p0)

∂pM

T∂f(p0)

∂pM

,

ou

¯∇φ(p0) ≈ 2

∂f(p0)

∂p1

T

∂f(p0)

∂p2

T

...

∂f(p0)

∂pM

T

[∂f(p0)

∂p1

∂f(p0)

∂p2

. . .∂f(p0)

∂pM

]= 2 ¯G(p0)T ¯G(p0) , (2.31)

em que ¯G(p0) é a matriz de sensibilidade (equação 2.18) avaliada em p0.

A partir das equações 2.17 e 2.31, o sistema de equações 2.28 pode ser escrito como

¯G(p0)T ¯G(p0)∆p = ¯G(p0)T[d o − f(p0)

], (2.32)

em que d o é o vetor de dados observados e f(p0) é o vetor de dados preditos avaliado

em p0. Este sistema de equações é o sistema de equações normais do problema inverso

não-linear.

A equação 2.32 descreve o cálculo da correção ∆p em uma determinada iteração

do método Gauss-Newton. Para o caso em que a função fi(p) é uma combinação

2. Formulação matemática do problema inverso 18

linear dos parâmetros (equação 2.19), a equação 2.32 se reduz ao sistema de equações

normais do problema inverso linear (equação 2.23) 1. Isto mostra que a solução do

problema inverso linear é um caso particular da solução do problema inverso não-

linear pelo método Gauss-Newton. Outra observação importante é a semelhança entre

as equações 2.23 e 2.32, o que evidencia que um problema inverso não-linear é uma

sucessão de problemas inversos lineares.

1Dica: Lembre que nos casos lineares f(p) = ¯Gp+ b.

19

Capítulo 3

Regularização

As equações 2.23 do problema inverso linear e 2.32 do problema inverso não-linear são

equações normais. Estas equações são sistemas lineares cuja matriz é quadrada. Para

que a solução de uma equação normal seja única, é necessário que esta matriz tenha

posto completo. Uma condição para que uma matriz tenha posto completo é que todas

as suas colunas (ou linhas) sejam linearmente independentes. Esta condição é equiva-

lente a dizer que a matriz possui determinante diferente de zero.

Em problemas geofísicos, é comum que a matriz do sistema normal possua determi-

nante próximo de zero. Ou seja, para fins práticos pode-se considerar que a matriz não

tenha posto completo. Isto faz com que o problema inverso geofísico seja um problema

mal posto.

Um problema mal-posto apresenta, principalmente, instabilidade e falta

de unicidade da solução.

A instabilidade é a alta variabilidade dos parâmetros perante peque-

nas variações nos dados.

Exemplo 3.1. Seja um conjunto de dados preditos A e outro conjunto de dados preditos

ligeiramente diferentes B. Se o problema inverso apresenta instabilidade, os parâme-

tros que produzem os dados preditos A são consideravelmente diferentes daqueles que

produzem os dados preditos B.

3. Regularização 20

A falta de unicidade é a existência de vários conjuntos de parâme-

tros que produzem os mesmos dados preditos.

Se o problema inverso apresenta falta de unicidade, os dados observados podem ser

explicados por vários conjuntos de parâmetros diferentes. Em termos práticos, existem

vários vetores de parâmetros diferentes que minimizam a função do ajuste (equação

2.11).

Os problemas inversos encontrados na geofísica são, em geral, mal-postos. Os princi-

pais motivos para isto são: a presença de ruído nos dados observados; e, principalmente,

a própria natureza do problema inverso. Neste último caso, mesmo se fôssemos capazes

de obter dados completamente isentos de ruído, ainda teríamos diversas combinações

de parâmetros capazes de explicar os dados observados (i.e., falta de unicidade). Por

este motivo, quando nos deparamos com problemas inversos na geofísica, é essencial a

utilização de regularização.

A regularização é um procedimento matemático que contorna os problemas de ins-

tabilidade e falta de unicidade em problemas inversos mal-postos. Esse procedimento

equivale a impor restrições aos parâmetros a serem estimados. Desta forma, em um pro-

blema inverso regularizado buscamos estimar um conjunto de parâmetros que ajustam

os dados observados e satisfaçam determinadas restrições. Estas restrições introdu-

zem informações a priori no problema inverso. As informações podem ser de natureza

geológica ou meramente matemática. Em geral, a introdução de informações a priori

é feita por meio de funções regularizadoras, que são funções escalares que dependem

dos parâmetros. Para incorporar informação a priori ao ajuste dos dados, formamos a

função objetivo

Ω(p) = φ(p) + µθ(p) , (3.1)

em que φ(p) é a função do ajuste (equação 2.11), θ(p) é uma função regularizadora e µ

é um escalar positivo denominado parâmetro de regularização. Desta forma, o problema

3. Regularização 21

inverso regularizado é definido como estimar um vetor de parâmetros p que minimiza

a função objetivo.

Neste ponto é importante ressaltar que

O parâmetro de regularização µ controla a importância relativa

entre o ajuste aos dados observados e a concordância com a infor-

mação a priori.

O valor de µ é definido pelo usuário da inversão e por isso é importante compreender

sua influência no resultado obtido (parâmetros estimados). Valores altos de µ tornam

o problema inverso bem posto e fazem com que os parâmetros estimados satisfaçam

quase completamente as informações a priori. Porém, isto geralmente faz com que haja

um desajuste entre os dados observados e preditos. Por outro lado, valores baixos de

µ fazem com que os parâmetros estimados ajustem os dados observados. No entanto,

a estimativa poderá ser não-única e/ou instável, dependendo de quanto o problema

inverso for mal-posto. Idealmente, deve-se encontrar um valor de µ que proporcione

um bom ajuste e satisfaça as informações a priori o suficiente para estabilizar a solução.

Procedimentos práticos para determinação do valor de µ serão discutidos no Capítulo

4.

No problema inverso regularizado, a linearidade não depende apenas da função fi(p)

(equação 2.4) que relaciona os parâmetros aos dados preditos, mas também da função

regularizadora. Nesse caso, para que o problema inverso seja linear, é necessário que

tanto fi(p) como o gradiente da função regularizadora sejam combinações lineares dos

parâmetros. Se ao menos uma destas funções for não-linear, o problema inverso deverá

ser resolvido como um problema inverso não-linear.

Como foi observado anteriormente (Seção 2.2), o problema inverso linear é um caso

particular do problema inverso não-linear. Por esta razão, a formulação geral para o

problema inverso regularizado será feita seguindo os procedimentos adotados para o

problema inverso não-linear. Assim sendo, iniciaremos expandindo a função objetivo

(equação 3.1) em série de Taylor até segunda ordem

3. Regularização 22

Ω(p0 + ∆p) ≈ Ω(p0) + ∇Ω(p0)T∆p+1

2∆pT ¯∇Ω(p0)∆p , (3.2)

sendo o vetor gradiente ∇Ω(p0) e a matriz Hessiana ¯∇Ω(p0) de Ω(p) dados por

∇Ω(p0) = ∇φ(p0) + µ∇θ(p0) (3.3)

e

¯∇Ω(p0) = ¯∇φ(p0) + µ ¯∇θ(p0) , (3.4)

em que µ é o parâmetro de regularização, ∇φ(p0) e ¯∇φ(p0) são o gradiente e a Hessiana

da função do ajuste (equações 2.17 e 2.31) e ∇θ(p0) e ¯∇θ(p0) são o gradiente e a

Hessiana da função regularizadora, todos avaliados em p0.

Seguindo a dedução feita para o problema inverso não-linear (Seção 2.2), a correção

∆p em uma determinada iteração do método Gauss-Newton é

¯∇Ω(p0)∆p = −∇Ω(p0) . (3.5)

Substituindo o gradiente e a Hessiana da função do ajuste (equações 2.17 e 2.31)

obtemos

[2 ¯G(p0)T ¯G(p0) + µ ¯∇θ(p0)

]∆p = 2 ¯G(p0)T

[d o − f(p0)

]− µ∇θ(p0) . (3.6)

Esta é a equação para o problema inverso regularizado. 1 A seguir, apresentaremos

diferentes funções regularizadoras comumente encontradas em problemas inversos na

geofísica. Também mostraremos como fica a equação 3.6 para cada um dos casos.1Tente obter a equação normal para o caso em que f(p) é linear.

3. Regularização 23

3.1 Norma mínima (Tikhonov de ordem 0)

A função regularizadora mais comumente usada é a chamada norma mínima (também

conhecida como ridge regression ou Tikhonov de ordem zero). Como seu nome sugere,

esta função é utilizada para incorporar a informação de que o vetor de parâmetros

deve ter a norma quadrática (`2 ou Euclidiana) mínima. Isto é, os parâmetros devem

ser assumir valores mais próximos possíveis a zero. De forma similar a equação 2.9, a

função regularizadora de norma mínima tem a seguinte forma

θNM(p) = pT p . (3.7)

O vetor gradiente ∇θNM(p) e a matriz Hessiana ¯∇θNM(p) (ver Apêndice A) desta

função são, respectivamente,

∇θNM(p) = 2 ¯Ip (3.8)

e

¯∇θNM(p) = 2 ¯I , (3.9)

em que ¯I é a matriz identidade de dimensão M ×M , lembrando que M é o número de

parâmetros. Note que o gradiente da função regularizadora de norma mínima é uma

combinação linear dos parâmetros.

Para o caso em que a função fi(p) que relaciona os dados preditos aos parâme-

tros também é linear (equação 2.19), a equação normal do problema inverso linear

regularizado, para o caso da regularização de norma mínima, é

(¯GT ¯G+ µ ¯I

)p = ¯GT

(d o − b

), (3.10)

em que µ é o parâmetro de regularização, d o é o vetor de dados observados, ¯G é a

matriz de sensibilidade, b é um vetor de constantes (equação 2.21) e p é a solução de

3. Regularização 24

norma mínima para o problema inverso linear. 2

Já para o caso em que fi(p) é não-linear, o problema inverso torna-se também não-

linear. Assim sendo, a equação normal do problema inverso não-linear regularizado,

para o caso da regularização de norma mínima, é

[¯G(p0)T ¯G(p0) + µ ¯I

]∆p = ¯G(p0)T

[d o − f(p0)

]− µp0 . (3.11)

em que f(p0) é o vetor de dados preditos avaliado em p0 e ∆p é a correção a ser aplicada

a p0.

3.2 Igualdade

Há casos em que se conhece um valor aproximado de um ou mais parâmetros. Esta

informação pode ser proveniente de furos de sondagem, afloramentos, etc. Nestes casos,

é desejável que o valor estimado para estes parâmetros seja o mais próximo possível dos

valores conhecidos (e preestabelecidos). Para tanto, utilizamos a função regularizadora

de igualdade

θIG(p) = (p− p a)T ¯AT ¯A (p− p a) , (3.12)

em que p a é um vetor cujo j-ésimo elemento é: o valor conhecido (preestabelecido) do

j-ésimo parâmetro; ou zero 3 caso não haja um valor conhecido do j-ésimo parâmetro.

A matriz ¯A é uma matriz diagonal quadrada de dimensão M ×M , lembrando que M é

o número de parâmetros. O j-ésimo elemento da diagonal de ¯A é 1 (um) caso haja um

valor conhecido (preestabelecido) para j-ésimo parâmetro, ou 0 (zero) caso contrário.

Exemplo 3.2. Digamos que em um determinado problema inverso a função fi(p)

depende de três parâmetros. Além disso, desejamos que o segundo parâmetro p2 seja2 Onde foi parar o termo ∇θNM (p0)? Dica: Mostre que, se o problema inverso é linear, a estimativa

p não depende da aproximação inicial p0.3Na prática, pode-se utilizar qualquer valor.

3. Regularização 25

o mais próximo possível do valor 26. Podemos então usar a função regularizadora de

igualdade para impor essa restrição. Neste caso, os vetores p e p a e a matriz ¯A serão

p =

p1

p2

p3

, p a =

0

26

0

e ¯A =

0 0 0

0 1 0

0 0 0

.

O vetor gradiente ∇θIG(p) e a matriz Hessiana ¯∇θIG(p) (ver Apêndice A) da função

regularizadora de igualdade são, respectivamente,

∇θIG(p) = 2 ¯AT ¯A (p− p a) (3.13)

e

¯∇θIG(p) = 2 ¯AT ¯A . (3.14)

Para o caso em que a função fi(p) que relaciona os dados preditos aos parâme-

tros também é linear (equação 2.19), a equação normal do problema inverso linear

regularizado, para o caso da regularização de igualdade, é

(¯GT ¯G+ µ ¯AT ¯A

)p = ¯GT

(d o − b

)+ µ ¯AT ¯Ap a, (3.15)

em que µ é o parâmetro de regularização, d o é o vetor de dados observados, ¯G é a

matriz de sensibilidade, b é um vetor de constantes (equação 2.21) e p é a solução do

problema inverso linear com regularização de igualdade.

Já para o caso em que fi(p) é não-linear, o problema inverso torna-se também não-

linear. Assim sendo, a equação normal do problema inverso não-linear regularizado,

para o caso da regularização de igualdade, é

[¯G(p0)T ¯G(p0) + µ ¯AT ¯A

]∆p = ¯G(p0)T

[d o − f(p0)

]− µ ¯AT ¯A [p0 − p a] . (3.16)

3. Regularização 26

em que f(p0) é o vetor de dados preditos avaliado em p0 e ∆p é a correção a ser aplicada

a p0.

3.3 Suavidade (Tikhonov de ordem 1)

Em certas situações é desejável que a distribuição dos parâmetros seja “suave”. Exis-

tem diversas interpretações para esta restrição, porém a mais comum é que parâmetros

espacialmente adjacentes devem ter valores mais próximos possível. Em outras pala-

vras, não devem haver variações abruptas entre parâmetros espacialmente adjacentes,

ou que a diferença entre estes parâmetros deve ser mínima.

A função regularizadora utilizada para incorporar a informação de suavidade tem

a seguinte forma

θSV (p) = vT v , (3.17)

em que v é um vetor com as diferêncas entre os parâmetros espacialmente adjacentes. O

vetor v é uma aproximação de diferenças finitas para a derivada espacial dos parâmetros

e pode ser escrito como

v = ¯R p , (3.18)

em que ¯R é uma matriz de diferenças finitas.

Exemplo 3.3. Seja o problema inverso de estimar o relevo de uma bacia sedimen-

tar. Uma possível parametrização seria discretizar a bacia em 7 prismas retangulares

justapostos com larguras fixas. Os parâmetros a serem estimados seriam então as 7

espessuras dos prismas (Figura 10). Impor suavidade ao relevo da bacia é equivalente

a minimizar as diferenças p1 − p2, p2 − p3, p3 − p4, p4 − p5, p5 − p6 e p6 − p7. Neste

caso, o vetor v das diferenças entre os parâmetros espacialmente adjacentes é dado por

3. Regularização 27

x

z

p1 p7p2

p3 p4

p5p6

Embasamento

Figura 10: Exemplo de parametrização do relevo de uma bacia sedimentar. Nestecaso, os 7 parâmetros utilizados são as espessuras de prismas retangulares justapostos.Impor suavidade ao relevo da bacia é equivalente a minimizar as diferenças p1 − p2,p2 − p3, p3 − p4, p4 − p5, p5 − p6 e p6 − p7.

v =

p1 − p2

p2 − p3

p3 − p4

p4 − p5

p5 − p6

p6 − p7

=

1 −1 0 0 0 0 0

0 1 −1 0 0 0 0

0 0 1 −1 0 0 0

0 0 0 1 −1 0 0

0 0 0 0 1 −1 0

0 0 0 0 0 1 −1

︸ ︷︷ ︸

¯R

p1

p2

p3

p4

p5

p6

p7

︸ ︷︷ ︸

p

. (3.19)

A função regularizadora de suavidade (também conhecida como Tikhonov de ordem

um) pode ser escrita como

θSV (p) = pT ¯RT ¯R p . (3.20)

O vetor gradiente ∇θSV (p) e a matriz Hessiana ¯∇θSV (p) (ver Apêndice A) desta função

são, respectivamente,

∇θSV (p) = 2 ¯RT ¯R p (3.21)

3. Regularização 28

e

¯∇θSV (p) = 2 ¯RT ¯R . (3.22)

O gradiente da função regularizadora de norma mínima é uma combinação linear

dos parâmetros. Para o caso em que a função fi(p) que relaciona os dados preditos aos

parâmetros também é linear (equação 2.19), a equação normal do problema inverso

linear regularizado, para o caso da regularização de suavidade, é

(¯GT ¯G+ µ ¯RT ¯R

)p = ¯GT

(d o − b

), (3.23)

em que µ é o parâmetro de regularização, d o é o vetor de dados observados, ¯G é a

matriz de sensibilidade, b é um vetor de constantes (equação 2.21) e p é a solução

suave para o problema inverso linear.

Já para o caso em que fi(p) é não-linear, o problema inverso torna-se também não-

linear. Assim sendo, a equação normal do problema inverso não-linear regularizado,

para o caso da regularização de suavidade, é

[¯G(p0)T ¯G(p0) + µ ¯RT ¯R

]∆p = ¯G(p0)T

[d o − f(p0)

]− µ ¯RT ¯R p0 . (3.24)

em que f(p0) é o vetor de dados preditos avaliado em p0 e ∆p é a correção a ser aplicada

a p0.

3.4 Variação total

Ao contrário do que vimos na Seção 3.3, há situações em que é desejável que hajam

algumas descontinuidades entre parâmetros espacialmente adjacentes. Nestes casos,

podemos utilizar a função regularizadora de variação total

θV T (p) =L∑k=1

|vk| , (3.25)

3. Regularização 29

em que vk, k = 1, 2, . . . , L, é o k-ésimo elemento do vetor de diferenças v (equação 3.18).

Como esta função não é diferenciável para valores de vk = 0, podemos aproximá-la por

θV T (p) ≈ θV Tβ (p) =L∑k=1

√v2k + β , (3.26)

sendo β um escalar positivo pequeno.

O vetor gradiente ∇θV Tβ (p) e a matriz Hessiana ¯∇θV Tβ (p) desta função são, respec-

tivamente, (Martins et al., 2011)

∇θV Tβ (p) = ¯RT q(p) (3.27)

e

¯∇θV Tβ (p) = ¯RT ¯Q(p) ¯R , (3.28)

sendo ¯R uma matriz de diferenças finitas (equação 3.18). O vetor q(p) e a matriz ¯Q(p)

são, respectivamente,

q(p) =

v1√v2

1 + β

v2√v2

2 + β...vL√v2L + β

(3.29)

e

¯Q(p) =

β

(v21 + β)

32

0 . . . 0

(v22 + β)

32

. . . 0

...... . . . ...

0 0 . . .β

(v2L + β)

32

. (3.30)

3. Regularização 30

Como o gradiente e a Hessiana da função regularizadora de variação total (equações

3.27 e 3.28) não são combinações lineares dos parâmetros, a utilização desta função

regularizadora torna o problema inverso não-linear. Assim sendo, a equação normal

do problema inverso não-linear regularizado, para o caso da regularização de variação

total, é

[¯G(p0)T ¯G(p0) + µ ¯RT ¯Q(p0) ¯R

]∆p = ¯G(p0)T

[d o − f(p0)

]− µ ¯RT q(p0) . (3.31)

em que f(p0) é o vetor de dados preditos avaliado em p0 e ∆p é a correção a ser aplicada

a p0.

31

Capítulo 4

Procedimentos práticos

Durante a solução de problemas inversos nos deparamos com diversos obstáculos, como:

• Por que a inversão não foi capaz de ajustar os dados observados?

• Como determinar um valor adequado para o parâmetro de regularização?

• Como analisar a estabilidade da solução?

Neste capítulo descreveremos alguns procedimentos práticos para abordar estas ques-

tões.

4.1 Ajustar os dados observados

Em situações reais, nem sempre é possível ajustar perfeitamente os dados observados.

Isso pode acontecer por três principais razões:

1. Os dados observados contêm ruído;

2. A função fi(p) escolhida para descrever a relação entre os parâmetros e os dados

preditos não é capaz de descrever o sistema físico de forma satisfatória;

3. Excesso de regularização;

4. Procedimentos práticos 32

O primeiro caso é esperado pois observações sempre contêm algum tipo de ruído.

Embora existam técnicas para a atenuação de ruído, é impossível que ele seja removido

completamente. Além disso, ruídos de caráter aleatório não são descritos pela relação

funcional entre os parâmetros e os dados preditos (equação 2.4). Sendo assim, não

ajustar perfeitamente os dados observados é aceitável, desde que o desajuste esteja

dentro do nível de ruído estimado dos dados.

Já no segundo caso, a função escolhida para descrever a relação entre os parâme-

tros e os dados preditos não representa o sistema físico de forma adequada. Nesse caso,

mesmo se os dados observados fossem isentos de ruído (situação impossível em condi-

ções reais), não seria possível estimar parâmetros que ajustassem os dados observados.

Isso é um indício de que a hipótese escolhida para representar o sistema físico não é vá-

lida. É crucial ressaltar que este é um resultado tão importante quanto encontrar uma

hipótese que descreva o sistema físico de forma satisfatória. Este tipo de resultado nos

permite descartar hipóteses (i.e., cenários geológicos) com base nos dados geofísicos.

O terceiro caso, em que há excesso de regularização, é devido a utilização de valores

demasiado elevados do parâmetro de regularização µ (equação 3.1). Valores elevados

de µ dão excessiva importância para a informação a priori (i.e., regularização), o que

causa um desajuste dos dados. Já valores baixos de µ priorizam o ajuste dos dados e

negligenciam a informação a priori. Neste caso, embora haja um bom ajuste dos dados,

o problema inverso se torna instável. Na Seção 4.2 trataremos da escolha de uma valor

adequado para µ.

4.2 Determinar um valor para o parâmetro de regu-

larização

O valor do parâmetro de regularização µ (equação 3.1) não pode ser determinado

para um caso geral. O valor de µ pode ser diferente para cada dado observado e

cada parametrização utilizada. Além disto, no caso em que haja mais de um função

4. Procedimentos práticos 33

regularizadora, pode ser desejável impor alguma delas mais fortemente que as outras.

Por exemplo, se utilizarmos ambas as funções regularizadoras de suavidade e igualdade

(equações 3.20 e 3.12, respectivamente), podemos desejar que a igualdade seja imposta

mais que a suavidade. No entanto, há procedimentos que permitem encontrar valores

de µ que sejam considerados adequados.

Um valor adequado de µ deve proporcionar um equilíbrio entre a estabilidade da

solução e o ajuste dos dados observados. Embora existam alguns métodos para escolha

automática do parâmetro de regularização (Aster et al., 2005), ainda não há uma

maneira de determinar um valor ótimo. Além disso, estes métodos automáticos são

aplicáveis somente a problemas inversos que utilizam uma única função regularizadora.

Sendo assim, apresentamos a seguir um procedimento prático que consiste em escolher

o menor valor possível para o parâmetro de regularização para que o problema inverso

seja bem-posto (i.e., com solução única e estável).

Seja um conjunto de N dados observados e M parâmetros, o procedimento prático

segue as seguintes etapas:

1. Gerar um vetor e de N valores aleatórios. Cada valor deve ser a realização de

uma variável aleatória com distribuição Gaussiana de média zero e desvio padrão

igual ao nível de ruído estabelecido para os dados;

2. Gerar um vetor d ′ de dados observados perturbados a partir da soma do vetor de

valores aleatórios e e o vetor de dados observados d 0

d ′ = d 0 + e ;

3. Repetir as etapas 1 e 2 para gerar um conjunto de Q vetores de dados observados

perturbados diferentes;

4. Estabelecer um valor pequeno para o parâmetro de regularização µ;

5. Utilizar este valor de µ em Q inversões para estimar Q vetores de parâmetros

4. Procedimentos práticos 34

diferentes. Cada vetor de parâmetros corresponde a um dos vetores de dados

observados perturbados;

6. Calcular o desvio padrão dos Q valores obtidos para o j-ésimo parâmetro;

7. Repetir a etapa 6 para cada um dos M parâmetros;

8. Se ao menos um dos M parâmetros apresentar um desvio padrão maior que um

valor máximo preestabelecido, considere que o valor de µ não é adequado. Nesse

caso, aumente ligeiramente o valor de µ e repita as etapas 5-7;

9. Se nenhum dos M parâmetros apresentar um desvio padrão maior que um valor

máximo preestabelecido, considere que o valor de µ é adequado e o problema

inverso está bem regularizado;

As etapas 1-3 e 5-8 acima são o procedimento adotado para a análise da estabilidade

da solução.

35

Capítulo 5

Leitura recomendada

Livros sobre problemas inversos:

• Tarantola (2005)

• Menke (1989)

• Aster et al. (2005)

Métodos potenciais:

• Silva et al. (2001)

• Medeiros and Silva (1996)

• Martins et al. (2011)

• Barbosa and Silva (2011)

Otimização:

• Kelley (1999)

36

Referências Bibliográficas

Aster, R. C., B. Borchers, and C. H. Thurber, 2005, Parameter estimation and inverseproblems: Elsevier Academic Press.

Barbosa, V. C. F., and J. B. C. Silva, 2011, Reconstruction of geologic bodies in depthassociated with a sedimentary basin using gravity and magnetic data: GeophysicalProspecting, 59, 1021–1034.

Kelley, C. T., 1999, Iterative methods for optimization: Raleigh: SIAM.Martins, C. M., W. A. Lima, V. C. F. Barbosa, and J. B. C. Silva, 2011, Total variation

regularization for depth-to-basement estimate: Part 1 - mathematical details andapplications: Geophysics, 76, I1–I12.

Medeiros, W. E., and J. B. C. Silva, 1996, Geophysical inversion using approximateequality constraints: Geophysics, 61, 1678–1688.

Menke, W., 1989, Geophysical data analysis: Discrete inverse theory: San Diego:Academic Press Inc.

Silva, J. B. C., W. E. Medeiros, and V. C. F. Barbosa, 2001, Potential field inversion:Choosing the appropriate technique to solve a geologic problem: Geophysics, 66,511–520.

Tarantola, A., 2005, Inverse problem theory and methods for model parameter estima-tion: Philadelphia: SIAM.

37

Apêndice A

Operações com matrizes

A seguir, demonstraremos como fazer algumas operações matemáticas com matrizes evetores. Mas antes, algumas definições.

A.1 Definições

Definição 1. Um vetor x de M elementos é uma matriz de uma coluna e M linhas

x =

x1

x2

...xM

(A.1)

Definição 2. Dado um conjunto de N funções f1(x), f2(x), . . . , fN(x), o vetor defunções f(x) é dado por

f(x) =

f1(x)

f2(x)...

fN(x)

(A.2)

Definição 3. A derivada do vetor f(x) de N elementos em relação ao i-ésimo elementode x, xi, é

A. Operações com matrizes 38

∂f(x)

∂xi=

∂f1(x)

∂xi

∂f2(x)

∂xi

...

∂fN(x)

∂xi

(A.3)

Definição 4. O vetor uNi de N elementos possui todos seus elementos iguais a zero,exceto o i-ésimo elemento que é igual a 1

uNi =

u1

...ui−1

ui

ui+1

...uN

=

0...0

1

0...0

(A.4)

Definição 5. O operador gradiente ∇ em relação ao vetor x de N elementos é igual a

∇ =

∂x1

∂x2

...

∂xN

(A.5)

Definição 6. O operador Hessiana ¯∇ em relação ao vetor x de N elementos é igual a

¯∇ = ∇∇T =

∂2

∂x21

∂2

∂x2∂x1

. . .∂2

∂xN∂x1

∂2

∂x1∂x2

∂2

∂x22

. . .∂2

∂xN∂x2

...... . . . ...

∂2

∂x1∂xN

∂2

∂x2∂xN. . .

∂2

∂x2N

(A.6)

A. Operações com matrizes 39

A.2 Derivadas

A seguir demonstramos como calcular ∂f(x)∂xj

e ∂f(x)∂xj

para diversos casos. Em todos oscasos, x é um vetor de M elementos e f(x) é um vetor de N elementos.

Exemplo A.1. Para o caso

f(x) = ¯Ax, (A.7)

em que ¯A é uma matriz de dimensão N ×M .Seja aj a j-ésima coluna de ¯A, podemos escrever a matriz ¯A em relação a suas M

colunas

¯A =[a1 . . . aj . . . aM

], (A.8)

e então

f(x) = x1a1 + · · ·+ xj aj + · · ·+ xM aM . (A.9)

Neste caso, a derivada de f(x) em relação a xj é

∂f(x)

∂xj=

0∂x1a1

∂xj+ · · ·+ ∂xj aj

∂xj+ · · ·+

>

0∂xM aM∂xj

= aj = ¯A uMj .

(A.10)

Exemplo A.2. Para o caso

f(x) = xT ¯AT , (A.11)

em que ¯A é uma matriz de dimensão N ×M .Seja aj a j-ésima coluna de ¯A, podemos escrever a matriz ¯AT em relação as M

colunas de ¯A

¯AT =

aT1...aTj...aTM

, (A.12)

e então

A. Operações com matrizes 40

f(x) = x1aT1 + · · ·+ xj a

Tj + · · ·+ xM a

TM . (A.13)

Neste caso, a derivada de f(x) em relação a xj é

∂f(x)

∂xj=

0∂x1a

T1

∂xj+ · · ·+

∂xj aTj

∂xj+ · · ·+

>

0∂xM aTM∂xj

= aTj = (uMj )T ¯AT .

(A.14)

Exemplo A.3. Para o caso

f(x) = x, (A.15)

a derivada de f(x) em relação a xj é

∂f(x)

∂xj= uMj . (A.16)

Exemplo A.4. Para o caso

f(x) = aT x = xT a, (A.17)

a derivada de f(x) em relação a xj é 1

∂f(x)

∂xj=

0∂x1a1

∂xj+ · · ·+ ∂xjaj

∂xj+ · · ·+

>

0∂xMaM∂xj

= aj = aT uMj . (A.18)

Exemplo A.5. Para o caso

f(x) = xT ¯AT ¯Ax = ( ¯Ax)T ( ¯Ax), (A.19)

a derivada de f(x) em relação a xj é (ver equação A.10)

∂f(x)

∂xj=

[∂( ¯Ax)

∂xj

]T( ¯Ax) + ( ¯Ax)T

[∂( ¯Ax)

∂xj

]

= ( ¯AuMj )T ( ¯Ax)︸ ︷︷ ︸escalar

+ ( ¯Ax)T ( ¯AuMj )︸ ︷︷ ︸escalar

.

(A.20)

Como o transposto de um escalar é igual a ele mesmo1O transposto de um escalar é igual a ele mesmo.

A. Operações com matrizes 41

∂f(x)

∂xj= 2( ¯AuMj )T ( ¯Ax) = 2(uMj )T ¯AT ¯A x . (A.21)

A.3 Gradientes

A seguir demonstramos como calcular ∇f(x) para diversos casos. Em todos os casos,x é um vetor de M elementos.

Exemplo A.6. Para o caso

f(x) = aT x = xT a, (A.22)

o gradiente de f(x) é (ver equação A.18) 2

∇f(x) =

∂f(x)

∂x1

∂f(x)

∂x2

...

∂f(x)

∂xM

=

a1

a2

...aM

= a . (A.23)

Exemplo A.7. Para o caso

f(x) = xT ¯AT ¯Ax, (A.24)

o gradiente de f(x) é (ver equação A.21)

∇f(x) =

∂f(x)

∂x1

∂f(x)

∂x2

...

∂f(x)

∂xM

=

2(uM1 )T ¯AT ¯A x

2(uM2 )T ¯AT ¯A x

...

2(uMM)T ¯AT ¯A x

= 2

(uM1 )T

(uM2 )T

...

(uMM)T

︸ ︷︷ ︸

¯I

¯AT ¯A x = 2 ¯AT ¯A x . (A.25)

2O transposto de um escalar é igual a ele mesmo.

A. Operações com matrizes 42

A.4 Hessianas

A seguir demonstramos como calcular ¯∇f(x) para diversos casos. Em todos os casos,x é um vetor de M elementos.

Exemplo A.8. Para o caso

f(x) = xT ¯AT ¯Ax, (A.26)

a Hessiana de f(x) é (ver equação A.25)

¯∇f(x) = ∇[∇f(x)]T = ∇(

2 ¯AT ¯A x)T

= 2∇(xT ¯AT ¯A︸︷︷︸

¯BT

)= 2∇

(xT ¯BT

). (A.27)

Sendo bj a j-ésima coluna de ¯B, segue que (ver equação A.14)

∇(xT ¯BT

)=

∂(xT ¯BT

)∂x1

∂(xT ¯BT

)∂x2

...

∂(xT ¯BT

)∂xM

=

bT1

bT2

...

bTM

= ¯BT . (A.28)

Logo,

¯∇f(x) = 2 ¯BT = 2 ¯AT ¯A . (A.29)

A.5 Expansão em série de Taylor

A seguir demonstramos como calcular a expansão em série de Taylor de f(x) até aprimeira e segunda ordem. Em todos os casos, x é um vetor de M elementos. Aexpansão será feita em torno de um ponto x0 tal que

x = x0 + ∆x . (A.30)

A. Operações com matrizes 43

A.5.1 Até primeira ordem

A expansão em série de Taylor de f(x) até a primeira ordem é

f(x) ≈ f(x0) +∂f

∂x1

(x0)∆x1 +∂f

∂x2

(x0)∆x2 + · · ·+ ∂f

∂xM(x0)∆xM

≈ f(x0) +

[∂f

∂x1

(x0)∂f

∂x2

(x0) . . .∂f

∂xM(x0)

]︸ ︷︷ ︸

∇f(x0)T

∆x1

∆x2

...∆xM

︸ ︷︷ ︸

∆x

≈ f(x0) + ∇f(x0)T∆x .

(A.31)

A.5.2 Até segunda ordem

Baseado na equação A.31, a expansão em série de Taylor de f(x) até a segunda ordemé

f(x) ≈f(x0) + ∇f(x0)T∆x

+1

2∆x1

∂2f(x0)

∂x21

∆x1 +1

2∆x1

∂2f(x0)

∂x2∂x1

∆x2 + · · ·+ 1

2∆x1

∂2f(x0)

∂xM∂x1

∆xM

+1

2∆x2

∂2f(x0)

∂x1∂x2

∆x1 +1

2∆x2

∂2f(x0)

∂x22

∆x2 + · · ·+ 1

2∆x2

∂2f(x0)

∂xM∂x2

∆xM

+ · · ·

+1

2∆xM

∂2f(x0)

∂x1∂xM∆x1 +

1

2∆xM

∂2f(x0)

∂x2∂xM∆x2 + · · ·+ 1

2∆xM

∂2f(x0)

∂x2M

∆xM .

(A.32)

Transformando as multiplicações de ∆xj do lado direito em um produto escalar com ovetor ∆x obtemos

A. Operações com matrizes 44

f(x) ≈f(x0) + ∇f(x0)T∆x

+1

2

[∆x1

∂2f(x0)

∂x21

∆x1∂2f(x0)

∂x2∂x1

. . . ∆x1∂2f(x0)

∂xM∂x1

]∆x

+1

2

[∆x2

∂2f(x0)

∂x1∂x2

∆x2∂2f(x0)

∂x22

. . . ∆x2∂2f(x0)

∂xM∂x2

]∆x

+ · · ·

+1

2

[∆xM

∂2f(x0)

∂x1∂xM∆xM

∂2f(x0)

∂x2∂xM. . . ∆xM

∂2f(x0)

∂x2M

]∆x .

(A.33)

Em seguida, colocamos ∆xj em evidência

f(x) ≈f(x0) + ∇f(x0)T∆x

+1

2∆x1

[∂2f(x0)

∂x21

∂2f(x0)

∂x2∂x1

. . .∂2f(x0)

∂xM∂x1

]∆x︸ ︷︷ ︸

escalar

+1

2∆x2

[∂2f(x0)

∂x1∂x2

∂2f(x0)

∂x22

. . .∂2f(x0)

∂xM∂x2

]∆x︸ ︷︷ ︸

escalar

+ · · ·

+1

2∆xM

[∂2f(x0)

∂x1∂xM

∂2f(x0)

∂x2∂xM. . .

∂2f(x0)

∂x2M

]∆x︸ ︷︷ ︸

escalar

.

(A.34)

Agora, transformamos a multiplicação de ∆xj a esquerda em um produto escalar como vetor ∆x

A. Operações com matrizes 45

f(x) ≈f(x0) + ∇f(x0)T∆x

+1

2∆xT

[∂2f(x0)

∂x21

∂2f(x0)

∂x2∂x1

. . .∂2f(x0)

∂xM∂x1

]∆x[

∂2f(x0)

∂x1∂x2

∂2f(x0)

∂x22

. . .∂2f(x0)

∂xM∂x2

]∆x

...[∂2f(x0)

∂x1∂xM

∂2f(x0)

∂x2∂xM. . .

∂2f(x0)

∂x2M

]∆x

. (A.35)

Colocando ∆x em evidência e lembrando da equação A.6 obtemos

f(x) ≈f(x0) + ∇f(x0)T∆x

+1

2∆xT

∂2f(x0)

∂x21

∂2f(x0)

∂x2∂x1

. . .∂2f(x0)

∂xM∂x1

∂2f(x0)

∂x1∂x2

∂2f(x0)

∂x22

. . .∂2f(x0)

∂xM∂x2...

... . . . ...∂2f(x0)

∂x1∂xM

∂2f(x0)

∂x2∂xM. . .

∂2f(x0)

∂x2M

︸ ︷︷ ︸

¯∇f(x0)

∆x. (A.36)

Finalmente, a expansão de f(x) em série de Taylor até segunda ordem é

f(x) ≈ f(x0) + ∇f(x0)T∆x+1

2∆xT ¯∇f(x0)∆x . (A.37)