CoCalc -- 2018-03-09.ipynb

| Download

This repository contains the course materials from Math 157: Intro to Mathematical Software.

Creative Commons BY-SA 4.0 license.

Project: Support and Testing

Path: math157 / assignments / 2018-03-09 / 2018-03-09.ipynb

Views: ³⁰³⁷
License: OTHER

Kernel: Python 3 (Ubuntu Linux)

Math 157: Intro to Mathematical Software

UC San Diego, winter 2018

Homework 8: due March 10, 2018

Note: Due to the cancellation of sections and rescheduling of office hours this week, I am preemptively moving the due date of this assignment to Saturday, March 10 at 8pm.

This problem set consists of 6 problems, each of equal value. As usual, please enter all answers within this notebook except as specified, and cite all sources and collaborators.

Throughout this problem set, use the Julia kernel.

Problem 1: Parallel computing

Grading criteria: correctness and thoroughness of explanations. No code required.

1a. Explain what the "GIL" is in Python and its effect on parallel computation.

In [0]:

1b. Explain what a "race condition" is and how the absence of a GIL pertains to this.

In [0]:

Problem 2: Euclidean algorithm

Grading criteria: correctness of code.

2a. Implement the Euclidean algorithm from scratch (not using any built-in functions other than arithmetic) for multiprecision integers.

In [0]:

2b. Implement the Chinese remainder theorem (for two moduli) using your answer to 2a. You should implement a function mimicking Sage's crt function, except that if the moduli are not coprime you may simply return an error. (By contrast, Sage handles this case more robustly: it returns a meaningful error when possible and returns an error otherwise.)

In [0]:

Problem 3: MATLAB vs. Julia

Grading criterion: correctness of code.

This exercise refers to the Math 18 MATLAB exercise set.

3a. Do MATLAB exercise 4.5 once more, this time using Julia.

3b. Repeat with MATLAB exercise 5.6, using the Moore-Penrose pseudoinverse to perform the least squares computation.

In [0]:

Problem 4: Hilbert matrices

Grading criterion: correctness of code and explanations.

4a. Define a Julia function f, which on the input of a positive integer $n$ , returns the $n \times n$ Hilbert matrix with floating-point entries. Your function should work directly from the definition; however, you may use Nemo's built-in function to check your answer.

In [0]:

4b. Compute the inverse of f(25).

In [0]:

4c. Is Julia's default algorithm for matrix inverse numerically stable? Justify your claim in terms of your answer to 4b.

In [0]:

4d. In light of the Julia mission statement:

Julia is a high-level, high-performance dynamic programming language for numerical computing.

explain why your answer to 4c is to be expected (that is, why the Julia designers chose as they did in this case).

In [0]:

Problem 5: Hadamard matrices

Grading criterion: correctness of code and explanations.

For this exercise, you will use the Nemo library.

In [4]:

using Nemo

  File "<ipython-input-4-bd0227f09cca>", line 1
    using Nemo
             ^
SyntaxError: invalid syntax

5a. Using Nemo's built-in function, construct a $28 \times 28$ Hadamard matrix. Call this matrix $H$ .

In [2]:

M = MatrixSpace(ZZ, 28, 28)
H = hadamard(M)

[1 -1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
[-1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1 1 -1]
[1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1]
[1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1]
[1 1 1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1]
[1 -1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1]
[1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1]
[1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1]
[1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1]
[1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1]
[1 1 1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1]
[1 -1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1]
[1 1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1]
[1 -1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1]
[1 1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1]
[1 -1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1]
[1 1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 1 -1 -1]
[1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1 -1 1]
[1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1 1 1]
[1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1 1 -1]
[1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1 1 1]
[1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1 1 -1]
[1 1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1 1 1 -1 -1]
[1 -1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1 -1 1]
[1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1 1 1]
[1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1 1 -1]
[1 1 1 1 -1 -1 1 1 1 1 -1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 -1 1 1 1 -1]
[1 -1 1 -1 -1 1 1 -1 1 -1 -1 1 -1 1 -1 1 -1 1 1 -1 1 -1 -1 1 1 -1 -1 -1]

In [7]:

A = [[Int(H[i,j]) for i in 1:28] for j in 1:28]

28-element Array{Array{Int64,1},1}:
 [1, -1, 1, 1, 1, 1, 1, 1, 1, 1  …  1, 1, 1, 1, 1, 1, 1, 1, 1, 1]          
 [-1, -1, 1, -1, 1, -1, 1, -1, 1, -1  …  1, -1, 1, -1, 1, -1, 1, -1, 1, -1]
 [1, 1, 1, -1, 1, 1, -1, -1, 1, 1  …  -1, -1, 1, 1, 1, 1, -1, -1, 1, 1]    
 [1, -1, -1, -1, 1, -1, -1, 1, 1, -1  …  -1, 1, 1, -1, 1, -1, -1, 1, 1, -1]
 [1, 1, 1, 1, 1, -1, 1, 1, -1, -1  …  -1, -1, -1, -1, 1, 1, 1, 1, -1, -1]  
 [1, -1, 1, -1, -1, -1, 1, -1, -1, 1  …  -1, 1, -1, 1, 1, -1, 1, -1, -1, 1]
 [1, 1, -1, -1, 1, 1, 1, -1, 1, 1  …  -1, -1, -1, -1, -1, -1, 1, 1, 1, 1]  
 [1, -1, -1, 1, 1, -1, -1, -1, 1, -1  …  -1, 1, -1, 1, -1, 1, 1, -1, 1, -1]
 [1, 1, 1, 1, -1, -1, 1, 1, 1, -1  …  -1, -1, -1, -1, -1, -1, -1, -1, 1, 1]
 [1, -1, 1, -1, -1, 1, 1, -1, -1, -1  …  -1, 1, -1, 1, -1, 1, -1, 1, 1, -1]
 [1, 1, 1, 1, 1, 1, -1, -1, 1, 1  …  1, 1, -1, -1, -1, -1, -1, -1, -1, -1] 
 [1, -1, 1, -1, 1, -1, -1, 1, 1, -1  …  1, -1, -1, 1, -1, 1, -1, 1, -1, 1] 
 [1, 1, -1, -1, 1, 1, 1, 1, -1, -1  …  1, 1, 1, 1, -1, -1, -1, -1, -1, -1] 
 ⋮                                                                         
 [1, 1, -1, -1, -1, -1, -1, -1, 1, 1  …  1, 1, -1, -1, 1, 1, 1, 1, -1, -1] 
 [1, -1, -1, 1, -1, 1, -1, 1, 1, -1  …  1, -1, -1, 1, 1, -1, 1, -1, -1, 1] 
 [1, 1, -1, -1, -1, -1, -1, -1, -1, -1  …  1, -1, 1, 1, -1, -1, 1, 1, 1, 1]
 [1, -1, -1, 1, -1, 1, -1, 1, -1, 1  …  -1, -1, 1, -1, -1, 1, 1, -1, 1, -1]
 [1, 1, 1, 1, -1, -1, -1, -1, -1, -1  …  1, 1, 1, -1, 1, 1, -1, -1, 1, 1]  
 [1, -1, 1, -1, -1, 1, -1, 1, -1, 1  …  1, -1, -1, -1, 1, -1, -1, 1, 1, -1]
 [1, 1, 1, 1, 1, 1, -1, -1, -1, -1  …  -1, -1, 1, 1, 1, -1, 1, 1, -1, -1]  
 [1, -1, 1, -1, 1, -1, -1, 1, -1, 1  …  -1, 1, 1, -1, -1, -1, 1, -1, -1, 1]
 [1, 1, -1, -1, 1, 1, 1, 1, -1, -1  …  1, 1, -1, -1, 1, 1, 1, -1, 1, 1]    
 [1, -1, -1, 1, 1, -1, 1, -1, -1, 1  …  1, -1, -1, 1, 1, -1, -1, -1, 1, -1]
 [1, 1, 1, 1, -1, -1, 1, 1, 1, 1  …  1, 1, 1, 1, -1, -1, 1, 1, 1, -1]      
 [1, -1, 1, -1, -1, 1, 1, -1, 1, -1  …  1, -1, 1, -1, -1, 1, 1, -1, -1, -1]

5b. Compute the determinant of $H$ and check that it achieves equality in the Hadamard bound.

In [0]:

5c. Compute the dot products between all pairs of rows in $H$ .

In [0]:

5d. Explain the relationship between your answers to 5b and 5c.

In [0]:

Problem 6: Singular value decomposition

Grading criterion: correctness and pertinence of code. (Correctness means it runs, pertinence means you followed the instructions.)

Consider the following example from the numpy documentation.

>>> a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
>>> b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
>>> u, s, vh = np.linalg.svd(a, full_matrices=True)
>>> u.shape, s.shape, vh.shape
((9, 9), (6,), (6, 6))
>>> np.allclose(a, np.dot(u[:, :6] * s, vh))
True
>>> smat = np.zeros((9, 6), dtype=complex)
>>> smat[:6, :6] = np.diag(s)
>>> np.allclose(a, np.dot(u, np.dot(smat, vh)))
True

Write Julia code to emulate this example as faithfully as possible.

In [6]:

import numpy as np

In [9]:

a = np.random.randn(9, 6) + 1j*np.random.randn(9, 6)
b = np.random.randn(2, 7, 8, 3) + 1j*np.random.randn(2, 7, 8, 3)
u, s, vh = np.linalg.svd(a, full_matrices=True)
u.shape, s.shape, vh.shape

((9, 9), (6,), (6, 6))

In [0]:

In [14]:

c = u[:,:6]*s

In [15]:

c.shape, vh.shape

((9, 6), (6, 6))

In [0]: