CSCI 2033 sections 001 and 010: Elementary Computational Linear Algebra (2024 Spring)

CSCI 2033 sections 001 and 010:  Elementary Computational Linear Algebra (2024 Spring)

Assignment 1

Problem 1 Vector Operations (5 points)

Create 3 random vectors u, v , w E R10000 as follows:

1    import   numpy   a s  np

2   rng  =  np . random . default _ rng (20232033)   #  fix   a   random   seed .   Please   do  not   modify   it

3   u  =  rng . random ((10000   ,1) )  #   generate   random   vector  u

4   v  =  rng . random ((10000   ,1) )  #   generate   random   vector  v

5   w  =  rng . random ((10000   ,1) )  #   generate   random   vector  w

we will use these vectors for all the following questions in Problem 1.

1.1 (1.5/5) Vector indexing and concatenation (textbook section 1.3)     Please obtain the following element or subvectors; we have provided some examples in the Prob0_Numpy_Tutorial file:

(a) The 2023rd element of vector u. NOTE: Python/NumPy indexing starts from 0 instead of 1;

(b) The 2023rd to 2033rd elements of vector v (including the 2023rd and 2033rd elements). NOTE: Python/NumPy indexing will not include the last element in indexing. Make sure that the size of the subvector you obtain is 11. You may want to use the built-in numpy.ndarray.shape to help you check the size of your subvector;

(c) Make a new vector by combining the first 30 elements of v and the last 100 elements of w. You need to use the Numpy built-in function numpy.concatenate.

Note: If want to learn more about this, you can go to this NumPy tutorial.

1.2 (1/5) Linear combinations (textbook section 1.3)    Calculate the following linear combinations:

u + v + w ,    2u + 3v + 3w.

1.3 (1.5/5) Inner products (textbook section 1.3)    Calculate the following inner products using the built-in function numpy.inner:

〈u, u〉,  〈u - 2v , w〉,  〈3u, 2v + w〉.

1.4 (1/5) Vector norms (textbook section 2.1)    Calculate the following vector norms using the NumPy built-in function numpy.linalg.norm:

IuI ,    Iv + 3wI .

Problem 2.  Matrix Operations (5 points)

Reminder about our notation and convention: Scalars are small letters (e.g., a,b, λ, α, β), vectors are boldface small letters (e.g., v , w), and matrices are boldface capital letters (e.g., A, B). Vectors by default are column vectors; they are matrices with single columns. Row vectors are matrices with single rows.

We start by generating a few random matrices and vectors:

1 import numpy as np

2 rng = np . random . default_rng (20232033) # fix a random seed . Please do not modify it

3 A = rng . random ((100 ,100) ) # generate random matrix A

4 B = rng . random ((100 ,200) ) # generate random matrix B

5 C = rng . random ((100 ,200) ) # generate random matrix C

6 D = rng . random ((100 ,100) ) # generate random matrix D

7 u = rng . random ((100 ,1) ) # generate random vector u

8 v = rng . random ((200 ,1) ) # generate random vector v

We will use these matrices and vectors for all the following questions in Problem 2. We also provided some examples in Prob0_Numpy_Tutorial file.

2.1 (0.5/5) Matrix norms (textbook section 2.1 & section 1.6)     The magnitude of a matrix can be measured similarly to that of vectors. For any matrix M E Rm×n, its (Frobenius) norm is defined as

IM IF  = ‘〈M , M〉,                                                          (1)

where F is for Frobenius (a famous German mathematician). Call the NumPy built-in function numpy.linalg.norm, and calculate the following

(a)  IAIF,

(b)  IB - CIF . This is the distance between B and C.

2.2 (0.5/5) Matrix indexing (Discussion Session)     Please obtain these submatrices:

(a) The top-left 50-by-50 submatrix of A;

(b) The bottom-right 30-by-25 submatrix of B.

Note: If want to learn more about this, you can go to this NumPy tutorial.

2.3 (0.5/5) Matrix-vector multiplication (textbook section 1.4)    Calculate the following matrix- vector multiplication using the built-in function numpy.matmul (for matrix multiplication) and numpy.transpose (for matrix transpose). NOTE: The @ operator can be used as a shorthand for NumPy.matmul on ndarrays; M.T can be used as a shorthand for NumPy.transpose of matrix M:

Au, C⊺u, Bv.

2.4  (0.5/5) Matrix-matrix multiplication (textbook section 1.4 & section 1.6)     Calculate the following matrix-matrix multiplication using the built-in function numpy.matmul (for matrix multiplication) and numpy.transpose (for matrix transpose). NOTE: The @ operator can be used as a shorthand for NumPy.matmul on ndarrays; M.T can be used as a shorthand for NumPy.transpose of matrix M:

AB , BC⊺ , C⊺B , uv ⊺ .

(a) Implement your own matrix power function mat_pow(): it should take any square matrix M and the integer power p ≥ 0, and output the values of the matrix Mp. NOTE: To debug, you are encouraged to test your implementation against the Numpy built-in matrix power function numpy.linalg.matrix_power. But, this is not required in your submission.

(b) Use your own mat_pow() function to calculate (A6)(A8) and A6+8, and also calculate the relative distance (see definition below) between (A6)(A8) and A6+8 — the relative distance should be very close to 0;

(c) Using your own mat_pow() function to calculate (A6)8 and A6∗8, and also calculate the relative distance between (A6)8 and A6∗8 — the relative distance should be very close to 0.

Definition: relative distance of matrices M and N of the same size equals  .

2.6 (1.5/5) Inverse and transpose of matrices (textbook section 1.6)    Complete the following calculations using the NumPy built-in function numpy.linalg.inv (for matrix inverse):

(a)  (AD)−1 and D−1 A−1, and the relative distance between them—the relative distance should be very close to 0;

(b)  (A−1)⊺ and  (A⊺)−1, and the relative distance between them—the relative distance should be very close to 0;

(c)  (AB)⊺ and B⊺A⊺, and the relative distance between them—the relative distance should be very close to 0.

Problem 3.  Gaussian Elimination and Back Substitution (5 points)

In this problem, we will implement Gaussian elimination and back substitution. In the end, we will solve a large linear system Ax = busing our implementation. The Gaussian elimination algorithm is largely based on Section 1.2 of the textbook; we make small necessary changes to ensure that it works reliably on computers. Check the Colab file Prob3_Gaussian_Elimination n Back_Substitution for code template.

3.0 (0/5) Preparation Gaussian elimination involves three types of row operations:

(a) Multiply a row by a non-zero factor. For example, multiplying λ (λ  0) on the i-row to produce the new i-th row can be written as

1                  M [[ i ] ,:]  =   lamb * M [[ i ] ,:]

(b) Subtract a multiple of a top row from a bottom row. For example, subtracting λ times the i-th row from the j-th row of M, where i < j, to produce the new j-th row, can be written as

1                  M [[ j ] ,:]  =  M [[ j ] ,:]  -  lamb * M [[ i ] ,:]

(c) Exchanging rows. For example, exchanging the i-th and j-th row of the matrix M can be written as

1                  M [[ i ,  j ] ,:]  =  M [[ j ,  i ] ,:]

3.1 (1.5/5) Gaussian elimination (Version 0) (textbook section 1.2)     Implement Gaussian elimi-

Algorithm 1 Gaussian Elimination Version 0

nation following the pseudocode in Algorithm 1. Your function should be called gauss_elim_v0 that:  (i) takes an square matrix A e Rn×n, a vector b e Rn, and a print flag print_flag that controls whether we print the intermediate augmented matrix after each row operation, and (ii) returns a matrix U e Rn×(n+1) where the left nx nsubmatrix of U is in the row echelon form. Hint: Suppose that two matrices M and N have the same number of rows. To concatenate them in the horizontal direction, we can call the built-in function numpy.concatenate:

1                  P  =  np . concatenate (( M , N ) , axis =1)

To test your implementation, let us take a test case

A =  [' 2(2) 1   3(3)l' , b =  [' 5(4)l' ,   and hence U =  [' 2(2) 1   3(3)   5(4)l' .                 (4)

Your Gaussian elimination should produce the following sequence of intermediate augmented matrices in the right order (Note: the elements marked red are the leading numbers that we are currently using to eliminate non-zeros below them):

2(1) - 1(1)   3(1)   4(1) R1=R1−2R喻0 0(1) 1   1(1)   2(1) R2=R2−2R喻0 0(1) 1   1(1)   2(1)

[2     0     3   5l                      [2     0     3   5l                      [0     2     1   3l

To get full credit, you need to print out the intermediate augmented matrix after each row operation.

3.2 (2/5) Back substitution (textbook section 1.2)    We first implement back substitution, and then combine Gaussian elimination and back substitution into a linear system solver for cases where A is square. Finally, we test our linear solver against the Numpy built-in.

Algorithm 2 Backward Substitution

(a) Implement back substitution following the pseudocode inAlgorithm 2. Your function should be called back_subs that: (i) takes an augmented matrix U e Rn×(n+1) in the row echelon form, and a print flag print_flag that controls whether we print the newly solved variable value after each substitution step, and (ii) returns an x e Rn as a solution to Ax = b. As a test, take our previous final augmented matrix in Eq. (5), back substitution should give us

R2 : x2  = ( — 1)/( — 1) = 1

R1 : x1  = (2 — 1 * 1)/1 = 1                                                          (6)

R0 : x0  = (1 — ( — 1) * 1 — 1 * 1)/1 = 1

as we move from bottom to top, row by row. To get full credit, you need to print out the intermediate newly solved variable after each substitution step (i.e., x2, x1, and x0 in our test).

(b) Implement a function my_solver_v0 by combining the gauss_elim_v0 and back_subs func- tions implemented above: this function takes a square matrix A e Rn×n and a vector b e Rn, and returns a vector x e Rn so that Ax = b. In other words, my_solver_v0 solves the linear system Ax = b when given A and b. To test your solver, in the code template, we provide a randomly generated A e R300×300 and b e R300. Please

(i) solve the given 300 x 300 linear system using your solver—we will denote this solution by x1;

(ii) validate your solution x1  by calculating the relative error IAx1 — bI / IAIF, which should be very close to 0 if your solver works well;