Tutorial en Octave para hacer A=QR por medio de Algoritmo de Householder
Учебник, чтобы сделать A = QR с помощью алгоритма домохозяина
Tutorial to do A = QR by means of Householder Algorithm
Autor: Alexander Arias
Fecha: 19 de abril de 2019
El método de Householder consiste en construir matrices unitariaselementales (reflexiones) que operan sucesivamente sobre las columnas de A para reducir esta matriz a una forma triangular superior.
https://www.inf.utfsm.cl/~parce/cc1/clase14-RP.html
http://www.ehu.eus/izaballa/Ana_Matr/Apuntes/lec7.pdf
Unitary Triangularization of a Nonsymmetric Matrix
https://www.stat.uchicago.edu/…/302/classics/householder.pdf
#***************************** *****************************
#A continuación empieza el código para Octave o matlab
#***************************** *****************************
disp("Tutorial para hacer A=Qr por medio de Algoritmo de Householder");
#Autor: Alexander Arias
#Fecha: 19 de abril de 2019
#Ejemplo de A=QR por medio de transformaciones de householder:
#Algunos items teóricos se tomarón de poujh:
#https://www.youtube.com/watch?v=oOHNgFtBFUo
#https://www.youtube.com/watch?v=6VTLJQ2v7So
#Obtener vectores base e:
#e = @(x) eye(size(A))(:,x);
#e(1)
disp('A=QR por medio de transformaciones de householder')
A=[2 1 2; 2 2 1; 1 2 2]
disp('-------------------------------------------------')
disp('Encontrar A2=H1*A')
a1=A(:,1)
e = eye(size(A))(:,1);
u=a1-norm(a1)*e
v=u/norm(u)
vvt=v*v'
E=eye(size(A));
H1=E-2*vvt
A2=H1*A
disp('-------------------------------------------------')
disp('Encontrar R=H2*A2')
A2p=A2(2:end,2:end)
a2p=A2p(:,1)
e = eye(size(A2p))(:,1)
u=a2p-norm(a2p)*e
v=u/norm(u)
vvt=v*v'
E=eye(size(A2p));
H2p=E-2*vvt
H2=eye(size(A));
H2(2:end,2:end)=H2p
disp('La R sin redondeo es:')
R=H2*A2
disp('-------------------------------------------------')
disp('Primera forma de redondeo con función round')
Rp=round(mymatrix .* 10000) ./ 10000
disp('-------------------------------------------------')
disp('Segunda forma de redondeo con archivos')
#Primero escribo la matriz en el archivo deseado con el formato que quiero
filename = "Householder_3x3_poujh.txt";
fid = fopen (filename, "a+");
fprintf(fid,"%4.4f %4.4f %4.4f\n", R);
fclose (fid);
#Luego la leo y la guardo en una segunda matriz
filename = "Householder_3x3_poujh.txt";
fid = fopen (filename, "r");
Rk=fscanf(fid,"%f",[3,3])
fclose (fid);
#Borrar el contenido completo de un archivo
#Lo abro con w y lo cierro
#filename = "Householder_3x3_poujh.txt";
#fid = fopen (filename, "w");
#fclose (fid);
disp('-------------------------------------------------')
disp('Encontramos Q=H1t * H2t');
Q=H1'*H2'
disp('-------------------------------------------------')
disp('Desplegamos A=QR calculado con el algoritmo de Householder paso a paso');
Q
Rk
disp('Rectificamos que el producto sea A=QR');
Q*Rk
disp('-------------------------------------------------')
disp('Ahora con el algoritmo de octave [Q,R]=qr(A)');
[Q,R]=qr(A)
disp('Rectificamos que el producto sea A=QR');
Q*R
#***************************** *****************************
#Aqui termina el código
#***************************** *****************************
#***************************** *****************************
#A continuación se muestra el resultado
#***************************** *****************************
https://www.stat.uchicago.edu/…/302/classics/householder.pdf
#***************************** *****************************
#A continuación empieza el código para Octave o matlab
#***************************** *****************************
disp("Tutorial para hacer A=Qr por medio de Algoritmo de Householder");
#Autor: Alexander Arias
#Fecha: 19 de abril de 2019
#Ejemplo de A=QR por medio de transformaciones de householder:
#Algunos items teóricos se tomarón de poujh:
#https://www.youtube.com/watch?v=oOHNgFtBFUo
#https://www.youtube.com/watch?v=6VTLJQ2v7So
#Obtener vectores base e:
#e = @(x) eye(size(A))(:,x);
#e(1)
disp('A=QR por medio de transformaciones de householder')
A=[2 1 2; 2 2 1; 1 2 2]
disp('-------------------------------------------------')
disp('Encontrar A2=H1*A')
a1=A(:,1)
e = eye(size(A))(:,1);
u=a1-norm(a1)*e
v=u/norm(u)
vvt=v*v'
E=eye(size(A));
H1=E-2*vvt
A2=H1*A
disp('-------------------------------------------------')
disp('Encontrar R=H2*A2')
A2p=A2(2:end,2:end)
a2p=A2p(:,1)
e = eye(size(A2p))(:,1)
u=a2p-norm(a2p)*e
v=u/norm(u)
vvt=v*v'
E=eye(size(A2p));
H2p=E-2*vvt
H2=eye(size(A));
H2(2:end,2:end)=H2p
disp('La R sin redondeo es:')
R=H2*A2
disp('-------------------------------------------------')
disp('Primera forma de redondeo con función round')
Rp=round(mymatrix .* 10000) ./ 10000
disp('-------------------------------------------------')
disp('Segunda forma de redondeo con archivos')
#Primero escribo la matriz en el archivo deseado con el formato que quiero
filename = "Householder_3x3_poujh.txt";
fid = fopen (filename, "a+");
fprintf(fid,"%4.4f %4.4f %4.4f\n", R);
fclose (fid);
#Luego la leo y la guardo en una segunda matriz
filename = "Householder_3x3_poujh.txt";
fid = fopen (filename, "r");
Rk=fscanf(fid,"%f",[3,3])
fclose (fid);
#Borrar el contenido completo de un archivo
#Lo abro con w y lo cierro
#filename = "Householder_3x3_poujh.txt";
#fid = fopen (filename, "w");
#fclose (fid);
disp('-------------------------------------------------')
disp('Encontramos Q=H1t * H2t');
Q=H1'*H2'
disp('-------------------------------------------------')
disp('Desplegamos A=QR calculado con el algoritmo de Householder paso a paso');
Q
Rk
disp('Rectificamos que el producto sea A=QR');
Q*Rk
disp('-------------------------------------------------')
disp('Ahora con el algoritmo de octave [Q,R]=qr(A)');
[Q,R]=qr(A)
disp('Rectificamos que el producto sea A=QR');
Q*R
#***************************** *****************************
#Aqui termina el código
#***************************** *****************************
#***************************** *****************************
#A continuación se muestra el resultado
#***************************** *****************************
octave:94> source("Householder_3x3_poujh.m") Tutorial para hacer A=Qr por medio de Algoritmo de Householder A=QR por medio de transformaciones de householder A = 2 1 2 2 2 1 1 2 2 ------------------------------------------------- Encontrar A2=H1*A a1 = 2 2 1 u = -1 2 1 v = -0.40825 0.81650 0.40825 vvt = 0.16667 -0.33333 -0.16667 -0.33333 0.66667 0.33333 -0.16667 0.33333 0.16667 H1 = 0.66667 0.66667 0.33333 0.66667 -0.33333 -0.66667 0.33333 -0.66667 0.66667 A2 = 3.0000e+00 2.6667e+00 2.6667e+00 -5.5511e-16 -1.3333e+00 -3.3333e-01 -3.3307e-16 3.3333e-01 1.3333e+00 ------------------------------------------------- Encontrar R=H2*A2 A2p = -1.33333 -0.33333 0.33333 1.33333 a2p = -1.33333 0.33333 e = Diagonal Matrix 1 0 u = -2.70770 0.33333 v = -0.99251 0.12218 vvt = 0.985071 -0.121268 -0.121268 0.014929 H2p = -0.97014 0.24254 0.24254 0.97014 H2 = 1.00000 0.00000 0.00000 0.00000 -0.97014 0.24254 0.00000 0.24254 0.97014 La R sin redondeo es: R = 3.0000e+00 2.6667e+00 2.6667e+00 4.5776e-16 1.3744e+00 6.4676e-01 -4.5776e-16 -1.6653e-16 1.2127e+00 ------------------------------------------------- Primera forma de redondeo con función round Rp = 3.00000 2.66670 2.66670 0.00000 1.37440 0.64680 -0.00000 -0.00000 1.21270 ------------------------------------------------- Segunda forma de redondeo con archivos Rk = 3.00000 2.66670 2.66670 0.00000 1.37440 0.64680 -0.00000 -0.00000 1.21270 ------------------------------------------------- Encontramos Q=H1t * H2t Q = 0.66667 -0.56592 0.48507 0.66667 0.16169 -0.72761 0.33333 0.80845 0.48507 ------------------------------------------------- Desplegamos A=QR calculado con el algoritmo de Householder paso a paso Q = 0.66667 -0.56592 0.48507 0.66667 0.16169 -0.72761 0.33333 0.80845 0.48507 Rk = 3.00000 2.66670 2.66670 0.00000 1.37440 0.64680 -0.00000 -0.00000 1.21270 Rectificamos que el producto sea A=QR ans = 2.0000 1.0000 2.0000 2.0000 2.0000 1.0000 1.0000 2.0000 2.0001 ------------------------------------------------- Ahora con el algoritmo de octave [Q,R]=qr(A) Q = -0.66667 0.56592 0.48507 -0.66667 -0.16169 -0.72761 -0.33333 -0.80845 0.48507 R = -3.00000 -2.66667 -2.66667 0.00000 -1.37437 -0.64676 0.00000 0.00000 1.21268 Rectificamos que el producto sea A=QR ans = 2.0000 1.0000 2.0000 2.0000 2.0000 1.0000 1.0000 2.0000 2.0000
No hay comentarios:
Publicar un comentario