Numerical computing
Question 1:
Copyright By PowCoder代写 加微信 powcoder
Question 2:
Explain what is wrong with A=QhatR and
As noted in the problem, the correct answer is to use 𝑄! instead of 𝑄h𝑎𝑡!. What is wrong is thatQhatisnotsquare,infactitismxn andwhereQismxm.Whenwemultiplybya orthogonal matrix, the norm is preserved. An orthogonal matrix has orthonormal columns and is a square. Thus ||b -Ax||=||𝑄!(𝑏 − 𝐴𝑥)||. As I stated before Qhat has orthonormal columns, but it is not square and thus it will not preserve the norm when multiplied onto(b-Ax).
||b -Ax|| ≠ ||𝑄!(𝑏 − 𝐴𝑥)||. As stated in the question, if used A=QhatR then we would have 0 norm but ||b-Ax|| has norm 𝑄,!b
Question 3:
function [A,q]=house3(A)
%we are modifying the house function in the textbook
% so we don’t have to have to waste computation computing % A column we will already know the result for [m,n]=size(A);q=zeros(1,n);
%this is from the house function
z=A(k:m,k);
el=[1;zeros(m-k,1)]; u=z+sign(z(1))*norm(z)*el; u=u/norm(u);
%above is computing the reflectors
%we need to add 1 to k so we will have(k+1)
%this will index us 1 further into the matrix and doing less %computation A(k:m,k+1:n)=A(k:m,k+1:n)-2*u*(u’*A(k:m,k+1:n));
%we need to make sure the sign is correct %u=z-(+/-||z||*e1)
%we choose the sign based sign(z(1))
%so sign is -sign(z(1))
A(k,k)=-sign(z(1))*norm(z);
q(k)=u(1);
A(k+1:m,k)=u(2:m-k+1);
testing on the input that prof gave us:
Comparing the results from house(A) provided in the textbook and house3(A) the code for this problem
The norm is about machine precision so that means that output generated from house3 is almost equal to the output generated from house
Question 4:
We are asked to change
Yes, it is the same mathematically. If we consider u=B,u’=C and A=D
Then we have B*(C*D) and by the associative property of matrix multiplication this is equal to
(B*C)*D: thus =
Now I will compare the results computationally.
Here is the modified code:
function [A,q]=house4(A)
%same code as house from the textbook [m,n]=size(A);q=zeros(1,n);
z=A(k:m,k);
el=[1;zeros(m-k,1)];
u=z+sign(z(1))*norm(z)*el;
u=u/norm(u);
%We are changing one line here to make it O(m^2n^2) flops A(k:m,k:n)=A(k:m,k:n)-2*(u*u’)*A(k:m,k:n);
q(k)=u(1);
A(k+1:m,k)=u(2:m-k+1);
and the results are:
So both ways of multiplying u u’ and A are equal mathematically and also produce the same result computationally.
the original way was u*(u’*A) A is mxn u’ is 1xm and u is mx1
𝑢1 𝐴11 ⋯ 𝐴1𝑛 𝑢2*{𝑢1 𝑢2 … 𝑢𝑚*2 ⋮ ⋱ ⋮ 7}
… 𝐴𝑚1 ⋯ 𝐴𝑚2 𝑢𝑚
after doing the arithmetic in the brackets we have a 1 x n vector and then we multiply a m x 1 by 1 x n and we have a m x n matrix that took o(mn^2) flops
After changing the order of arithmetic
𝑢1 𝐴11 ⋯ 𝐴1𝑛 {𝑢2*𝑢1 𝑢2 … 𝑢𝑚}*2 ⋮ ⋱ ⋮ 7
… 𝐴𝑚1 ⋯ 𝐴𝑚2 𝑢𝑚
u*u’is mx1*1xm=mxm
then we must multiply a m x m matrix onto a m x n matrix. Each time we multiply the first row of u*u’ onto A we do O(m) flops now we have to multiply n more columns of A which is O(m) flops each and we have multiply m more rows of u*u’ onto A each of which takes O(m) flops
So, (u*u’)*AèO(𝑚”𝑛#) and we must remember that this code is in a loop that runs from 1 to n
So the final flops for computing u*u’)*A in house is O(𝑚”𝑛”)
Graphing the timings:
Here is the code I used to graph, it is a script so there is not function signature:
%this is a matlab script for plottting og=[];%save the values for the original mod=[];%save the values for the modified x=[];%save the values for x axis
for i=1:100
g = @() house(mat);
h = @() house4(mat);
mat=rand(i*10,i);
og=[og timeit(g)]; %get the time and add it to list mod=[mod timeit(h)];%get the time and add it to list x=[x i];
semilogy(x,og,’r-*’)
title(“timings between House and Modified House of matrix sizes m=n*10,n”) xlabel (‘size of n ‘)
ylabel(‘time’)
semilogy(x,mod,’b-*’)
legend(‘Orginal’,’Modified’)
Here is the graph:
In part B I stated that the modified flops is O(𝑚”𝑛”) and the original flops were O(𝑚#𝑛”)
When we created the matrix, we set m=10*n
𝑚”𝑛”/ m𝑛”=m so I am expecting the modified version to be 10 time more flops because m= 10.
If we look at the graph, we the two points differ by a factor of 10 which supports the claim I made in B
Question 5:
Part A) is to write your qr code called myqr
I am acknowledging that professor Overton told me in office hours about doing the multiplication efficiently.
function [Q,R]=myqr(A,info)
%take 1 or 2 inputs but the if there is a second input it has to be 0 economy=1; %this will be the size of m or n depending on if info=0 or not %calling the house function to get all of the reflectors [R,us]=house(A);
[m,n]=size(R);
Q=eye(m,m);
%code to have QR or reduced QR
if nargin == 1
economy=m;
elseif nargin==2
if info==0
Q=eye(m,n);
economy=n;
error(“2nd arg must be 0”) end
%creating R to be returned, setting the values below the diagonal to 0 %we will produce an economy sized version or a full sized version of R for i=1 :n
R(i+1:economy,i)=0;
if(economy==n)
R=R(1:n,:);
%computing Q
%first we have to compute u which is adding the corresponding element from %us (the vector returned by house) to the remaining part of the vector %stored below the diagonal of R
%in order to not compute the full Q matrix each time we will be operating %on a subset of the ID matrix Q each time so we don’t have to compute Q %until we reach the last householder reflector
for k=n :-1:1
u=[us(k);vec(k+1:m,k)]; Q(k:m,k:economy)=Q(k:m,k:economy)-(2*(u*(u’*Q(k:m,k:economy))));
Testing the code from part a)
For full size:
For economy size:
R is upper triangular
Both economy size and full size myqr produce machine precision results which means it is correct.
I will test it on a very large case as well:
I will not go through the other tests because I already proved it above.
For full QR flops:
When we are doing the full QR we have this line:
(2*(u*(u’*Q(k:m,k:economy))))
And in full QR economy = m so,
(2*(u*(u’*Q(k:m,k:m)))). u and u’ both have m elements when k=1 Q is an m x m matrix
So, u’*Q(1:m,1:m) we have each row times col is O(m) flops and we have to do this m times so u’*Q(1:m,1:m) is O(𝑚!) flops and this is in a for loop that goes n times so full QR has O(𝑚!𝑛) flops.
Now for the reduced size QR flops:
When we are doing the reduced QR we have this line:
(2*(u*(u’*Q(k:m,k:economy))))
And in reduced QR economy = n so,
(2*(u*(u’*Q(k:m,k:n)))). u and u’ both have m elements when k=1 Q is an m x n matrix
So, u’*Q(1:m,1:n) we have each row times col is O(m) flops and we have to do this n times so u’*Q(1:m,1:n) is O(𝑚𝑛) flops and this is in a for loop that goes n times so reduced QR has O(𝑚𝑛!) flops.
Code for plotting full QR and when I want reduced QR I change the code slightly :
%this is a script that I created for plotting qr1=[];%storing qr times
mqr=[];%storing myqr times
og=[];%storing times that is known to be O(mn^2) flops x=[];
for i=1:10:150
g = @() qr(mat); % this for for full qr
h = @() myqr(mat); % this for for full qr mat=randn(i*10,i);
qr1=[qr1 timeit(g)]; %getting times qr mqr=[mqr timeit(h)]; %getting times myqr f = @() house(mat);
og=[og timeit(f)]; %getting times house which is O(mn^2) flops
semilogy(x,qr1,’r-*’)
title(“timings between qr and myqr House of matrix sizes m=n*10 x n”) xlabel (‘size of n ‘)
ylabel(‘time’)
semilogy(x,mqr,’b-*’)
semilogy(x,og,’g-*’)
legend(QR’, myqr’,’O(m*n^2) flops’) hold off;
Here is a plot where I had qr myqr and a plot that is known to have O(m*n^2) flops
Here is a plot where I had qr myqr
Looking at both graphs we see that myqr is certainly slower than matlab’s qr. Looking at the first graph, we see the full QR is slower {more flops} than the known time that is O(mn^2)
Taking the ratio of myQR algo with the O(mn^2), is about 26.7518 which supports the claim about the flops of full QR
Now looking at the ratio of the times
Above are the ratios of time and they stay about the same as n increases. So the ratio is about constant.
Now looking at the reduced QR graphs:
Clearly looking at both graphs matlab’s QR is faster. But, looking at the first graph the reduced QR is about the same as the graph that is known to be O(mn^2) flops so reduced myqr is O(mn^2) which I said in part C.
Looking at the ratios now of myqr and qr for reduced size:
Looking at the values, they are about constant as the size of n increases.
And the avg ration of myqr to the O(mn^2) plot is about 2.76 so that definitely supports my claim.
function [A,q]=house6(A)
%this is the same code as house in the book [m,n]=size(A);q=zeros(1,n);
z=A(k:m,k);
el=[1;zeros(m-k,1)];
u=z-sign(z(1))*norm(z)*el; %this is the only line that is different % instead of a + i made it a –
u=u/norm(u);
A(k:m,k:n)=A(k:m,k:n)-2*u*(u’*A(k:m,k:n));
q(k)=u(1);
A(k+1:m,k)=u(2:m-k+1);
I will be using the same code as myqr but I am going to call house6(A) instead of house (A). THAT IS THE ONLY DIFFERENCE
I need to craft a matrix that will give me cancellation.
When crafting the matrix, lets look at u=z-sign(z(1))*norm(z)*el;
That line which is in house. I want to create matrix that has really large diagonal values and small off diagonal values
Here some code to achieve this matrix
A=rand(20,11)/10;B=eye(20,11)*1000;A=A+B
After running the code: lets look at the residual
This is using code that will give us cancelation. And the residual for both full and economy size produced poor results. Interestingly, Q still has orthonormal columns.
Now lets compare this myqr
Q still has orthonormal columns
And using myqr, which uses the original house has much better results.
Question 7:
function [A,q] = houseQ7(A,p)
%we are given information that there p-1 entries below the diagnol are %non-zero and everything else is 0
[m,n]=size(A);q=zeros(1,n);
%this is a modification of the house code from the book
%k+p-1 is non-zero so I set index to be k+p-1
% use index to retrive last row of non-zero for each col
%: what this means A(k:index,k:n)
%we also have to change when we create el vector because we don’t want %operate on elements that will already be 0
%we calculate el by taking the first entry of z and putting p-1 zeros after for k=1:n
index=(k+p-1);
z=A(k:index,k); %using index to retrive last row of non-zero el=[1;zeros(p-1,1)]; % making sure to operate on values we know are 0 u=z+sign(z(1))*norm(z)*el;
u=u/norm(u);
A(k:index,k:n)=A(k:index,k:n)-2*u*(u’*A(k:index ,k:n));%using index again q(k)=u(1);
A(k+1:index,k)=u(2:p); %there are p non-zero elements
This is the matrix A that prof wants us to use
Confirming answer is correct. I compared my results with the originally house function and the norm is machine precision so it is correct.
In order to answer this question, lets look at this line of code:
el=[1;zeros(p-1,1)] A(k:index,k:n)=A(k:index,k:n)-2*u*(u’*A(k:index ,k:n))
U and u’ will always have p elements now and when k=1, u’*A(1:p ,1:n))
then we see that we have to do n multiplications that take O(p) flops each so the multiplication take O(pn) flops and then we have to remember that this code is in a loop that goes n times. So the totally number of flops are O(pn^2)
%this is a matlab script for graphing
og=[];% storing the time for orignal house mod=[];% storing the time for mod house x=[];
for i=1:200
g = @() house(mat);
h = @() houseQ7(mat,p); mat=triu(randn(10*i,i),-(p-1));
og=[og timeit(g)]; %computing time and storing it mod=[mod timeit(h)]; %computing time and storing it x=[x i];
semilogy(x,og,’r-*’)
title(“timings between House and House_P of matrix sizes m=n*10,n”) xlabel (‘size of n ‘)
ylabel(‘time’)
semilogy(x,mod,’b-*’)
legend(‘Orginal House’,’House ModP’)
We see that when use the modified house here it is faster than the originally house which supports my claim above.
Here is another graph with a larger n value.
We can see that as n increases the gap between the modified house and the original house increases. This is because flops for modified house is pn^2 and the flops for originally house depends in mn^2. P is constant, m is changing, thus as n gets large we will have a larger and larger gap.
程序代写 CS代考 加微信: powcoder QQ: 1823890830 Email: powcoder@163.com