Homework Problems for Course Numerical Methods for CSE

1 Homework Problems for CourseNumerical Methods for CSER. Hiptmair, G. Alberti, F. LeonardiVersion of July 16, 2016A steady and persistent effort spent on Homework Problems is essential for success inthe should expect to spend 4-6 hours per week on trying to solve the Homework many involve small coding projects, the time it will take an individual student toarrive at a solution is hard to predict. The assignment sheets will be uploaded on the Course webpage on Thursday everyweek. Some or all of the Problems of an assignment sheet will be discussed in the tutorialclasses on Monday 112weeks after the problem sheet has been published. A few Problems on each sheet will be marked ascore Problems .

2 Every participantof the Course is strongly advised to try and solveat leastthe core Problems . If you want your tutor to examine your solution of the current problem sheet, pleaseput it into the plexiglass trays in front of HG G 53/54 by the Thursday after thepublication. You should submit your codes using the online submission is voluntary, but feedback on your performance on Homework Problems can beimportant. You are encouraged to hand-in incomplete and wrong solutions, since you can re-ceive valuable feedback even on incomplete attempts. Please clearly mark the Homework Problems that you want your tutor to R. HiptmairG. Alberti,F. LeonardiAS 2015 Numerical Methods for CSEETH Z richD-MATHP roblem Sheet 0 These Problems are meant as an introduction to EIGENin the first tutorial classes of thenew 1 Gram-Schmidt orthogonalization with EIGEN[1, Code ] presents a MATLAB code that effects the Gram-Schmidt orthogonalizationof the columns of an argument matrix.

3 (1a)Based on the C++ linear algebra library EIGEN implement a functiont e m p l a t e<c l a s sMatrix>Matrix gramschmidt(c o n s tMatrix that performs the same computations as [1, Code ]. (1b)Test your implementation by applying it to a small random matrix and checkingthe orthonormality of the columns of the output 2 Fast matrix multiplication[1, Rem. ] presents Strassen s algorithm that can achieve the multiplication of twodense square matrices of sizen=2k,k N, with an asymptotic complexity better thanO(n3).(2a)Using EIGEN implement a function1 MatrixXd strassenMatMult(c o n s tMatrixXd & A,c o n s tMatrixXd & B)that uses Strassen s algorithm to multiply the two matricesAandBand return the resultas :See Listing 1.)

4 (2b)Validate the correctness of your code by comparing the result with EIGEN sbuilt-in matrix :See Listing 1.(2c)Measure the runtime of your functionstrassenMatMultfor random matri-ces of sizes2k,k=4,..,10, and compare with the matrix multiplication offered by the -operator of :See Listing 1: EIGENI mplementation of the Strassen s algorithm and runtime # i n c l u d e<Eigen / Dense>2# i n c l u d e< i o s t r e a m >3# i n c l u d e< v e c t o r >45# i n c l u d e" t i m e r . h "67u s i n g namespaceEigen ;8u s i n g namespaces t d ;910//! \brief Compute the Matrix productA BusingStrassen's ! \param[in] A Matrix2k 2k12//! \param[in] B Matrix2k 2k13//! \param[out] Matrix product of A and B of dim2k 2k14M a t r i x X d s t r a s s e n M a t M u l t (c o n s tM a t r i x X d & A ,c o n s tM a t r i x X d & B )15{16i n tn=A.}

5 Rows ( ) ;17M a t r i x X d C( n , n ) ;18219i f( n ==2)20{21C<< A ( 0 , 0 )*B ( 0 , 0 ) + A ( 0 , 1 )*B ( 1 , 0 ) ,22A ( 0 , 0 )*B ( 0 , 1 ) + A ( 0 , 1 )*B ( 1 , 1 ) ,23A ( 1 , 0 )*B ( 0 , 0 ) + A ( 1 , 1 )*B ( 1 , 0 ) ,24A ( 1 , 0 )*B ( 0 , 1 ) + A ( 1 , 1 )*B ( 1 , 1 ) ;25r e t u r nC ;26}2728e l s e29{M a t r i x X dQ0 ( n / 2 , n / 2 ) ,Q1 ( n / 2 , n / 2 ) ,Q2 ( n / 2 , n / 2 ) ,Q3 ( n / 2 , n / 2 ) ,30Q4 ( n / 2 , n / 2 ) ,Q5 ( n / 2 , n / 2 ) ,Q6 ( n / 2 , n / 2 ) ;3132M a t r i x X d A11=A . t o p L e f t C o r n e r ( n / 2 , n / 2 ) ;33M a t r i x X d A12=A . t o p R i g h t C o r n e r ( n / 2 , n / 2 ) ;34M a t r i x X d A21=A . b o t t o m L e f t C o r n e r ( n / 2 , n / 2 ) ;35M a t r i x X d A22=A.}

6 B o t t o m R i g h t C o r n e r ( n / 2 , n / 2 ) ;3637M a t r i x X d B11=B . t o p L e f t C o r n e r ( n / 2 , n / 2 ) ;38M a t r i x X d B12=B . t o p R i g h t C o r n e r ( n / 2 , n / 2 ) ;39M a t r i x X d B21=B . b o t t o m L e f t C o r n e r ( n / 2 , n / 2 ) ;40M a t r i x X d B22=B . b o t t o m R i g h t C o r n e r ( n / 2 , n / 2 ) ;4142Q0= s t r a s s e n M a t M u l t ( A11+A22 , B11+B22 ) ;43Q1= s t r a s s e n M a t M u l t ( A21+A22 , B11 ) ;44Q2= s t r a s s e n M a t M u l t ( A11 , B12 B22 ) ;45Q3= s t r a s s e n M a t M u l t ( A22 , B21 B11 ) ;46Q4= s t r a s s e n M a t M u l t ( A11+A12 , B22 ) ;47Q5= s t r a s s e n M a t M u l t ( A21 A11 , B11+B12 ) ;48Q6= s t r a s s e n M a t M u l t ( A12 A22 , B21+B22 ) ;4950C<< Q0+Q3 Q4+Q6 ,51Q2+Q4 ,52Q1+Q3 ,53Q0+Q2 Q1+Q5 ;54r e t u r nC ;55}356}5758i n tmain (v o i d)59{60sra nd ( (u n s i g n e d i n t) t i m e ( 0 ) ) ;6162//check if strassenMatMult works63i n tk = 2.

7 64i n tn=pow ( 2 , k ) ;65M a t r i x X d A= M a t r i x X d : : Random ( n , n ) ;66M a t r i x X d B= M a t r i x X d : : Random ( n , n ) ;67M a t r i x X d AB ( n , n ) , AxB ( n , n ) ;68AB= s t r a s s e n M a t M u l t ( A , B ) ;69 AxB=A*B ;70c o u t <<" U s i n g S t r a s s e n ' s method , A*B= "<<AB<< e n d l ;71c o u t <<" U s i n g s t a n d a r d method , A*B= "<<AxB<< e n d l ;72c o u t <<" The norm o f t h e e r r o r i s"< <(AB AxB ) . norm ( ) << e n d l ;7374//compare runtimes of strassenMatMult and of directmultiplication7576u n s i g n e d i n tr e p e a t s = 1 0 ;77t i m e r <> tm_x , t m _ s t r a s s e n ;78s t d : : v e c t o r <i n t> t i m es _ x , t i m e s _ s t r a s s e n ;7980f o r(u n s i g n e d i n tk = 4 ; k 1 0 ; k ++) {81tm_x.}

8 R e s e t ( ) ;82t m _ s t r a s s e n . r e s e t ( ) ;83f o r(u n s i g n e d i n tr = 0 ; r < r e p e a t s ; ++ r ) {84u n s i g n e d i n tn = pow ( 2 , k ) ;85A = M a t r i x X d : : Random ( n , n ) ;86B = M a t r i x X d : : Random ( n , n ) ;87M a t r i x X d AB ( n , n ) ;8889tm_x . s t a r t ( ) ;90AB=A*B ;91tm_x . s t o p ( ) ;49293t m _ s t r a s s e n . s t a r t ( ) ;94AB= s t r a s s e n M a t M u l t ( A , B ) ;95t m _ s t r a s s e n . s t o p ( ) ;96}97s t d : : c o u t <<" The s t a n d a r d m a t r i x m u l t i p l i c a t i o nt o o k :"<< tm_x . min ( ) . c o u n t ( ) /1000000. <<" ms "<< s t d : : e n d l ;98s t d : : c o u t <<" The S t r a s s e n ' s a l g o r i t h m t o o k :"<< t m _ s t r a s s e n.

9 Min ( ) . c o u n t ( ) /1000000. <<" ms "<< s t d : : e n d l ;99100t i m e s _ x . push_back ( tm_x . min ( ) . c o u n t ( ) ) ;101t i m e s _ s t r a s s e n . push_back (t m _ s t r a s s e n . min ( ) . c o u n t ( ) ) ;102}103104f o r(a u t oi t = t i m e s _ x . b e g i n ( ) ; i t ! = t i m e s _ x . end ( ) ;++ i t ) {105s t d : : c o u t <<*i t <<" ";106}107s t d : : c o u t << s t d : : e n d l ;108f o r(a u t oi t = t i m e s _ s t r a s s e n . b e g i n ( ) ; i t ! =t i m e s _ s t r a s s e n . end ( ) ; ++ i t ) {109s t d : : c o u t <<*i t <<" ";110}111s t d : : c o u t << s t d : : e n d l ;112113}Problem 3 Householder reflectionsThis problem is a supplement to [1, Section ] and related to Gram-Schmidt orthog-onalization, see [1, Code ].

10 Before you tackle this problem, please make sure thatyou remember and understand the notion of a QR-decomposition of a matrix, see [1,Thm. ]. This problem will put to the test your advanced linear algebra skills.(3a)Listing 2 implements a particular 2: MATLAB implementation for Problem 3 in u n c t i o nZ = houserefl(v)2% Porting of to Matlab code3% v is a column vector4% Size of v5n =s i z e(v,1);67w = v/norm(v);8u = w + [1;z e r o s(n-1,1)];9q = u/norm(u);10X =e y e(n) - 2*q*q';1112% Remove first column X(:,1) \in span(v)13Z = X(:,2:end);14endWrite a C++ function with declaration:void houserefl(const VectorXd &v, MatrixXd that is equivalent to the MATLAB functionhouserefl().)