Algoritmo de Levenberg-Mardquardt

El algoritmo de Levenberg-Marquardt (LM) es un algoritmo iterativo

de optimizacin en el que el mtodo de iteracin presenta una ligera
modificacin sobre el mtodo tradicional de Newton. Las ecuaciones
normales N=J
T
J=J
T
(J representa el jacobiano de la funcin, los
incrementos de los parmetros y el vector de errores residuales del
ajuste), son reemplazadas por las ecuaciones normales aumentadas
N=J
T
, donde N
ii
=(1+
i
) N
ii
y N
ii
= N
ii
para ij. El valor de es
inicialmente puesto a algn valor, normalmente =10
-3
. SI el valor de
obtenido resolviendo las ecuaciones aumentadas conduce a una
reduccin del error, entonces el incremento es aceptado y es
dividido por 10 para la siguiente iteracin. Por otro lado si si el valor
de conduce a in aumento del error, entonces es multiplicado por
10 y se resuelven de nuevo las ecuaciones normales aumentadas,
este proceso contina hasta que el valor de encontrado da lugar
aun decremento del error. Este proceso de resolver repetidamente las
ecuaciones normales aumentadas para diferentes valores de hasta
encontrar un valor aceptable de es lo que constituye una iteracin
del algoritmo de LM.

Nota de implementacin: Para usar el algoritmo de LM en la
minimizacin de una funcin solo es necesario dar una rutina que
calcule la funcin a minimizar, un vector de observaciones o valores
deseados de la funcin y una estimacin inicial. El clculo de la matriz
jacobiana se puede llevar a cabo de forma tanto de forma numrica
como dando una funcin que la calcule..

En el aso de diferenciacin numrica, cada variable independiente x
i

se incrementa por turnos en x
i
+, se calcula es valor de la funcin en
el nuevo punto y la derivada se calcula como un cociente. Buenos
resultados han sido encontrados poniendo el valor de al mximo
entre |10
-4
x
i
| y 10
-6
. En la practica no se aprecia ventaja en usar un
mtodo de diferenciacin numrica o dar una rutina de clculo de la
derivada.

Versin para matrices dispersas
La versin antes descrita es adecuada para la minimizacin con
respecto a un nmero pequeo de parmetros. Por ejemplo en el
caso de estimacin de una homografa, cuando el vector de
parmetros son los elementos de la matriz de la homografa. Sin
embargo, para minimizar funciones con un alto nmero de
parmetros el algoritmo clsico no es muy adecuado. LA razn
principal es que la resolucin de las ecuaciones normales tiene
complejidad N
3
en el nmero de parmetros y este paso se repite
muchas veces. Sin embargo en la solucin de muchos de los
problemas que se estudian en loe problemas de la geometra de
mltiples vistas, la matriz de las ecuaciones normales presenta un
cierto estructura de bloques dispersos de la cual podemos ahcer uso
para alcanzar un ahorro considerable en el tiempo de clculo.

Un ejemplo de esta situacin es el caso de la estimacin de una
homografa en donde suponemos errores en amabas imgenes y
minimizamos el error de mxima verosimilitud. En este caso la
parametrizacin del problema esta constituida por 2n+9 parmetros,
siendo n el nmero de correspondencias usado.

Otro ejemplo donde estos mtodos son tiles es en le probelam de la
reconstruccin en el cual tenemos correspondencias a travs de dos o
ms imgenes (vistas) y deseamos estimar los parmetros de las
cmaras de todas las cmaras al mismo tiempo que la posicin 3D de
todos los puntos. Podemos asumir cmaras proyectivas cualesquiera
o cmaras parcial o totalmente calibradas. En este caso el problema
puede ser parametrizado por las entradas de todas las cmaras ( es
deir 11m o 12m dependiendo de la parametrizacin de las cmaras),
mas 3n parmetros para las coordenadas de los puntos 3D.

A continuacin se dan algunos resultados que permiten una
implementacin simple de este algoritmo haciendo uso de una librera
de clculo matricial.

Ahora un vector representado por la notacin regruesada, a, ser un
vector fila.

Particin del vector de parmetros

Como ya hemos visto los problemas que presentan matrices
dispersas por bloques hacen uso de una parametrizacin en la cual se
pueden distinguir dos tipos distintos de parmetros ( las matrices y
los puntos). As pues cabe establecer una versin del algoritmo
pensando en una divisin del vector de parmetros en dos bloques.

Sea PR
M
el vector de parmetros que suponemos se puede dividir
en dos vectores a y b tal que P=(a
T
,b
T
)
T
. Suponemos que tenemos
un vector de medidas X en el espacio R
N
. Adems sea
X
la matriz de
covarianzas para el vector de medidas. Consideramos una funcin
general f: R
M
R
N
que define la estimacin del vector de medidas
=f(P). Notando por = X- el error entre las observaciones y las
estimaciones lo que buscamos es el conjunto de parmetros que
minimizan el cuadrado de la distancia de Mahalanobis
X

1
X
T

=
2
X
.

Correspondiendo a la divisin de parmetros P=(a
T
,b
T
)
T
, la matroiz
jacobiana J=[ /P] tiene una estructura de bloque de la forma
J=[A|B] donde las submatrices jacobianas estan definidas por
X

A=[ /a], B=[ /b] X

El conjunto de ecuaciones J= del paso central del algoritmo de LM
toma la forma
J=[A|B] (
b
a

)=
Entonces las ecuaciones normales toman la forma

|
|
.
|

\
|
=
|
|
.
|

\
|
(

B
A

B B A B
B A A A
1
X
1
X a
1
X
1
X
1
X
1
X
T
T
T T
T T
b


Justo en este paso es donde los bloque diagonales son aumentados
multiplicando sus valores diagonales por un factor 1+. Este cambio
modifica las matrices de la diagonal principal. Las matrices
resultantes se notan por (A
X
-1
A)
*
y (B
X
-1
B)
*
.

El conjunto de ecuaciones anterior se puede escribir como
|
|
.
|

\
|
=
|
|
.
|

\
|
|
|
.
|

\
|
b
a
b
a

V W
W U
*
*
T


Eliminando el bloque superior izquierdo a travs de la multiplicacin
por la matriz

|
obtenemos
|
|
.
|

\


I 0
WV I
1 *

|
|
.
|

\
|
=
|
|
.
|

\
|
|
|
.
|

\
|

b
b a
b
a

WV

V W
0 W WV U
1 1 *
*
* *
T
T


De esta ecuacin puede despejarse y calcularse el valor de
a

posteriormente el de
b
.

Una vez que hemos calculado los incrementos de los parmetros
actuamos igual que en el algoritmo clsico tratando de averiguar si el
nuevo vector P=((a+
a
)
T
, (b+
b
)
T
)
T
reduce o no el error y en cada caso
actuando como ya se ha establecido.

Aunque en este mtodo se resuelve primero para
a
y luego para
b

basndonos en el nuevo valor de a, hemos de insistir que esto no
quiere decir que el mtodo conduzca a iterar sobre a y b de forma
independiente. Si se desea se podran hacer iteraciones
independientes para a y b pero la convergencia sera sin duda ms
lenta.
Matriz de Covarianzas

Hemos visto que la matriz de covarianzas de los parmetros
estimados esta dada por

P
=(J
T

X
-1
J)
+

En el caso de sobre-parametrizacin, esta matriz de covarianzas es
singular, y en particular, a los parmetros no se les permite variar en
las direcciones perpendiculares a la restriccin a la superficie, la
varianza es cero en esas direcciones.

En el caso de parmetros separados P=(a
T
,b
T
)
T
la matriz J
T

X
-1
J la
escribimos en la forma

(

=
(

V W
W U
B B A B
B A A A
J J
1
X
1
X
1
X
1
X
T T T
T T
1 T
X


Por tanto la matriz de covarianzas es la pseudoinversa de esta
matriz. Bajo la hiptesis de que V es invertible, redefinimos Y=WV
-1
.
Entonces la matriz puede ser diagonalizada de acuerdo a

(

=
(

I Y
0 I
V 0
0 W WV U
I 0
Y I
V W
W U
J J
T
1 T
T
1 T
T
X


Asumiendo que se verifica la propiedad (GHG
T
)
+
=G
-T
H
+
G
-1
para
matrices G y H con G invertible, se puede dar una expresin final
para la pseudo-inversa

(

= =

1 T T
1 T
V XY Y X Y
XY X
J J
P X

donde X=(U - WV
-1
W
T
)
+
.

Este resultado nos permite tambin obtener las matrices de
covarianza de los parmetros a y b de forma separada. (
a
=X ,

b
=Y
T

a
Y+V
-1
).

Mtodo de LM-Disperso General

A continuacin vamos a ver que la generalizacin de la tcnica que
acabamos de describir puede ser usada con ventaja en un caso
general sin la matriz jacobiana presenta una condicin de dispersin
particular que analizaremos.

Supongamos que el vector de medidas XR
N
puede ser partido en
trozos como X=( X
1
T
, X
2
T
, ..., X
n
T
)
T
, de forma similar suponemos que el
vector de parmetros PR
M
, puede dividirse en P=(a
T
, b
1
T
,b
2
T
,...,b
n
T
)
T
.
Los valores estimados de X
i
correspondientes a un conjunto de
parmetros dado se notar por X .
i

Haremos la hiptesis de dispersin que cada solo depende de a y
i
X

b
i
pero de ninguno de los otros b
j
para ji. En este caso /b
i
X

i
=0
para ij. No se hace ninguna hiptesis sobre /a. Esta situacin se
presenta en los problemas de reconstruccin ya comentados con
anterioridad, ya que en este caso b
i
X

i
es el vector de parmetros del i-
simo punto, y es el vector de medidas de la imagen de este
punto en todas las vistas. En este caso ya que la imagen de un punto
no depende de ningn otro punto se verifica la hiptesis de
dispersin. Definimos las matrices jacobianas por
i
X

A
i
=[ /a], B
i
X

i
=[ /b
i
X

i
]

Dado un vector de errores de la forma =(
1
T
,
2
T
, ...,
n
T
)
T
=X- el
conjunto de ecuaciones J= tiene ahora la forma
X

|
|
|
.
|

\
|
=
|
|
|
|
|
.
|

\
|
(
(
(
(

=
n
1

B A
B A
B A
J M
M O M
n
b
b
a
n n
1 2 2
1 1

Suponemos adems que todas las medidas X
i
son independientes con
matrices de covarianza
X
= diag( , ..., ).
1
X

2
X

n
X

En la notacin del algoritmo anterior tenemos

A=[A
1
T
, A
2
T
,...,A
n
T
]
T
B=diag(B ,B ,...,B )

X
= diag( , ..., )
1
X

2
X

n
X

b
=( , , , )
T
1
b
T
2
b

T
n
b

T
=(
1
T
,
2
T
, ...,
n
T
)
T


Ahora es una tarea directa sustituir estas formulas en el algoritmo
anterior. El resultado es el algoritmo que se presenta a continuacin.
Es importante resaltar que de esta forma cada paso del algoritmo
requiere tiempo lineal en n. Sin la ventaja resultante de la estructura
dispersa el algoritmo tendra una complejidad de n
3
.





Algoritmo LM para matrices dispersas (2 vistas)

1. Calcular las matrices derivadas A
i
=[ /a], B
i
X

i
=[ /b
i
X

i
] y el
vector de errores
i
=X
i
- .
i
X

2. Calcular los valores intermedios

i
i
X
T
i
i
A A U


=
1

( )
n 2 1
V , , V , V V L diag = donde
i X
T
i
i
B B V
i
1
=
( )
n 2 1
W , , W , W W L = donde
i X
T
i
i
B A W
i
1
=
i
i
X
T
i
i
A
A

=
1

( )
T
T
B
T
B
T
B
n 2 1

B
, , , L = donde
i X
T
i
i
B
1
B
i

=
Y
i
=W
i
V
i
*-1


3. Calcular
a
a partir de la ecuacin
i
i
i a
i
T
i i B A
Y W Y U

=
|
.
|

\
|

*


4. Calcular por turnos cada
bi
a partir de la ecuacin
) (
1
a
T
i B b
W V
i i
=
*
i


Covarianza

1.- Redefinir Y
i
=W
i
V
i
-1

2.-
a
=(U-

)
i
T
i i
W Y
+

3.- =Y
j
b b
i
i
T

a
-1
Y
j
+
ij
V
i
-1
4.-
i
Y
a ab
i
=



Aplicaciones a la estimacin de una homografa y la matriz
fundamental

El algoritmo dado puede ser aplicado de forma directa al caso de la
estimacin por mxima verosimilitud de una homografa o de la
matriz fundamental sin ms que calcular los valores de las matrices
derivadas A
i
y B
i
y vista su estructura tratar de obtener la mxima
ventaja de la misma a la hora de realizar los clculos.

En el caso de la matriz fundamental en el que la parametrizacin se
realiza a travs de la matriz de segunda cmara y los valores 3D, es
posible obtener la matriz de covarianza de F a travs de la expresin
que relaciona F y P=[M|m] . Supuesto que ||F||=1, podemos expresar
F=[m]
x
M/(||[m]
x
M||) lo que nos permite definir F como una funcin
sencilla de P y calcular el jacobiano de la transformacin. Entonces
obtenemos

F
=J
P
J
T

y
P
dada por el algoritmo anterior.

Algoritmo LM para matrices dispersas (n vistas)

Supongamos ahora que disponemos de n>2 vistas simultaneas y
correspondencias de puntos a travs de ellas. Por simplicidad en la
notacin supondremos que los puntos aparecen en todas las vistas
(cosa que no tiene que ser verdad). Esta extensin es valida tanto
para la estimacin de los tensores de orden superior (trifocal y
cuadrifocal) como para llevar a cabo un ajuste global de n imgenes

Ahora tomaremos ventaja no solo de la ausencia de interaccin de
unos puntos con otros si no tambin de la ausencia de interaccin de
los parmetros de unas cmaras con otras. Usando la misma
notacin del algoritmo anterior, el vector de medidas X puede
descomponerse en subvectores X
i
=(x
i1
T
, x
i2
T
, ..., x
im
T
)
T
donde x
ij
es la
imagen del i-simo punto en la j-sima imagen. El vectyor de los
parmetros de la cmara puede tambin partirse en m sub-vectores
representando los parmetros de cada una de las cmaras a=(a
1
T
,
a
2
T
, ...,a
m
T
)
T
donde a
j
son los parmetros de la j-sima cmara.

Ahora teniendo en cuanta las hiptesis de dispersin asumida
/a
ij
x
k
=0 a menos que j=k y /b
ij
x
k
=0 a menos que i=k.

Las matrice jacobianas definidas en el algoritmo de 2-vistas ahora
son A
i
=[ /a] es una matriz diagonal por bloques
i
X

A
i
=diag(A
i1
, A
i2
,..., A
im
) donde A
ij
= /a
ij
x
j

De forma similar, la matriz B
i
=[ /b
i
X

i
] se descompone como
B
i
=[B
i1
T
, B
i2
T
,..., B
im
T
]
T
donde B
ij
= /b
ij
x
i

Respecto de la matriz de covarianzas puede suponerse que tiene una
estructura diagonal

Xi
= diag( , ..., )
1 i
X

2 i
X

in
X

dando por bueno que las medidas de distintas imgenes son

incorreladas.

Usando estas expresiones el algoritmo de dos vistas puede adaptarse
fcilmente al caos de n-vistas.