1. Ciencia

Universidad Politécnica de Madrid–Escuela Técnica Superior de Ingenieros Industriales
Matemáticas de Especialidad
Ingeniería Eléctrica
Mínimos cuadrados
Lineales
José Luis de la Fuente O’Connor
[email protected]
[email protected]
Clase_mincua_lineal_2015.pdf
1/97
Índice
Introducción
Fundamentos teóricos
Sistemas incompatibles. Ecuaciones normales
Sistemas indeterminados
Resolución numérica del problema
Método de Gram-Schmidt
Método de Gram-Schmidt modificado
Factorización QR
Descomposición numérica en valores singulares
Comparación de los métodos
Matlab y la solución de problemas de mínimos cuadrados
2/97
Introducción
El conocimiento y control del estado de funcionamiento de muchos sistemas de
operación, control, –eléctricos y mecánicos– y predicción de la vida real exigen
la toma de medidas, toma de decisiones y posteriores acciones, muchas veces en
tiempo real.
Cualquier medida siempre está sujeta a errores, por pequeños que sean, por el
desajuste de la calibración del aparato que la realiza, su propia imperfección, las
condiciones ambientales, vibraciones, envejecimiento, etc.
Par mitigar el efecto de esos errores, aislarlos, identificarlos y filtrarlos se toma
un número de medidas de los parámetros del sistema bastante mayor del
estrictamente necesario –redundante–.
3/97
La redundancia de medidas conduce normalmente, en los modelos matemáticos
que determinan los parámetros de funcionamiento de un sistema, a sistemas de
ecuaciones incompatibles: con muchas más ecuaciones que incógnitas.
La falta de suficientes medidas, a sistemas indeterminados.
Para obtener la solución más probable que represente un sistema, o
pseudosolución, y que más se aproximase a la ideal si no se diesen los errores o
perturbaciones citados, se proyecta el problema en otro de menor dimensión
–otro subespacio– para suprimir así, filtrar o aislar los datos irrelevantes o
erróneos.
La proyección más común es la que determina el método de los mínimos
cuadrados: proyección ortogonal, (ver teorema de la proyección) aunque hay
otras.
4/97
Casos posibles
sistemas
lineales, una
vez más
Casos de
posibles
de sistemas
lineales,
una vez más
m=n
m=n
·
=
m=n
m=n
rango(A) = m = n
·
=
rango(A) = m = n
·
=
m=n
·
=
·
=
·
=
·
=
rango(A) < n < m
rango(A) < n < m
·
=
rango(A) < m = n
rango(A) = m = n
rango(A) < m = n
1a
1b
1a
m>n
m>n
1b
·
=
·
=
rango(A) = n < m
rango(A) = n < m
·
=
m>n
m>n
m>n
m>n
rango(A) = n < m
rango(A) < n < m
2a
2b
2a
2b
m<n
m<n
rango(A) = m < n
m<n
rango(A) = m < n
rango(A) = m < n
3a
3a
m<n
=
m<n
·
=
rango(A) < m < n
·
=
rango(A) < m < n
·
m<n
rango(A) < m < n
3b
3b
·
=
·
=
·
=
5/101
5/97
Estudiaremos pues problemas sin solución, en los que rango.Ajb/ ¤ rango.A/,
y a los que, sin embargo, se les puede encontrar una pseudosolución que cumpla
un criterio concreto: por ejemplo, el de minimizar la norma kAx bk2.
b
a8
01234
a6
6/97
También problemas con muchas soluciones, de las que:
Se escoge aquella x cuya norma euclídea, kxk2, es mínima.
Se estudia otro tipo de solución; por ejemplo:
Que minimice
m ˇ
X
ˇ T
ˇaj x
j D1
Que minimice
nˇ
ˇ T
mKax ˇaj x
j
ˇ
ˇ
bj ˇ
ˇo
ˇ
bj ˇ
7/97
El hecho de que se utilice el criterio de minimizar, de una manera u otra, la
norma euclídea –raiz cuadrada positiva de la suma de los cuadrados de las
desviaciones entre dos vectores de valores reales– es lo que da nombre a los
procedimientos para resolver esos problemas: mínimos cuadrados.
El problema lineal de mínimos cuadrados se plantea formalmente así:
Dada una matriz A
encontrar un vector
2 Rm n, de rango k Ä mKın.m; n/, y un vector b 2 Rm,
x 2 Rn que minimice kAx bk2.
8/97
El ejemplo por excelencia de las técnicas que estudiamos lo constituye el tratar
de ajustar a un conjunto de m pares de puntos .ti ; bi / una función f .x; t / de n
parámetros independientes x1; x2 : : : xn: los coeficientes del vector x.
Los pares de puntos los pueden definir unas mediciones, bi , obtenidas en
unos tiempos, ti .
Si la función es lineal en x1; : : : ; xn se tiene un problema de mínimos
cuadrados lineales en el que:
i. Si los n parámetros se disponen como los coeficientes de un vector
n-dimensional, x ; y
ii. los datos obtenidos, en otro vector m-dimensional
m n),
se llega a una relación de la forma
b 2 Rm .
b (usualmente
Ax D b, donde A 2 Rm n, x 2 Rn y
9/97
Ejemplo
Supongamos que queremos ajustar al conjunto de pares de puntos f.ti ; bi /g =
f.1; 2/; .2; 3/; .3; 5/; .4; 6/g la función
f .x0; x1; x2; t / D x0 C x1t C x2t 2;
según representa la figura.
b
f (x 0 , x 1 , x 2 , t ) = x 0 + x 1 t + x 2 t 2
5
4
3
2
1
1
2
3
4
5
6
7
t
10/97
Para los datos y parámetros de este ejemplo el sistema Ax D b tiene la forma
siguiente:
2
3
2 3
1 2 3
2
x0
637
47
7 4x15 D 6 7 :
455
95
x2
16 ˜
6
š x
˜
b
A
1
61
6
41
1
1
2
3
4
Resolver este problema en la zona de trabajo de Matlab sería trivial:
>> Am=[1 1 1;1 2 4;1 3 9;1 4 16];
>> b=[2 3 5 6];
>> Am\b
ans =
0.50000000000000
1.40000000000000
-0.00000000000000
11/97
Con una táctica de “ajuste” por mínimos cuadrados, también con Matlab:
>> x=[1 2 3 4];
>> y=[2 3 5 6];
>> p=polyfit(x,y,2)
p =
-0.0000
1.4000
0.5000
12/97
Si se quiere profundizar un poco en el problema dibujando los datos y la función
ajustada, habría que hacer:
>> x1=linspace(0,5,150);
>> y1=polyval(p,x1);
>> plot(x,y,’o’,x1,y1,’-’)
Los resultados son los de la figura.
8
7
6
5
4
3
2
1
0
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
La línea recta y D 0; 5 C 1; 4t, en verde, y los puntos a los que se ajusta, en
azul.
13/97
Otros muchos sistemas de la ciencia, ingeniería, economía, etc. recurren a
modelos de mínimos cuadrados.
14/97
Sistemas de predicción a partir de datos masivos
15/97
Sistemas de navegación
VOR, DME, ADF, RMI, MLS
16/97
Índice
Introducción
Fundamentos teóricos
Sistemas incompatibles. Ecuaciones normales
Sistemas indeterminados
Resolución numérica del problema
Método de Gram-Schmidt
Método de Gram-Schmidt modificado
Factorización QR
Descomposición numérica en valores singulares
Comparación de los métodos
Matlab y la solución de problemas de mínimos cuadrados
17/97
Fundamentos teóricos
Teorema
Descomposición en valores singulares Si A 2 Rm
matrices ortogonales U 2 Rm m y V 2 Rn n tales que
n
es una matriz de rango r, existen
A D U †V T ;
donde
†D
† 2 Rm
n
y † r D diag. 1 ,
2; : : : ;
r /,
Ä
(1)
†r 0
;
0 0
con
1
2
r
> 0:
Si las matrices U y V se escriben como
U D Œu1 ; : : : ; um 
y V D Œv1 ; : : : ; vn  ;
los ui y vi son los vectores singulares izquierdos y derechos, respectivamente, correspondientes a los valores
singulares i , i D 1; : : : ; r.
18/97
Teorema El vector
Ä
†r 1 0
xDV
UTb
0 0
es la solución del problema
minimizar
kAx
x 2Rn
que hace mínima kxk2, donde A 2 Rm
n
bk2
y rango.A/ D r Ä mKın.m; n/.
Definición A la matriz
Ä 1
†r 0
AŽ D V
U T 2 Rn
0 0
n
se la denomina matriz pseudoinversa o inversa generalizada Moore-Penrose
de A.
19/97
Obsérvese que .A T /Ž D .A Ž/T ; en general,
.AB/Ž ¤ B ŽA Ž:
De acuerdo con esa definición, la solución de mKınx 2Rn kAx
bk2 es
x D A Žb.
b
Además, cumple que
a8
x ? ker.A/ y que Ax D P Im.A /b;
01234
a6
donde P Im.A /b es la proyección ortogonal sobre Im.A/ (paralelamente a
ker.A T /) de b.
20/97
La matriz pseudoinversa satisface las denominadas condiciones de Penrose:
AA ŽA
A ŽAA Ž
.AA Ž/T
.A ŽA/T
D
D
D
D
A
AŽ
AA Ž
A ŽA:
La pseudoinversa también proporciona unas fórmulas para la obtención de las
matrices de proyección ortogonal sobre los cuatro subespacios fundamentales de
A:
P Im.A / D AA Ž
P ker.A T / D I AA Ž
P Im.A T / D A ŽA
P ker.A /
D I
A ŽA
21/97
Dos casos de matriz pseudoinversa son de especial interés:
a) Si A 2 Rm n, m
n y rango.A/ D n,
A Ž D .A T A/ 1A T .
b) Si A 2 Rm n, m Ä n y rango.A/ D m,
A Ž D A T .AA T / 1.
a) es el problema de mínimos cuadrados resultante de un sistema de
ecuaciones incompatible, con matriz de rango completo.
b) es el de un sistema de ecuaciones compatible indeterminado con matriz
de rango completo: resuelve este problema:
minimizar kxk2;
x 2S
donde
S D fx W Ax D bg:
22/97
Sistemas incompatibles. Ecuaciones normales
Si se tiene una ecuación Ax D b, A 2 Rm n, que no tiene solución pues
b … Im.A/, se puede buscar una pseudosolución, x, que acerque Ax lo más
posible a b en el sentido de la k k2, es decir,
mKın kAx
x 2Rn
bk2:
Teorema Sean X e Y dos espacios vectoriales de dimensiones finitas n y
m sobre el cuerpo R y A una transformación lineal representada en dos bases
de X e Y por la matriz A. Para un vector dado b 2 Y , el vector x 2 X
minimiza kAx bk2 si y sólo si A T Ax D A T b.
I D EMOSTRACIÓN . Sean Im.A/ D fAx W x 2 Rn g y ker.A/ D fx W Ax D 0g.
El complemento ortogonal del conjunto Im.A/ será:
.Im.A//? D fr W r T z D 0; 8z 2 Im.A/g D fr W r T A D 0T g
D fr W A T r D 0g D ker.A T /:
O 2 , donde bO 2 Im.A/.
El problema planteado es obviamente equivalente a minimizar kb bk
Por el teorema de la proyección, bO es un vector que minimiza la norma anterior si y sólo si b bO 2
O D A T b A T Ax.
.Im.A//? ; es decir, si b bO 2 ker.A T /, o de forma equivalente, 0 D A T .b b/
23/97
La representación geométrica en tres dimensiones es esta.
7
73
2012344 3273
4
5898
01234
5696
58
3
56
El vector de residuos, r D b Ax, es ortogonal a Im.A/ y a los vectores que
lo definen: a1 y a2; es decir, A T .Ax b/ D 0.
24/97
Si la matriz A T A es invertible,
rDb
Ax D .I
P Im.A //b;
donde P Im.A / D A.A T A/ 1A T es la matriz de proyección ortogonal sobre
Im.A/ paralelamente a ker.A T /.
Al sistema de ecuaciones que define la relación
A T Ax D A T b
se le denomina
ecuaciones normales.
El vector solución x es único si A T A es invertible (si y sólo si la transformación
lineal A es inyectiva: rango.A/ D n); en este caso
x D .A T A/ 1A T b:
25/97
Sistemas indeterminados
Si el sistema tiene más de una solución, siempre se puede calcular aquella de
menor norma euclídea.
Teorema Sean X e Y dos espacios vectoriales de dimensiones n y m sobre
el cuerpo R y A una transformación lineal representada en dos bases de X
e Y por la matriz A. El vector x de norma euclídea mínima que satisface la
ecuación Ax D b es el dado por x D A T z, donde z es una solución de la
ecuación AA T z D b.
I D EMOSTRACIÓN . Si x 1 es una solución de la ecuación Ax D b, cualquier solución de la misma se
puede expresar como x D x 1 C u, donde u 2 ker.A/; es decir, estará en la variedad lineal x 1 C ker.A/.
El teorema de la proyección garantiza la existencia en esta variedad lineal de un único x tal que su
norma kxk2 es mínima y además pertenece a .ker.A//? .
Como x 2 .ker.A//? , pertenecerá a Im.A T /, es decir, se podrá expresar como x D A T z para algún
z 2 Y . Como Ax D b, entonces
AA T z D b:
26/97
Cuando la matriz AA T es invertible, la solución óptima es
x D A T .AA T / 1b:
La interpretación geométrica de este resultado en R3 se esquematiza así:
69
678012345
67
6
012345
27/97
Ejemplo
Ä
x1
D 3, con norma euclídea mínima. El
x2
sistema Ax D b es evidentemente indeterminado.
Se quiere resolver el problema 1 2
Cualquier solución se podrá expresar como x 1 C ker.A/, donde x 1 es cualquier
T
75
en la figura.
vector solución (por ejemplo 1 1 ) y ker.A/ es el que se ve
SOLUTION FOR A SYSTEM OF LINEAR EQUATIONS
2
Im.A T /
x2
1.5
(3/5, 6/5) = x
(1, 1) = x 1
1
subespacio de soluciones
(3, 0)
1
2
x1
3
ker.A/
La solución que se busca es Figure 2.1
A minimum-norm solution.
28/97
La solución que se busca es
x DA T AA T
Ä
Ä
Ä Â
Ä Ã 1
3
1
0;6
1
1
1
D
:
3D
bD
1 2
1;2
2
2
5 2
Se ve que está en Im.A T /.
29/97
Índice
Introducción
Fundamentos teóricos
Sistemas incompatibles. Ecuaciones normales
Sistemas indeterminados
Resolución numérica del problema
Método de Gram-Schmidt
Método de Gram-Schmidt modificado
Factorización QR
Descomposición numérica en valores singulares
Comparación de los métodos
Matlab y la solución de problemas de mínimos cuadrados
30/97
Resolución numérica del problema
Las ecuaciones normales se pueden resolver con cualquier método que sea
adecuado para sistemas en los que la matriz es cuadrada y simétrica,
aplicándoselo a
A TAx D A T b, en el caso de que el sistema fuese incompatible, o a
AA Tz D b, cuando se diese un sistema indeterminado.
Los números de condición, Ä2, de AA T y A T A, son el cuadrado del de la
matriz A, por lo que si el problema originalmente no está bien condicionado, las
dificultades numéricas pueden resultar insalvables al resolver el sistema
correspondiente.
31/97
Ejemplo
Consideremos la matriz
2
3
1 1 1 1 1
6"
7
6
7
6
7
6 "
7
AD6
7:
"
6
7
6
7
4
" 5
"
El rango de A es 5, para
2 " ¤ 0. La matriz
3
2
1C"
1
1
1
1
6 1
1 C "2
1
1
1 7
6
7
6
7
T
A AD6 1
1
1 C "2
1
1 7
6
7
1
1
1 C "2
1 5
4 1
1
1
1
1
1 C "2
también es de rango 5, para " ¤ 0.
32/97
El número de condición Ä2.A T A/ D Ä2.A/2 D .5 C "2/="2.
Si " es mayor que la precisión de la máquina pero "2 no (por ejemplo, si
" D 0,5 10 5, "2 D 0,25 10 10 y la precisión de la máquina
D 1,0 10 10), la representación interna de la matriz A T A será
2
3
1 1 1 1 1
61 1 1 1 17
6
7
6
7
61 1 1 1 17
6
7
41 1 1 1 15
1 1 1 1 1
por lo que, a efectos numéricos en esa máquina, esta matriz será singular y de
rango 1: las ecuaciones normales no servirían.
33/97
Otro aspecto importante que aconseja tener mucho cuidado al utilizar A T A ó
AA T , nace del hecho de que aun cuando la matriz original A tenga muchos
elementos cero, A T A o AA T pueden ser totalmente densas.
Un ejemplo sería
2
3
1 1 1 1
60
7
6
7
6
7
AD6 0
7;
6
7
0 5
4
0
y
2
1
6
61
AT A D 6
41
1
1
1
1
1
1
1
1
1
3
1
7
17
7:
15
1
34/97
Método de Gram-Schmidt
El procedimiento clásico de Gram-Schmidt obtiene una base ortonormalizada
del subespacio Im.A/.
Comienza normalizando el primer vector columna de A: e 1 D a1=ka1k2.
A continuación se sustrae del vector a2 su coeficiente en la dirección de e 1,
ha2je 1ie 1, resultando un vector ortogonal a e 1, el cual a su vez se
normaliza: : :
El proceso continúa con los demás vectores columna de A.
El número de operaciones del método es O.mn2/ sumas+restas y
multiplicaciones+divisiones y O.n/ raíces cuadradas.
35/97
Los diversos vectores ortonormales de la base se obtenien así:
a1
e1 D
I
ka1k2
a2 ha2je 1ie 1
e2 D
I
ka2 ha2je 1ie 1k2
a3 ha3je 1ie 1 ha3je 2ie 2
I
e3 D
ka3 ha3je 1ie 1 ha3je 2ie 2k2
:::
a
e3
−
3
a3
e2
a
e1
3|
e1
a3 − a3 |e1 e1 − a3 |e2 e2
a3 |e2 e2
e1
a3 |e1 e1
36/97
El algoritmo, para una matriz general A m n, es el siguiente.
Ortogonalización de A por Gram-Schmidt
for j D 1 to n
e.1 W m; j /
a.1 W m; j /
for i D 1 to j 1
u.i; j /
e.1 W m; i /T a.1 W m; j /
e.1 W m; j /
e.1 W m; j / u.i; j / e.1 W m; i /
end
v
uX
u m
t
u.j; j /
e.k; j /2
kD1
e.1 W m; j /
end
e.1 W m; j /=u.j; j /
37/97
El algoritmo hace A D E U , donde E m n es la matriz de columnas e i y U n
la matriz triangular superior de los productos interiores auxiliares uij .
n
Sustituyendo esta expresión de A en las ecuaciones normales,
A T Ax D A T b, resulta que
U T ET EU x D U T ET b
y, por fin, dado que E T E D I,
Sistema triangular superior.
U x D E T b:
En condiciones adecuadas, por consiguiente, el método de Gram-Schmidt podría
valer para resolver un problema lineal de mínimos cuadrados.
38/97
Método de Gram-Schmidt modificado
En la práctica, la perdida de ortogonalidad de los vectores e i que se van
calculando en el proceso de Gram-Schmidt es más que evidente, debido a errores
de redondeo y de cancelación si, por ejemplo, alguno de los vectores columna aj
está próximo al subespacio generado por los vectores anteriores e 1; : : : ; ej 1.
Pj 1
En ese caso, los sumandos de la expresión aj
i D1 haj je i ie i pueden llegar a
ser muy pequeños o muy distantes unos de otros, si bien el resultado final puede
ser muy pequeño, por lo que el error numérico que se va produciendo es
relativamente grande. Al dividir el resultado por su norma (también muy
pequeña) los errores se amplificarán aún más.
39/97
En 1966 Rice modificó el método haciendo que en una etapa k, en vez de
sustraer del vector ak sus coeficientes sobre todos los k 1 vectores e i ya
calculados, e k se hace igual a ak al principio y luego se le van sustrayendo su
proyección en e 1, pasando a ser el nuevo e k el resultado, el cual se proyecta
luego en e 2, y así sucesivamente en cada uno de los k 1 e i anteriores. El
resultado con esta simple nueva disposición de los cálculos es indiscutiblemente
mejor numéricamente.
Otras versiones, una vez calculado un vector ortonormal e k en una etapa k, van
sustrayendo de cada ai , i D k C 1; : : : ; n, su proyección sobre e k . El resultado
es el mismo.
40/97
Los dos algoritmos se listan a continuación.
for j D 1 to n
e.1 W m; j /
a.1 W m; j /
1
u.i; j /
e.1 W m; i /T a.1 W m; j /
e.1 W m; j /
e.1 W m; j / u.i; j / e.1 W m; i /
for i D 1 to j
Clásico
end
u.j; j /
v
uX
u m
t
e.k; j /2
kD1
e.1 W m; j /
e.1 W m; j /=u.j; j /
end
for j D 1 to n
e.1 W m; j /
a.1 W m; j /
1
u.i; j /
e.1 W m; i /T e.1 W m; j /
e.1 W m; j /
e.1 W m; j / u.i; j / e.1 W m; i /
for i D 1 to j
Modificado
end
u.j; j /
v
uX
u m
t
e.k; j /2
kD1
e.1 W m; j /
end
e.1 W m; j /=u.j; j /
41/97
La versión clásica y modificada en Matlab son estas.
function [x r2 e]=Grmsch_3(A,b)
% Se resuelve Ax=b mediante el método de Gram-Schmidt modificado
[m,n]=size(A); x=zeros(n,1); e=zeros(m,n); u=triu(zeros(n,n));
for j=1:n
e(:,j)=A(:,j);
for i=1:j-1
function [x r2 e]=Grmsch_2(A,b)
% Se resuelve Ax=b mediante el método de Gram-Schmidt clásico
u(i,j)=e(:,i)’*e(:,j);
[m,n]=size(A); x=zeros(n,1); e=zeros(m,n); u=triu(zeros(n,n));
e(:,j)=e(:,j)-u(i,j)*e(:,i);
for j=1:n
e(:,j)=A(:,j);
end
for i=1:j-1
u(j,j)=norm(e(:,j));
u(i,j)=e(:,i)’*A(:,j);
e(:,j)=e(:,j)-u(i,j)*e(:,i);
e(:,j)=e(:,j)/u(j,j);
end
end
u(j,j)=norm(e(:,j));
for i=n:-1:1
% Rx=b
e(:,j)=e(:,j)/u(j,j);
x(i)=(e(:,i)’*b-u(i,i+1:n)*x(i+1:n))/u(i,i); end
for i=n:-1:1
% Rx=b
end
x(i)=(e(:,i)’*b-u(i,i+1:n)*x(i+1:n))/u(i,i);
r2=norm(abs(A*x-b),2)^2;
%end
Residuos
r2=norm(abs(A*x-b),2)^2;
% Residuos^2
end
end
42/97
El cara a cara, en Matlab supercompacto, del clásico y el modificado en la
versión alternativa es este.
function [Q, R] = gs_m_VS_c(A)
[m, n] = size(A);
Q = zeros(m,n);
R = zeros(n);
[m, n] = size(A);
Q = zeros(m,n);
R = zeros(n);
for k=1:n
R(k,k) = norm(A(:,k));
Q(:,k) = A(:,k)/R(k,k);
R(k,k+1:n) = Q(:,k)’*A(:,k+1:n);
A(:,k+1:n) = A(:,k+1:n) - Q(:,k)*R(k,k+1:n);
end
for j=1:n
R(1:j-1,j) = Q(:,1:j-1)’*A(:,j);
temp = A(:,j) - Q(:,1:j-1)*R(1:j-1,j);
R(j,j) = norm(temp);
Q(:,j) = temp/R(j,j);
end
43/97
Si resolvemos el ejemplo inicial.
>> A=[1 1 1;1 2 4;1 3 9;1 4 16]
A =
1
1
1
1
2
4
1
3
9
1
4
16
>> b=[2;3;5;6]
b =
2
3
5
6
>> [x r2]=Grmsch_3(A,b)
x =
0.5000
1.4000
-0.0000
r2 =
0.2000
>> A\b
ans =
0.5000
1.4000
-0.0000
44/97
% Script_GRSCH_1.m - Script de Ortogonalidad con Gram Schmidt clásico y modificado
format short
n=7; A=hilb(n);
% Matriz de Hilbert (muy mal condicionada)
cond(A), pause
b=A*ones(n,1);
% Término independiente para sol. x=1.
disp(’Clásico:’), [x r2 e]=Grmsch_2(A,b);
% Gram Schmidt clásico
x, pause, e, pause, norm(abs(x-ones(n,1)),2), pause
ortogonalidad=norm(e’*e-eye(n)), pause
% Ortogonalidad real de la matriz e
disp(’Modificado:’), [x r2 e]=Grmsch_3(A,b);% Gram Schmidt modificado
x, pause, e, pause, norm(abs(x-ones(n,1)),2), pause
ortogonalidad=norm(e’*e-eye(n))
% Ortogonalidad real de la matriz e
45/97
Factorización QR
Recordemos que una propiedad importante de las transformaciones ortogonales
es que conservan la norma euclídea; esto es, si Qn n es una matriz ortogonal y
x un vector n-dimensional, se cumple que
kQxk2 D kxk2:
En efecto:
kQxk2 D
p
hQxjQxi D
q
T
x T Q Qx D
p
x T x D kxk2:
Según esto, si Q es una matriz ortogonal, al premultiplicar el vector Ax
por ella, su norma euclídea queda igual:
kQAx
Qbk2 D kQ.Ax
b/k2 D kAx
b
bk2:
46/97
Unas transformaciones ortogonales adecuadas pueden convertir el problema de
mínimos cuadrados en otro más sencillo de resolver numéricamente.
Si A 2 Rm n, m > n, b 2 Rm, rango.A/ D n y se ha calculado una matriz
ortogonal Q 2 Rm m tal que
Ä
R1
n
QA D R D
0 m n
tiene R 1 triangular superior, si se hace
Ä
c
Qb D
d
n
m
entonces
kAx
bk2 D kQAx
q
kR 1x
D
Qbk2 D
n
Ä
;
R 1x c
d
2
ck22 C kdk22; para cualquier x 2 Rn.
47/97
La solución de mKınx 2Rn kAx
bk2 será aquella que haga mínimo
kR 1x
ck22 C kdk22:
Como kdk22 es constante, la solución será la que haga mínimo el otro
sumando: cuando R 1x D c .
Resolviendo este sistema por sustitución inversa se llega a la solución del
problema de mínimos cuadrados.
La suma de residuos al cuadrado será kdk22 y el vector de residuos
Ä
0
r D QT
:
d
El proceso de reducción de A a R se denomina
triangularización ortogonal.
factorización QR o
48/97
Teorema Sea la matriz A 2 Rm n de rango n y su factorización A D QR.
El factor R tiene todos los elementos de su diagonal principal positivos y es
igual al que resulta de la factorización de Cholesky, G T G , de A T A.
I D EMOSTRACIÓN . Si rango.A/ D n, de acuerdo con un teorema anterior, la factorización de Cholesky
de A T A es única.
El teorema de la proyección garantiza la existencia en esta variedad lineal de un único x tal que su
norma kxk2 es mínima y además pertenece a .ker.A//? .
Por otro lado,
Ä
h
i
R1
A T A D R T1 ; 0 QQT
D R T1 R 1 :
0
49/97
Transformaciones de Householder
Definición Se denomina transformación o reflexión de Householder a una
transformación lineal de Rn en Rn caracterizada por una matriz H n n –matriz
de Householder– de la forma
H DI
2wwT ;
donde w 2 Rn; kwk2 D 1, es el vector de Householder.
Teorema Toda transformación de Householder es simétrica y ortogonal.
I D EMOSTRACIÓN . Por definición
HT D I
2.wwT /T D I
2.wT /T wT D I
2wwT D H :
Como además wT w D kwk22 D 1,
HTH D H2
D .I 2wwT /.I 2wwT /
D I 4wwT C 4w.wT w/wT D I:
50/97
Aplicar una transformación de Householder a un vector cualquiera equivale a
“reflejarlo” en el subespacio .Im.w//?.
6
3
01203445
8 7 396 3
8 9
6 7 3 6 3
En efecto
H a D .I
2wwT /a D a
2wwT a D a
2.wT a/w:
El vector .wT a/w es la proyección de a sobre w; es decir, H a es igual al
vector a menos dos veces su proyección sobre w.
51/97
La importancia fundamental de estas transformaciones radica en su capacidad
para modificar ortogonalmente –eventualmente hacer cero– determinados
coeficientes de un vector dado.
En efecto, si x e y son dos vectores no nulos de igual norma euclídea y w
se define como
wD
entonces
.I
1
kx
yk2
.x
y/;
2wwT /x D y:
52/97
Comprobémoslo:
Â
Ã
.x y x
x y
p
2 p
D
.x y /T .x y /
.x y /T .x y /
Á
xT x y T x
= x 2 .x y /T .x y / .x y/
Á
xT x y T x
= x 2 2.x T x y T x / .x y/ D y.
x
/T
Esto es así pues, al tener x e y la misma norma euclídea,
.x
y/T .x
y/ D x T x y T x x T y C y T y
D 2.x T x y T x/;
pues x T x D y T y y y T x D x T y.
53/97
Este resultado, geométricamente, se deduce inmediatamente de la reflexión
antes mencionada.
El vector w es colineal con el vector x
y.
Como x e y tienen la misma longitud, la reflexión de x respecto a
.Im.w//? es y.
6
678
3
01203445
8
54/97
En un sistema Ax D b se pueden construir transformaciones de Householder
que anulen los coeficientes que se deseen de cada vector columna de A dejando
los demás como estaban.
La figura representa los cuatro pasos del proceso de reducir una matriz A 6
a una triangular superior R 6 4.
××××
××××
××××
××××
××××
××××
A0
✷✷✷✷
0 ✷✷✷
0 ✷✷✷
0 ✷✷✷
0 ✷✷✷
0 ✷✷✷
A1
✷✷✷✷
0
0 0
0 0
0 0
0 0
A2
✷✷✷✷
0
0 0
0 0 0
0 0 0
0 0 0
A3
4
✷✷✷✷
0
0 0
0 0 0✸
0 0 0 0
0 0 0 0
A4
La matriz A 1 resultaría de la transformación H 1A 0; la A 2 sería
H 2A 1 D H 2H 1A 0; y así cuatro veces.
55/97
A una matriz m n se aplicará una sucesión k D 1; : : : ; n de transformaciones
que hagan cero los coeficientes k C 1; : : : ; m del vector columna k.
Se pretende que la transformación k-ésima haga:
H k ak D y D
aik para i D 1; 2; : : : ; k 1
0 para i D k C 1; : : : ; m:
Los coeficientes del vector y deberán ser:
y1
y2
yk
ykC1
ym
D a1k
D
::: a2k
q
2
2
C akC1k
C
D ˙ akk
D
::: 0
D 0:
Como y12 C y22 C
2
2
C ym2 D a1k
C a2k
C
2
C amk
¡OJO!
2
C amk
, jjyjj2 D jjak jj2.
La transformación de Householder que hay que construir tendrá como vector
w D .ak y/=kak yk2.
56/97
Como se observa, el coeficiente k-ésimo del vector y puede adoptar dos signos.
Para evitar errores, se escoge el signo de la raíz cuadrada opuesto al del
coeficiente k-ésimo de ak .
x2
w
w alt
y
a1 − y
a
a1
y
a1 + y
x1
En esta figura hay dos posibilidades, y o y 0, de transformar a en uno de ellos
con a2 D 0 y la misma magnitud:
Si la transformación de Householder la define w, que convierte a en y 0, su
primer coeficiente es w1 D a1 . y1/ D a1 C y1; el segundo, el de a en x2.
La transformación alternativa, walt convierte a a y: su primer coeficiente es
a1 y1; el segundo el mismo de w.
57/97
En definitiva, en la transformación de Householder k-ésima que se aplica a la
matriz A, los valores numéricos del vector w son:
2
3
0
::
6
7
6
7
6
1
akk C s signo.akk /7
6
7;
wDp
6
7
akC1k
2s.s C jakk j/ 6
7
:
4
5
:
amk
q
2
2
C akC1k
C
donde s D akk
2
C amk
.
58/97
Caso numérico 1: Resolución de Ax D b, A m n, m > n y
rango completo
Mediante transformaciones
de Householder seÄreduce la matriz A a una
Ä
R1
c
triangular superior
y el vector b a otro
.
0
d
La solución de mKınx 2Rn kAx
sustitución inversa).
bk2 sería la del sistema R 1x D c (por
La suma de residuos al cuadrado será kdk22.
El algoritmo numérico completo es el que describe la tabla que sigue.
59/97
Resolución de Ax D b por transf. de Householder
Transformación columnas de la Matriz A m
for j D 1 to n
if mKa0
x fja.j; j /j; : : : ; ja.m;
j /jg D 0 then stop
1
v
uX
m
u
D @t
a.k; j /2 A signo.a.j; j //
n
kDj
w.j W m/
a.j W m; j /; w.j /
w.j / C ; ˇ D 2
m
X
w 2 .k/; a.j; j /
kDj
for l D j C 1 to n
a.j W m; l/
a.j W m; l/ w.j W m/ w T .j W m/ a.j W m; l/
end
Transformación del vector b.
b.j W m/
b.j W m/ w.j W m/ w T .j W m/ b.j W m/ ˇ
end
Resolución del sistema Rx D b.
for j D n to01
1,
n
X
@b.j /
x.j /
a.j; k/ x.k/A
a.j; j /
ˇ
kDj C1
end
Residuos al cuadrado.
m
X
rescua
b 2 .k/
kDnC1
60/97
En Matlab, con la posibilidad de obtener Q y R:
function [x r2 Q R]=Qrdes_3(A,b)
% Resolución de Ax=b mediante transformaciones de Householder; calcula Q y R
[m n]=size(A); x=zeros(n,1); Q=eye(m);
for j=1:n
w=Housv(A(j:m,j));
% Householder de a(j:m,j)
A(j:m,j:n)=A(j:m,j:n)-2*w*(w’*A(j:m,j:n));
b(j:m)=b(j:m)-2*w*(w’*b(j:m));
Qk=eye(m);
Qk(j:m,j:m)=eye(m+1-j)-2*(w*w’);
Q=Qk*Q;
end
for i=n:-1:1
% Rx=b
x(i)=(b(i)-x(i+1:n)’*A(i,i+1:n)’)/A(i,i);
end
r2=norm(b(n+1:m))^2;
% Residuos al cuadrado
R=triu(A); Q=Q’;
% Matrices R y Q
function w = Housv(x)
% Transformación de Householder del vector x.
m=max(abs(x)); w=x/m;
sw=1; if w(1)<0, sw=-1; end
w(1)=w(1)+sw*norm(w);
w=w/norm(w); w=w(:);
61/97
Si utilizamos este script para resolver el ejemplo.
>> A=[1 1 1;1 2 4;1 3 9;1 4 16]
A =
1
1
1
1
2
4
1
3
9
1
4
16
>> b=[2;3;5;6]
b =
2
3
5
6
>> [x r]=Qrdes_3(A,b)
x =
0.5000
1.4000
-0.0000
r =
0.2000
>> A\b
ans =
0.5000
1.4000
-0.0000
62/97
El vector de residuos,
2
6
6
rD6
4
3
0,1
7
0,37
7;
0,35
0,1
Ä
c
se puede obtener, si el algoritmo ha transformado
b
en
, sin más que hacer:
Ä
d
0
r
d
for k D n to 1
r
H kr
end
El número de operaciones de este método es:
O.mn2 n3=3/ sumas+restas y multiplicaciones+divisiones, para
transformar la matriz A en R;
n raíces cuadradas y las de la sustitución inversa, O.n2=2/ .
63/97
Caso numérico 2: Resolución de Ax D b, A m n, n > m y
rango completo
Este problema, si tiene solución, es indeterminado: tiene muchas.
La de menor norma euclídea se puede calcular mediante estos pasos:
Paso 1 Se aplica el algoritmo QR a la matriz A T , en vez de a A. Resultará
Ä
R
QT A T D
;
0
Ä
R
o, lo que es lo mismo, A T D Q
, donde Q es una matriz ortogonal
0
n n y R una triangular superior m m.
64/97
Paso 2 La matriz original A será
A D R T ; 0T QT :
Si se sustituye en la ecuación Ax D b, se tendrá que
R T ; 0T QT x D b:
Si se hace el cambio de variable z D QT x, la última ecuación queda
R T ; 0T z D b:
Como zT z D .QT x/T .QT x/ D x T QT Qx D x T x, las normas
euclídeas de x y z serán iguales.
Ä
z
Estructurando el vector z en R y b de igual manera, la solución de
z0
R T ; 0T z D b saldrá de resolver
R T zR D bR ;
siendo los demás coeficientes del vector z, z0, nulos.
65/97
Paso 3 El vector solución x que se busca resultará de deshacer el cambio de
variable introducido; es decir:
Ä
z
xDQ R :
0
66/97
Caso 3: Resolución numérica de Ax D b, A m n, m > n ó
m < n y rango incompleto
Caso más general que se puede dar en mínimos cuadrados de ecuaciones lineales.
Paso 1 Se transforma la matriz A mediante transformaciones de
Householder y permutaciones de columnas para llegar a:
r
0
m−r
En cada etapa k se calcula la norma euclídea de ak , k; : : : ; n,
limitándólos a sus coeficiente k; : : : ; m. Se intercambia la columna k con
la de mayor norma.
67/97
Paso 2 Se ha llegado a Ä
Ä
R 11 R 12
r
c
y Qb D
0 R 22
m r
d
r n r
donde kR 22k2 Ä 1kAk2. A partir de aquí hay dos opciones:
QAP D R D
r
m
r
;
Que r D n (rango completo). La solución sale de resolver R 11x D c.
Que r < n (rango incompleto). Se deberán construir unas
transformaciones ortogonales, Qn1 n, tales que
R 11; R 12 Q1 D W ; 0 , donde W r r es triangular superior.
¿Cómo hacerlo? Se actúa sobre ŒR 11, R 12T , a) en la figura, y se llega a b).
a)
b)
r
0
n−r
r
Esto se hace en r etapas. En una de ellas, k, se premultiplica por una
transformación de Householder que haga cero los elementos r C 1 a n de
la columna k y que deje inalterados del 1 al k 1 y del k C 1 a r.
×
×
×
×
×
×
×
×
×
×
×
×
× ×
× ⊗
× ⊗
1
×
×
×
×
×
×
×
×
×
×
×
×
× ×
⊗ 0
⊗ 0
2
×
×
×
×
×
×
×
×
×
⊗
⊗
×
× ×
0 0
0 0
3
×
×
×
×
⊗
⊗
×
×
×
0
0
×
× ×
0 0
0 0
4
×
×
×
×
0
0
×
×
×
0
0
×
× ×
0 0
0 0
El óvalo indica el elemento que se utiliza para definir cada
transformación; los que se hacen cero con el signo ˝.
69/97
Paso 3 De los dos pasos anteriores se tendrá que
kAx
bk2 D k.QAP/.P T x/
Qbk2:
Ahora bien, .QAP/P T x se puede escribir .QAP/Q1QT1 P T x y
también,
Ä
Ä
W 0
c
QT1 P T x D
:
0 0
d
Si se hace QT1 P T x D y y se resuelve W y 1 D c, el vector solución que
se busca, x, resultará de
Ä
y1
x D PQ1
:
0
70/97
function [x r res]=Mincua_QR(A,b)
% Resolución de Ax=b general mediante transformaciones de Householder
%
Posible rango incompleto r
[m n]=size(A); x=zeros(n,1); tol=sqrt(eps); W=zeros(n,m); ip=1:n; r=n;
for j=1:n
jm=j; c=0;
for k=j:n
h=norm(A(j:m,k));
if h>c, c=h; jm=k; end
end
if jm~=j, A(:,[j jm])=A(:,[jm j]); ip([j jm]) = ip([jm j]); end
if j==m, break, end
w=Housv(A(j:m,j));
% Householder de A(j:m,j); luego a A y b
A(j:m,j:n)=A(j:m,j:n)-2*w*(w’*A(j:m,j:n)); b(j:m)=b(j:m)-2*w*(w’*b(j:m));
end
for j=1:n
% Ver rango
if abs(A(j,j))<=tol, r=j-1; break, end
end
res=norm(b(r+1:m))^2; w1=zeros(r,n-r+1);
W(1:n,1:r)=A(1:r,1:n)’;
% Trasp. de A = W
if r~=n
for i=r:-1:1
w1(i,1:n-r+1)=Housv([W(i,i);W(r+1:n,i)]);
% Householder hacia W
W([i r+1:n],i:-1:1)=W([i r+1:n],i:-1:1)-2*w1(i,:)’*(w1(i,:)*W([i r+1:n],i:-1:1));
end
end
for i=r:-1:1
% Resol. Wx=c
x(i)=(b(i)-x(i+1:r)’*W(i+1:r,i))/W(i,i);
end
if r~=n
x(r+1:n)=0;
% Aplicar a x ultimas Householder hacia W
for i=1:r, x([i r+1:n])=x([i r+1:n])-2*w1(i,:)’*(w1(i,:)*x([i r+1:n])); end
end
x(ip)=x(:);
% Deshacer perm. col.
71/97
Apliquemos este programa en una sesión de Matlab:
» a=rand(200,4);
» b=a*ones(4,1);
» A=[a a(:,1)*2+a(:,2)*0.5 a(:,3)*2+a(:,4)*0.5 a(:,2)*2+a(:,3)*0.5];
» size(A)
ans = 200
7
» format long
» [x r res]=Mincua_QR(A,b)
x =0.168704156479218
0.156479217603912
0.009779951100245
0.792176039119804
0.415647921760391
0.415647921760391
0.317848410757946
r = 4
res = 1.205973193713402e-029
>> x-pinv(A)*b
ans = 1.0e-015 *
-0.111022302462516
0.333066907387547
0.194289029309402
-0.333066907387547
0.166533453693773
-0.943689570931383
0.499600361081320
72/97
Si lo utilizamos para resolver el ejemplo.
>> A=[1 1 1;1 2 4;1 3 9;1 4 16]
A =
1
1
1
1
2
4
1
3
9
1
4
16
>> b=[2;3;5;6]
b =
2
3
5
6
>> [x r res]=Mincua_QR(A,b)
x =
0.500000000000000
1.400000000000000
0
r =
3
res =
0.200000000000001
>> A\b
ans =
0.5000
1.4000
-0.0000
73/97
Transformaciones de Givens
Definición Se denomina transformación de Givens a una transformación
lineal ortogonal de Rn en Rn caracterizada por una matriz de Givens
2
3
1
6 :::
7
6
7
1
6
7
6
7
c
s
i
6
7
:
7
::
G .i; j / D 6
6
7
j
6
7
s
c
6
7
1
6
7
4
::: 5
1
donde c 2 C s2
D 1.
74/97
Al aplicar a x 2 Rn una transformación, o rotación, de Givens,
G .i; j / W Rn ! Rn, con c D cos  y s D sen  , producirá lo siguiente:
2
x:1
::
6
6
xi 1
6
6 xi cos  C xj sen Â
6
xi:C1
6
6
::
G .i; j /x D 6
xj 1
6
6 x sen  C x cos Â
j
6 i
6
xj:C1
4
::
xn
3
7
7
7
7
7
7
7
7
7
7
7
7
5
i
j
Se rota el vector x un ángulo  en el subespacio que generan los vectores e i y
ej de Rn.
75/97
Si se desea hacer cero alguno de los coeficientes i ó j de un vector x,
concretamente el j , se deberá escoger un  tal que xi sen  C xj cos  D 0,
es decir, habrá que hacer
xj
tan  D ;
xi
o, lo que es equivalente,
c D cos  D q x2 i 2
xi Cxj
y
s D sen  D
x
q j
:
xi2 Cxj2
76/97
Ejemplo
En la figura se describe, en el espacio euclídeo tridimensional, la rotación del
vector
2 3
1
x D 415
1
en el plano z
y para anular el tercer coeficiente.
z
....
....
....
....
....
....
...
...
...
...
..
..
...
...
...
..
..
..
..
..
..
..
..
..
..
..
..
...
...
...
...
...
..
..
..
x
1
1
............................
.....
...
..
...
...
..
.
...
..
..
..
.
.
.
...
.
.
.
.
.
.
.
.
.
....
x
y
1
x
77/97
Como el ángulo que hay que rotar x es 45ı, la matriz de Givens que hay que
utilizar es,
2
3
1 p0
0
p
6
7
G .2; 3/ D 40 p2=2 p2=25 :
2=2 2=2
0
El nuevo vector será
2
3
1
p
G x D x 0 D 4 25 :
0
p
La norma euclídea de éste y del original es 3.
78/97
Las transformaciones de Givens se pueden utilizar en un problema de mínimos
cuadrados para transformar la matriz A m n, en n etapas, en una triangular
superior R. En cada una de esas etapas, j , se harían cero, uno a uno, los
coeficientes j C 1 a m.
Por ejemplo, las operaciones
necesarias
para transformar la matriz
2
3
2
AD4
3
6 0
7
RD4 0 0
5, son las que siguen.
0 0 0
3
2✷ ✷ ✷3
2
0 ✷ ✷
0 ✷ ✷5
A 1 D G .1; 2/A D 4
; A 2 D G .1; 3/A 1 D 4 0 ✷ ✷ 5 ;
5
en
2
5
60
A 3 D G .1; 4/A 2 D 4 0
0
2
5
60
A 5 D G .2; 4/A 4 D 4 0
0
5 5
✷
✷
✷
5
5
0
0
3
2
5
✷7
60
✷ 5 ; A 4 D G .2; 3/A 3 D 4 0
✷
0
3
2
5
5
57
60
5 ; A 6 D G .3; 4/A 5 D 4 0
0
5 5
0
✷
5
5
0
0
✷
3
7
5;
3
5
57
55 :
0
Los símbolos ✷, y 5 indican que el coeficiente al que se hace referencia ha experimentado 1, 2 ó 3 transformaciones
desde su valor inicial .
79/97
Resolución de Ax D b por transf. de Givens
Transformación de la Matriz A m n
for i D 1 to n
for k D i C 1 to m
Hacer nulo el elemento .k; i /.
if a.k; i / ¤ 0 then
if ja.k; i /j ja.i; i /j then p
t D a.i; i /=a.k; i /I s D 1= 1 C t 2 I c D s t
else
p
t D a.k; i /=a.i; i /I c D 1= 1 C t 2 I s D c t
end
a.i; i /
c a.i; i / C s a.k; i /
for j D i C 1 to n
aux D c a.i; j / C s a.k; j /; a.k; j /
s a.i; j / C c a.k; j /; a.i; j /
end
Transformación del vector b.
aux D c b.i / C s b.k/; b.k/
s b.i / C c b.k/; a.i /
aux
end
end
end
Resolución del sistema Rx D b.
for j D n to01
1,
n
X
@b.j /
x.j /
a.j; k/ x.k/A
a.j; j /
aux
kDj C1
end
Residuos al cuadrado.
m
X
rescua
b 2 .k/
kDnC1
80/97
En Matlab:
function [x r2]=Givens(A,b)
% Resolución de Ax=b mediante transformaciones de Givens
[m,n]=size(A); x=zeros(n,1);
for i=1:n
% Factorización de A
for k=i+1:m
if 1+abs(A(k,i))==1, continue, end
if abs(A(k,i))>=abs(A(i,i))
t=A(i,i)/A(k,i); s=1/sqrt(1+t*t);
c=s*t;
else
t=A(k,i)/A(i,i); c=1/sqrt(1+t*t);
s=c*t;
end
A(i,i)=c*A(i,i)+s*A(k,i);
q(i+1:n)=c*A(i,i+1:n)+s*A(k,i+1:n);
A(k,i+1:n)=-s*A(i,i+1:n)+c*A(k,i+1:n);
A(i,i+1:n)=q(i+1:n);
q1=c*b(i)+s*b(k);
% Transformar b
b(k)=-s*b(i)+c*b(k);
b(i)=q1;
end
end
for i=n:-1:1
% Sustitución inversa
x(i)=(b(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);
end
r2=norm(b(n+1:m))^2;
% Residuos al cuadrado
81/97
Si utilizamos este nuevo script para resolver el ejemplo.
>> A=[1 1 1;1 2 4;1 3 9;1 4 16]
A =
1
1
1
1
2
4
1
3
9
1
4
16
>> b=[2;3;5;6]
b =
2
3
5
6
>> [x,r]=Givens(A,b)
x =
0.5000
1.4000
-0.0000
r =
0.2000
>> A\b
ans =
0.5000
1.4000
-0.0000
82/97
El número de operaciones que requiere este algoritmo para transformar A es
O.2mn2 2n3=3/ sumas+restas y multiplicaciones+divisiones, y O.mn=2/
raíces cuadradas; además, O.n2=2/ sumas+restas y multiplicaciones+divisiones
para efectuar la sustitución inversa.
Con precisión semejante, Givens es el doble de caro que Householder. ¿Cuándo
utilizar Givens y Householder entonces?
La clave está en considerar la estructura de la matriz A del problema: si
ésta es densa, es decir, muchos de sus coeficientes son distintos de cero, el
método de Householder es el más aconsejable; si la estructura de A es
dispersa, convendrá centrarse en hacer cero sólo aquellos elementos no nulos
en las columnas correspondientes, por lo que, a priori, el método de Givens
deberá ser más ventajoso.
83/97
Transformaciones rápidas e Givens
Obsérvese que al aplicar a A una transformación de Givens, definida por una
matriz G .i; j /, los únicos elementos que se ven modificados son los de las filas
i y j , de tal forma que
aik
aj k
c aik C s aj k
s aik C c aj k ;
(2)
para k D 1; : : : ; n. Para cada k, se hacen cuatro multiplicaciones (4n en total).
Si alguno de los c o s se pudiese reemplazar por un 1, el número de
multiplicaciones se reduciría a la mitad.
Esto es lo que llevan a efecto las transformaciones rápidas de Givens mediante
una reordenación de los cálculos.
84/97
Descomposición en valores singulares
Existen diversos métodos –iterativos– muy sofisticados para calcular la
descomposición numérica en valores singulares
A D U †V T :
Según el teorema de la descomposición en valores singulares, la solución de
mKınx 2Rn kAx bk2 es
Ä 1
†r 0
xDV
U T b:
0 0
Para una matriz A de cualquier forma y rango, la solución del problema de
mínimos cuadrados de menor norma euclídea es
X uT b
i
vi :
xD
i ¤0
i
85/97
Ejemplo
Se desea resolver, utilizando la descomposición U †V T de A, el problema
2
1
62
6
63
6
44
5
6
7
8
9
10
3
11
2 3
127
x1
7
4x2 5
137
7
145
x3
15
x
™
A
˜
2 3
5
657
6 7
7
D 6
657 :
455
5
–
b
La descomposición es
A D
D
U †V T
"
0;3546
0;3987
0;4428
0;4870
0;5311
h
0;2017
0;5168
0;8320
0;6887
0;3756
0;0624
0;2507
0;5638
0;8903
0;2573
0;3757
0;5700
0;7455
0;1702
0;2966
0;0490
0;4082
0;8165
0;4082
0;1764
0;2235
0;3652
0;7652
0;4472
i
0;2096
0;3071
0;7985
0;1627
0;4445
# "35;1272
0
0
0
0
0
2;4654
0
0
0
0
0
0;0000
0
0
#
:
86/97
La solución de norma euclídea mínima se obtiene de
x D
D
uT1 b
1
v1 C
11;0709
35;1272
Ä
uT2 b
v2
2
0;2017
0;5168
0;8320
1;5606
C 2;4654
Ä
0;8903
0;2573
0;3757
D
Ä
0;5
0;0
0;5
:
✉
Este programa de Matlab resuelve el ejemplo llamando a svd( ) para obtener
la descomposición en valores singulares de A.
function [x S r] = Svdre(a,b)
% Resolución ||Ax-b|| mediante la desc. en valores singulares de A
[m,n] = size(a); tol=sqrt(eps); tmp=zeros(m); x=zeros(n,1);
[U S V]=svd(a); S=diag(S);
r=0;
for j=1:n
if S(j)>=tol
r=r+1;
tmp(r)=dot(U(:,j),b)/S(j);
end
end
for j=1:r
x=x+tmp(j)*(V(:,j));
end
87/97
Comparación de los diversos métodos
Método
Operaciones
Ecuaciones Normales
mn2
2
Transformaciones de Householder
mn2
Transformaciones de Givens
2mn2
Método de Gram Schmidt
mn2
Método de Gram Schmidt Modificado
mn2
Método de Golub-Reinsch (SVD)
Método de Golub-Reinsch-Chan (SVD)
3
C n6
n3
3
2 3
n
3
2mn2 C 4n3
mn2 C 17
n3
3
Los métodos basados en transformaciones ortogonales son los más precisos y
habituales.
El basado en SVD es el más robusto, aunque más caro.
88/97
Índice
Introducción
Fundamentos teóricos
Sistemas incompatibles. Ecuaciones normales
Sistemas indeterminados
Resolución numérica del problema
Método de Gram-Schmidt
Método de Gram-Schmidt modificado
Factorización QR
Descomposición numérica en valores singulares
Comparación de los métodos
Matlab y la solución de problemas de mínimos cuadrados
89/97
Matlab y el problema de mínimos
cuadrados
Para calibrar las posibilidades de Matlab con los problemas de mínimos
cuadrados, vamos a utilizarlo para ajustar a unos puntos la función
y D c1xe c2x
que se utiliza en prospección de hidrocarburos y minerales.
Como no disponemos de datos reales, vamos a generar unos sintéticos a partir
de hacer, por ejemplo, c1 D 5 y c2 D 3.
Generaremos 300 puntos y los perturbaremos con un ruido aleatorio
normalizado de media 0.
90/97
Luego utilizaremos como métodos para obtener los parámetros c1 y c2, en un
modelo linealizado:
El operador n
Las ecuaciones normales
La descomposición QR
El método de Gram-Schmidt
La descomposición en valores singulares y
La matriz pseudoinversa,
todos con el auxilio de Matlab.
Para linealizar el modelo, haremos los cambios de variable
v D ln.y=x/; u D x; ˇ D ln c1 y ˛ D c2;
resultando
v D ˛u C ˇ.
91/97
En lo que sigue se lista el “diary” de Matlab que se ha seguido.
>> x0=0.01;
>> x=linspace(x0,2,300);
>> y=5*x.*exp(-3*x);
>> yn=abs(y+0.05*(rand(size(x))-0.5));
>> v=log(yn./x);
>> x=x(:);
>> v=v(:);
>> A=[ones(size(x)) x];
>> c=A\v
c =
1.6443
-3.0579
>> G=chol(A’*A);
>> c1=G\(G’\(A’*v));
>> c1
c1 =
1.6443
-3.0579
>> [Q,R]=qr(A,0);
>> c2=R\(Q’*v)
c2 =
1.6443
-3.0579
%primer punto de muestra
% 300 puntos
%nube de puntos
% + ruido: valores pos.
%cambio de variable
% Matriz del sistema
% Prim. respuesta: con \
% Ecuaciones normales
% Segunda respuesta
% Factorización QR
% Tercera respuesta
92/97
>> [Q,R]=gs_m(A);
>> c3=R\(Q’*v)
c3 =
1.6443
-3.0579
>> format long
>> c4=pinv(A)*v
c4 =
1.64428682050583
-3.05786465731645
>> [U S V]=svd(A);
>> c5=V\(S\(U\v))
c5 =
1.64428682050583
-3.05786465731645
% Gram-Schmidt modi.
% Cuarta respuesta
% Matriz pseudoinversa
% Quinta respuesta
% Descomposición val. sing.
% Sexta respuesta
93/97
Fichero script.m con todas estas instrucciones.
function demoXexp(n)
% demoXexp
Ajuste datos sintéticos a y = c(1)*x*exp(c(2)*x)
%
% Synopsis: demoXexp(n)
%
% Dato: n = (opcional) número de puntos sintéticos a generar.
%
defecto=200
%
if nargin<1, n=200; end
x0=0.01;
x=linspace(x0,2,n);
y=5*x.*exp(-3*x);
yn=abs(y+0.05*(rand(size(x))-0.5));
v=log(yn./x);
x=x(:);
v=v(:);
A=[ones(size(x)) x];
c=A\v;
% Construcción de los datos
%
%
con cambio de variable
%
%
para linealizar el
%
%
modelo
%
fprintf(’Parámetros ajustados:\ncon A\\b:
exp(c(1)),c(2));
c1 = %18.15f c2 = %18.15f\n’,...
94/97
% --- Plot datos
xfit = linspace(min(x),max(x));
yfit = exp(c(1))*xfit.*exp(c(2)*xfit);
if n<30, s = ’v’; else s = ’-’; end
% Símbolo para datos originales
plot(x,y,s,x,yn,’o’,xfit,yfit,’--’);
xlabel(’x’);
ylabel(’y’); legend(’original’,’+ruido’,’ajustado’);
xmax = max(x);
ymax = max(y);
text(0.5*xmax,0.7*ymax,sprintf(’c1 = %6.4f c2 = %6.4f’,exp(c(1)),c(2)));
text(0.5*xmax,0.6*ymax,sprintf(’%d puntos "sintéticos"’,n));
% --- Plot funciones ajustadas
% Cholesky
G=chol(A’*A);
c1=G\(G’\(A’*v));
fprintf(’con chol(A\’’*A);
c1 = %18.15f c2 = %18.15f\n’,exp(c1(1)),c1(2));
% Factorización QR
[Q,R]=qr(A,0);
c2=R\(Q’*v);
fprintf(’con [Q,R]=qr(A);
c1 = %18.15f c2 = %18.15f\n’,exp(c2(1)),c2(2));
95/97
% Gram.Schmidt modificado
[Q,R]=gs_m(A);
c3=R\(Q’*v);
fprintf(’con [Q,R]=gr_m(A);
c1 = %18.15f c2 = %18.15f\n’,exp(c3(1)),c3(2));
% Matriz pseudoinversa
c4=pinv(A)*v;
fprintf(’con pinv(A)*b;
c1 = %18.15f c2 = %18.15f\n’,exp(c4(1)),c4(2));
% Descomposición en valores singulares
[U S V]=svd(A);
c5=V\(S\(U\v));
fprintf(’con [U S V]=svd(A); c1 = %18.15f c2 = %18.15f\n’,exp(c5(1)),c5(2));
96/97
Gráfico con los valores reales y las curvas ajustadas, para 100 puntos sintéticos.
0.7
original
+ruido
ajustado
0.6
0.5
c1 = 5.6363 c2 = −3.1949
0.4
y
100 puntos "sintéticos"
0.3
0.2
0.1
0
0
0.2
0.4
0.6
0.8
1
x
1.2
1.4
1.6
1.8
2
97/97