M´ınimos quadrados e extensoe˜ s - NACAD/COPPE-UFRJamit/alglin/aulals.pdf · 2010. 4. 15. · Invers˜aoMQ:Escolhercomoestimativa bx queminimiza kAbx−bk: i.e., o desvio entre

Ḿınimos quadrados e extensões

Amit Bhaya,Programa de Engenharia Elétrica

COPPE/UFRJUniversidade Federal do Rio de Janeiro

[email protected]

http://www.nacad.ufrj.br/̃ amit

May 9, 2006

MPC 2006, Amit Bhaya Aula MQ e extensões

Métodos de ḿınimos quadrados

• Soluções aproximadas (min. quad.) de sistemassobredeterminados (A AM)

• Estimação no sentido de MQ

• Ajuste de dados no sentido de MQ

• Soluções de problemas de MQ via fatoração QR

• Aplicação a identificação de sistemas

1


Sistemas lineares sobredeterminados

Considere Ax = b, onde A ∈ Rm×n é alta e magra(m > n). Este é um sistema sobredeterminado (i.e.,mais equações do que incógnitas).

Para maioria de b, inexiste solução x (porque b 6∈R(A))

Um enfoque: Solução aproximada:• Define reśıduo or erro r := b − Ax• Objetivo: Achar x̂ = xls que minimiza ‖r‖

x̂ = xls é denominado solução (aproximada) no sentidode ḿınimos quadrados de Ax = b.

interpretação geométrica: Axls é o ponto em R(A)mais próximo a b – Axls é a projeção de b sobreR(A).

2


Equações normais via cálculo matricial

Hipótese A AM, com posto completo por colunas.Queremos achar xls: formar reśıduo ao quadrado.

‖r‖2 = xTATAx − 2bTAx + bTb

e depois derivar em relação a x e igualar a zero:

2xTATA − 2bTA = 0

levando à equação normal :

ATAxls = ATb

Pela hipótese (A AM posto completo por col.) ATAé inverśıvel, portanto xls = (A

TA)−1ATb• xls é função linear de b• xls = A−1b se A for quadrada.• xls é solução de b = Axls se b ∈ R(A)• A† := (ATA)−1AT é denominada pseudoinversa deA e é uma inversa à esquerda de A AM, p.c.c. (postocompleto por colunas).

3


Matriz projetora, prinćıpio de

ortogonalidade

Axls, a projeção de y sobre R(A), é linear, dada por

Axls = A(ATA)−1ATb

• A(ATA)−1AT é denominada matriz projetora.

Prinćıpio de ortogonalidade: Reśıduo ótimo

r = Axls − b = (A(ATA)−1AT − I)b

é ortogonal a R(A):

〈r,Az〉 = bT (A(ATA)−1AT − I)TAz = 0

para todo z ∈ Rn.

4


Estimação no sentido de MQ

Muitas aplicações em problemas de inversão,estimação, e reconstrução do tipo:

b = Ax + v

• x é o que se deseja estimar ou reconstruir• b é o vetor de medidas dos sensores• v é um rúıdo desconhecido ou erro de medida(suposto pequeno)• i-ésima linha de A caracteriza o i-ésimo sensor

Inversão MQ: Escolher como estimativa x̂ que minimiza‖Ax̂ − b‖: i.e., o desvio entre• o que realmente observamos (b), e• o que observaŕıamos se x = x̂, e não existisse rúıdo(v = 0)

Estimativa MQ é simplesmente x̂ = (ATA)−1ATb

5


Propriedade BLUE

Medidas lineares com rúıdo: b = Ax + v, com A AM, p.c.c.

Considere um estimador linear da forma x̂ = Bb.

É denominado unbiased se x̂ = x sempre que v = 0 (i.e., não há

erro de estimação quando não há ruido) coincide com BA = I,

i.e., B é inversa à esquerda de A

Erro de estimação de estimador linear unbiased (sem viés) é

x − x̂ = x − B(Ax + v) = −Bv

Evidentemente, gostaŕıamos de B ‘pequena’.

Fato: A† := (ATA)−1AT é a menor inversa a esquerda de A:Para qualquer B satisfazendo BA = I, temos

∑

i,j

b2ij ≥

∑

i,j

a†2ij

i.e., ḿınimos quadrados leva ao melhor estimador linear sem v́ıcio

(= best linear unbiased estimator (BLUE))

6


Regressores

Considere faḿılia de problemas de MQ

minimize

∥∥∥∥∥

p∑

i=1

xiai − b∥∥∥∥∥

para p = 1, . . . , n.Em estat́ıstica, os ai são denominados regressores.

Jargão equivalente:

1. Aproxime b utilizando combinação linear de a1, . . . , ap.

2. Projete b no subespaço EG{a1, . . . ap}3. Faça a regressão de b sobre a1, . . . , ap

Quando p aumenta, o ajuste é melhor, portanto reśıduo ótimo

diminui

Solução para cada p ≤ n é dada por xls = R−1p QTp ,onde Rp consiste na submatriz ĺıder p × p de R, andQp = [q1 · · ·qp] consiste nas primeiras p colunas de Q:i.e.. uma fatoração QR resolve uma faḿılia crescentede problemas MQ.

7


Identificação de sistemas discretos

utilizando MQ

Medimos entrada u(t) e sáıda y(t) para t = 0, . . . , Nde um sistema (SISO) desconhecidoProblema de identificação de sistema: Encontremodelo razoável para o sistema baseado nos dadosE/S u e y.Exemplo SISO (fácil extensão ao caso MIMO): ajustarmodelo MA (média móvel) com n atrasos aos dadosE/S:

ŷ(t) = h0u(t) + h1u(t − 1) + · · · + hnu(t − n)onde h0, . . . , hn ∈ RA sáıda prevista (= sáıda do modelo) pode ser escritocomo:

ŷ(n)ŷ(n + 1)

...ŷ(N)

=

u(n) u(n − 1) · · · u(0)u(n + 1) u(n) · · · u(1)

......

...u(N) u(N − 1) · · · u(N − n)

h0h1...

hn

Erro de predição do modelo é:

e = [y(n) − ŷ(n) · · · y(N) − ŷ(N)]T

8


identificação MQ: Escolher modelo (i.e., h) queminimiza norma do erro de predição ‖e‖ ... é umproblema MQ sobredeterminado...

9


Conjuntos crescentes de medidas

Problema de ḿınimos quadrados em termos de ‘linhas’:

minimize‖Ax − b‖2 =m∑

i=1

(aTi x − bi)2

onde aTi são as linhas de A.

• x ∈ Rn vetor a ser estimado.

• cada par (ai, bi) corresponde a uma medida.

• Solução é

xMQ =

(m∑

i=1

aiaTi

)−1 m∑

i=1

biai

• supondo que ai e bi se tornam dispońıveissequencialmente (isto é, m aumenta com tempo).

10


Ḿınimos quadrados recursivos

Podemos computar xMQ(m) =(∑m

i=1 aiaTi

)−1 ∑mi=1 biai

recursivamente.

• Inicialize P(0) = 0 ∈ Rn×n, q(0) = 0 ∈ Rn.

• Para m = 0, 1, . . . ,

P(m + 1) = P(m) + am+1aTm+1

q(m + 1) = q(m) + bm+1am+1.

• Se P(m) for inverśıvel, temos xMQ =P(m)−1q(m).

• P(m) inverśıvel ⇐⇒ EG{a1, . . . ,am} = Rn(portanto, uma vez que P(m) se torna inverśıvelele continua inverśıvel).

11


Atualização rápida para ḿınimos

quadrados

Podemos calcular

P(m + 1)−1 =(P(m) + am+1a

Tm+1

)−1

eficientemente, a partir de P−1 utilizando a fórmulade atualização de posto um

(P + aaT

)−1= P−1 − 1

1 + aTP−1a(P−1a)(P−1a)T

válido quando P = PT , e P e P + aaT ambasinverśıveis.

• fornece método O(n2) para computar P(m + 1)−1a partir de P(m)−1.

• Métodos padrão para computar P(m+1)−1 a partirde P(m)−1 são O(n3).

12


Verificação da fórmula de atualização de

posto um

(P + aaT )(P−1 − 1

1+aTP−1a(P−1a)(P−1a)T

)

= I + aaTP−1 − 11+aTP−1aP(P

−1a)(P−1a)T

− 11+aTP−1aaa

T (P−1a)(P−1a)T

= I + aaTP−1 − 11+aTP−1aaa

TP−1 − aTP−1a1+aTP−1aaa

TP−1

= I

13


Ḿınimos quadrados multiobjetivo

Em muitos problemas temos dois (ou mais) objetivosdefinidos como funções da variável de projeto x ∈ Rn:

• queremos J1 = ‖Ax − b‖2 pequeno.

• queremos, ao mesmo tempo, J2 = ‖Fx − g‖2pequeno.

• geralmente, objetivos são conflitantes.

• Podemos diminuir um em detrimento do outro.Exemplo comum: F = I, g = 0 queremos erropequeno com solução x de norma pequena.

14


Gráfico de pares atinǵıveis

J1

J2

x(1)

x(2)x(3)

Note que x ∈ Rn, mas este gráfico está em R2, o pontocom rótulo x(1) na verdade é (J2(x

(1)), J1(x(1))).

• área hachurada mostra (J2, J1) atinǵıvel por algum x ∈ Rn.

• área em branco mostra (J2, J1) não atinǵıvel por qualquer x ∈ Rn.

• fronteira da região é chamada curva de compromisso ótima

• Os x correspondentes são chamados de ótimos no sentido de Pareto (ótimos n.s.P.) (paraos objetivos J1, J2).

Três exemplos de escolha de x

• x(3) é pior do que x(2) em ambos os ı́ndices.

• x(1) é melhor do que x(2) em J2, porém pior em J1.

15


Função Objetivo de soma ponderada

• para achar pontos ótimos (i.e., x em cima dacurva de compromisso ótima) minimizamos a funçãoobjetivo de soma ponderada

J1 + µJ2 = ‖Ax − b‖2 + µ‖Fx − g‖2

• o parâmetro µ ≥ 0 especifica a ponderação relativaentre objetivo J1 e J2.

• Pontos onde a soma ponderada é constante (J1 +µJ2 = α) correspondem a uma reta com inclinação−µ no gráfico de J2 contra J1.

16


Minimização da objetivo soma ponderada

J1

J2

x(1)

x(2)x(3)

J1 + µJ2 = α

• x(2) minimiza objetivo de soma ponderada para o µescolhido no gráfico.

• Varrendo µ de 0 a +∞, podemos traçar a curva decompromisso ótima inteira.

17


Minimização de objetivo de soma

ponderada

Podemos expressar objetivo de soma ponderada comoMQ padrão:

‖Ax − b‖2 + µ‖Fx − g‖2 =∥∥∥∥

[A√µ F

]x −

[b√µg

]∥∥∥∥2

=∥∥∥Ãx − b̃

∥∥∥2

onde

Ã =

[A√µF

], b̃ =

[b√µg

].

Portanto, se Ã tiver posto completo, a solução é:

xMQ =(ÃT Ã

)−1ÃT b̃

= (ATA + µFTF)−1(ATb + µFTg)

18


Ḿınimos quadrados regularizados

Quando F = I e g = 0, os objetivos são:

J1 = ‖Ax − b‖2, J2 = ‖x‖2

O minimizador da função objetivo ponderada é:

x = (AAT + µI)−1ATb

também denominada solução (aproximada) regularizadano sentido de ḿınimos quadrados ou regularização deTikhonov.

Existe µ > 0 que pode ser utilizado para qualquerA (sem restrições sobre forma (AM/BG) e posto).

19


Aplicações de regularização

Estimação/inversão

• Ax − b é reśıduo de sensor

• Informação a priori: x pequeno

• ou, por exemplo, modelo confiável (preciso) apenaspara x pequeno

• Solução regularizada busca equilibrar ajuste desensor e tamanho de x.

20


Equações subdeterminados: soluções de

norma ḿınima

• Solução de norma ḿınima de sistema linearsubdeterminado

• Solução de norma ḿınima via fatoração QR

• derivação via cálculo (multiplicador de Lagrange)

• Relação com MQ regularizado

21


Solução de norma ḿınima

Considere Ax = b, onde A ∈ Rm×n BG (m < n),i.e.,

• há mais variáveis do que equações

• x é sub-especificado, i.e., muitas escolhas de xlevam ao mesmo b.

Supomos A de posto completo (= m), então para cada b ∈ Rm,o conjunto de todas as soluções é dado por:

{x : Ax = b} = {xp + z : z ∈ N (A)},

sendo xp qualquer solução “particular” (i.e., Axp = b).

• z caracteriza escolhas dispońıveis na solução• Solução possui dim N (A) ‘graus de liberdade’• Podemos escolher z para satisfazer outros critérios ou como

parâmetro a ser otimizado

22


Solução de norma ḿınima

Uma solução particular é:

xnm = AT (AAT )−1b

(AAT é inverśıvel, pois A possui posto completo.

De fato, xnm é a solução de b = Ax que minimiza‖x‖. Em outras palavras, xnm é a solulção do problemade minimização

minimize ‖x‖sujeito a Ax = b, x ∈ Rn

Seja Ax = b, então A(x − xnm) = 0, e

(x − xnm)T xnm = (x − xnm)T AT (AAT )−1b= (A(x − xnm))T (AAT )−1b= 0

i.e., (x − xnm) ⊥ xnm, então

‖x‖2 = ‖xnm + x − xnm‖2 = ‖xnm‖2 + ‖x − xnm‖2 ≥ ‖xnm‖2

i.e., xnm é, de fato, a solução de norma ḿınima.

23


Solução de norma ḿınima: geometria

{x : Ax = b}

N (A) = {x : Ax = 0}

xnm

0

• condição de ortogonalidade: xnm ⊥ N (A)

• Interpretação de projeção: xnm é a projeção de 0no conjunto solução {x : Ax = b.

24


Solução de norma ḿınima, ḿınimos

quadrados: fórmulas

• A† = AT (AAT )−1 é chamado pseudo-inversa deA (BG, posto completo por linha)

• AT (AAT )−1 é uma inversa a direita de A.

• I − AT (AAT )−1A realiza projeção sobre N (A).

• A† = (ATA)−1AT é a fórmula análoga para A(AM, posto completo por coluna)

• (ATA)−1AT é uma inversa a esquerda de A.

• A(ATA)−1AT realiza projeção sobre R(A)

25


Solução de norma ḿınima via fatoração

QR

Ache fatoração QR de AT , i.e., AT = QR tal que

• Q ∈ Rn×m, QTQ = Im

• R ∈ Rm×m triangular superior, não singular.

Então

• xnm = AT (AAT )−1b = QR−Tb

• ‖xnm‖ = ‖R−Tb‖

Obs.: R−T denota (R−1)T

26


Derivação via multiplicadores de

Lagrange

• xnm é solução do problema de otimização minxTx,sujeito a Ax = b

• Introduzindo multiplicadores de Lagrange: L(x, λ) =xTx + λT (Ax − b)

• Condições de otimalidade são:

∇xL = 2x + ATλ = 0, ∇λL = Ax − b = 0

• Da 1a. condição, x = −(1/2)ATλ

• Substituindo na segunda equação, obtemos: λ =−2(AAT )−1b

• Portanto, x = AT (AAT )−1b

27


Solução de norma ḿınima via

regularização

• Seja A ∈ Rm×n BG, posto completo por linha.

• Seja J1 = ‖Ax − b‖2, J2 = ‖x‖2

• Solução de norma ḿınima minimiza J2 com J1 = 0.

• minimizador da soma ponderada J1+µJ2 = ‖Ax−b‖2 + µ‖x‖2 é

xµ = (ATA + µI)−1ATb

• Fato (mostre!): xµ → xnm quando µ → 0, i.e.,solução regularizada converge a solução de normaḿınima quando µ → 0

• Em termos matricias, quando µ → 0, para A BG ede posto completo

(ATA + µI)−1AT → AT (AAT )−1

28

Documents

M´ınimos quadrados e extensoe˜ s - NACAD/COPPE-UFRJamit/alglin/aulals.pdf · 2010. 4. 15. · Invers˜aoMQ:Escolhercomoestimativa bx queminimiza kAbx−bk: i.e., o desvio entre