Re: duda
jros
> Esta es una inversa generalizada ?
>
> [(L'*L)^(-1)]*L'
>
>
> gracias
>
Si tienes el sistema A*x=b, y suponemos que el sistema es incompatible,
por ejemplo porque A tiene más filas que columnas, y b no es una
combinación lineal de las columnas de A (es decir no tiene solución).
Lo que puedes hacer es transformarlo en el sistema
A'*A*x=A'*b (notación de MatLab)
Este sistema si que tiene solución, aunque evidentemente no es solución
del sistema anterior, ahora si la solución de este sistema minimiza la
norma cuadrática de la del vector A*x-B (es decir el error de la
solución del sistema) como veremos más adelante.
La solución del sistema anterior es
x=(A'*A)^(-1)*A'*b
Si el sistema hubiese sido compatible determinado, la solución hubiese
sido
x=A^(-1)*b
Así por comparación de las dos últimas ecuaciones se decide llamar a
(A'*A)^(-1)*A'
Es importante notar que
(A'*A)^(-1)*A' * A = I (matriz identidad)
pseudoinversa. Es más pseudoinversa por la izquierza, ya que puede
plantearse también una por la derecha que creo que es para la solución
del mismo problema expresado en forma transpuesta y*B=c, donde B tiene
má columnas que filas e y y c son vectores horizontales. Está
pseudoinversas por la izquierda y por la derecha se denominan "inversas
por un lado" (one-sided inverse).
El porque se minimiza el error cuadrático (A*x-b)'*(A*x-b) es porque
La solución que minimiza el error es evidentemente
d (A*x-b)'*(A*x-b)/ d x = 0
A'*(A*x-b)=0 => A'*A*x=A'*b es decir el problema que hemos resuelto.
octave:1> A=[1 4; 2 5; 3 6]
A =
1 4
2 5
3 6
octave:2> AplusL=(A'*A)^(-1)*A'
AplusL =
-0.94444 -0.11111 0.72222
0.44444 0.11111 -0.22222
octave:3> AplusL*A
ans =
1.0000e+00 3.1086e-15
-7.7716e-16 1.0000e+00
octave:4>
En todo caso si el problema es compatible determinado, la pseudoinversa
por la izda coincide con la inversa. Si el problema es compatible
indeterminado (deficiente en rango), la pseudoinversa por la izda no
existe ya que A'*A no tiene inversa.
octave:4> A=A'
A =
1 2 3
4 5 6
octave:5> A'*A
ans =
17 22 27
22 29 36
27 36 45
octave:6> inv(A'*A)
warning: inverse: matrix singular to machine precision, rcond =
8.07128e-18
ans =
1.4340e+14 -2.8680e+14 1.4340e+14
-2.8680e+14 5.7359e+14 -2.8680e+14
1.4340e+14 -2.8680e+14 1.4340e+14
En el caso anterior algoritmo QR podría utilizarse para buscar una de
las infinitas soluciones.
También podría utilizarse la inversa por la derecha para buscar la
solución a x', del sistema transpuesto x'*A'=b':
Así si B=A', B'*(B*B')^(-1) es la pseudoinversa por la derecha y por
tanto
x'=b'*B'*(B*B')^(-1)
será la solución, que aunque no lo he mirado tendrá alguna
interpretación (mirar que minimiza queda propuesto por ejercicio yo no
lo he pensado).
octave:7> A*A'*(A*A')^(-1)
ans =
1.0000e+00 -1.7764e-15
4.5519e-15 1.0000e+00
Por último en la bibliografía se definen varias pseudoinversas. Por
ejemplo la pseudoinversa de Moore-Penrose
(http://en.wikipedia.org/wiki/One-sided_inverse), coincide con la
pseudoinversa por la izda cuando hay más filas que columnas y con la
pseudoinversa por la derecha cuando hay más columnas que filas. Además
cuando numero de filas y columnas es igual coincide con la inversa
estandar (supongo matrices con rango máximo).
Puedes leerlo leerlo en
http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse#Full_rank
Nota que el * es equivalente a la transpuesta si estamos con números
reales.
La notación para la pseudoinversa de $A$ es normalmente $A^+$
(utilizando código LaTeX) independientemente del tipo de pseudoinversa
de que se trate, así que es importante el contexto.
Cuando existe la inversa por la izda (y la matriz no es inversible), no
existe la inversa por la derecha (y viceversa).
La inversa de Moore-Penrose existe siempre, y coincide con la inversa
por la izquierda inversa o inversa por la derecha
Es muy interesante ver como puede calcularse la pseudoinversa de
Moore-Penroseutilizando la SVD.
Si A=U*S*V' entonces AplusMP=V*Splus*U' donde Splus es la matriz
diagonal cuyos elementos son los inversos de los de S salvo para los que
son cero.
En el ejemplo anterior:
octave:8> A=A'
A =
1 4
2 5
3 6
octave:9> [U, S, V] = svd (A)
U =
-0.42867 0.80596 0.40825
-0.56631 0.11238 -0.81650
-0.70395 -0.58120 0.40825
S =
9.50803 0.00000
0.00000 0.77287
0.00000 0.00000
V =
-0.38632 -0.92237
-0.92237 0.38632
octave:10> A*(V*[1/S(1,1),0,0;0,1/S(2,2),0]*U')*A
ans =
1.00000 4.00000
2.00000 5.00000
3.00000 6.00000
Que es justamente A.
La pseudoinversa de Moore-Penrose en MatLab y octave puede calcularla
con
octave:11> pinv(A)
ans =
-0.94444 -0.11111 0.72222
0.44444 0.11111 -0.22222
octave:12> V*[1/S(1,1),0,0;0,1/S(2,2),0]*U'
ans =
-0.94444 -0.11111 0.72222
0.44444 0.11111 -0.22222
Un saludo,
Javier
PD: Creo que esto que te escribo podrías añadirlo como artículo al Blog,
ya que es de interés para todos.