CSCI 2033 sections 001 and 010: Elementary Computational Linear Algebra (2024 Spring) Assignment 1
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit
CSCI 2033 sections 001 and 010: Elementary Computational Linear Algebra (2024 Spring)
Assignment 1
Due 11:59pm, February 27, 2024 on Gradescope. A 25% penalty will be applied for submissions that are up to 24 hours late, and a 50% penalty for submissions up to 48 hours late. Any later submission will not be graded and will get zero automatically.
Notation 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.
Instruction
. This homework set totals 20 points (final grade percentage will be either 8% or 17% depending on the other homework scores);
. We assume that you know basic concepts in a high-level programming language such as Python, C++, Java, Matlab, Julia—this is a prerequisite for this course. But we are using Python throughout this course, because it is the No. 1 language used in modern scientific computing and industrial applications, especially in areas related to modern artificial intelligence, data science, machine/deep learning, computer vision, AR/VR where most modern applications and job positions gravitate. Please find resources to pick up Python yourself; there are tons of options online, for example https://www.pythonlikeyoumeanit.com/index.html.
. Problems 0–2 are designed to familiarize you with NumPy—the de-facto standard for scien- tific computing in Python. Problems 3–4 are about applications using NumPy functions.
. We assume that you are using the Google Colab environment (https://colab.research. google.com/), which provides a convenient and free Jupyter notebook environment ready for computing. Please watch this video tutorial https://youtu.be/oCngVVBSsmAor search and watch tons of similar video tutorials to get started. If you are advanced and comfortable with local installation and running of Jupyter Notebook or JupyterLab (https://jupyter.org/), feel free to do so. But we will not provide support for these and you will need to resolve your own installation and running issues.
. Please show all your work in the 4 Colab files (.ipynb) we release with this homework. Do not modify any provided code and only write your code in regions marked "YOUR CODE STARTS HERE". In your final submission, submit the 4 files separately for their corresponding problems in Gradescope.
You will need to work through the Prob0_Numpy_Tutorial file to master the minimal background necessary to proceed. We will point you to additional tutorial materials as we move on; they are mostly linked from the clickable words and phrases that are in blue.
The problems in this homework are closely related to the textbook of this course — Linear Algebra: Step by Stepby Kuldeep Singh, 2013. In the following problems, we will simply call it the textbook.
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 ⊺ .
2.5 (1.5/5) Matrix power (textbook section 1.5) For any square matrix M ∈ Rn×n, its p-th power is defined naturally as
Mp = MMM...M .
p times (2)
We have two identities for matrix power parallel to those for scalar power:
(Mp )(Mq ) = Mp+q , (Mp )q = Mpq . (3)
Follow the following steps to numerically verify the two identities:
(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;
(iii) call the NumPy built-in functionnumpy.linalg.solveto solve the given linear system to give a solution x2. Ideally, x1 and x2 should be the same. Please calculate the relative distance between x1 and x2, i.e., ∥x1 − x2 ∥ / ∥x2 ∥ . The relative distance should be very close to 0 if your solver works well.
Congratulations! Now you have a simple solver for large linear systems!
2024-03-01