We provide signatures of the functions that you have to implement. Make sure you follow the signatures defined, otherwise your coding solutions will not be graded.
Read homework rules carefully. If you do not follow it you will likely be penalized.
In [1]:
from scipy.linalg import toeplitz
import numpy as np
import math
import scipy.io.wavfile as wav
import matplotlib.pyplot as plt
from IPython.display import Audio
%matplotlib notebook
In [2]:
# reading
rate, audio = wav.read("TMaRdy00.wav")
# plotting
plt.plot(audio)
plt.ylabel("Amplitude")
plt.xlabel("Time")
plt.title("You wanna piece of me, boy?")
plt.show()
# playing
Audio(audio, rate=rate)
Out[2]:
Our next goal is to process this signal by multiplying it by a special type of matrix (convolution operation) that will smooth the signal.
signal[::p].
In [3]:
N = len(audio) ** 2 * 8 / 1024 / 1024
print(N)
audio = audio[::4]
def gen_toeplitz(N, alpha):
return T
In [4]:
# INPUT: N - integer (positive), alpha - float (positive)
# OUTPUT: T - np.array (shape: NxN)
def gen_toeplitz(N, alpha):
i = np.linspace(1, N, N)
j = np.linspace(1, N, N)
im, jm = np.meshgrid(i, j, indexing="ij")
T = np.sqrt(alpha / np.pi) * np.exp(-alpha * (im - jm)**2)
return T
In [5]:
# INPUT: signal - np.array (shape: Nx1), N - int (positive), alpha - float (positive)
# OUTPUT: convolved_signal - np.array (shape: Nx1)
def convolution(signal, N, alpha):
T = gen_toeplitz(N, alpha)
convolved_signal = np.matmul(T, signal)
return convolved_signal
(3 pts) Plot the first $100$ points of the result and the first $100$ points of your signal on the same figure. Do the same plots for $\alpha = \frac{1}{5}$, $\alpha = \frac{1}{100}$ using plt.subplots in matplotlib. Each subplot should contain first $100$ points of initial and convolved signals for some $\alpha$. Make sure that you got results that look like smoothed initial signal.
(2 pts) Play the resulting signal. In order to do so you should also scale the frequency (rate), which is one of the inputs in Audio.
Note that you cannot play a signal which is too small.
In [6]:
audio_alpha1 = convolution(audio, len(audio), 0.2)
audio_alpha2 = convolution(audio, len(audio), 0.01)
# plotting
fig = plt.figure(figsize=(15,10))
fig.tight_layout()
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)
ax1.plot(audio[:100:], label="Original")
ax1.plot(audio_alpha1[:100:], label="Convolved")
ax1.set_ylabel("Amplitude", size=13)
ax1.set_xlabel("Time", size=13)
ax1.set_title("Original and convolved signals (alpha = 0.2)", size=13)
ax1.legend(loc="upper right")
ax2.plot(audio[:100:], label="Original")
ax2.plot(audio_alpha2[:100:], label="Convolved")
ax2.set_ylabel("Amplitude", size=13)
ax2.set_xlabel("Time", size=13)
ax2.set_title("Original and convolved signals (alpha = 0.01)", size=13)
ax2.legend(loc="upper right")
Out[6]:
In [7]:
# As we put every 4 sample from original audio track we need to descrease sample rate in 4 times
Audio(audio_alpha1, rate=rate/4)
Out[7]:
In [ ]:
Audio(audio_alpha2, rate=rate/4)
Given a convolved signal $y$ and an initial signal $x$ our goal now is to recover $x$ by solving the system $$ y = Tx. $$ To do so we will run iterative process $$ x_{k+1} = x_{k} - \tau_k (Tx_k - y), \quad k=1,2,\dots $$ starting from zero vector $x_0$. There are different ways how to define parameters $\tau_k$. Different choices lead to different methods (e.g. Richardson iteration, Chebyshev iteration, etc.). This topic will be covered in details later in our course.
To get some intuition why this process converges to the solution of $Tx=y$, we can consider the following. Let us note that if $x_k$ converges to some limit $x$, then so does $x_{k+1}$. Taking $k\to \infty$ we arrive at $x = x - \tau (Tx - y)$ and hence $x$ is the solution of $Tx = y$.
Another important point is that iterative process requires only matrix-vector porducts $Tx_k$ on each iteration instead of the whole matrix. In this problem we, however, work with the full matrix, but keep in mind, that convolution can be done efficiently without storing the whole matrix.
First of all, let's work out the expression $Tx_{k+1} - y$: $$ Tx_{k+1} - y = T(x_k - \tau_k(Tx_k-y))-y = T(x_k - \tau_k T x_k + \tau_k y) - y = Tx_k - \tau_k T^2 x_k + \tau_kTy - y = Tx_k - y + \tau_k(Ty - T^2x_k) = r_k - \tau_k Tr_k $$ So $$ \|r_k - \tau_k Tr_k\|_2 \to \min_{\tau_k} $$ Then expression can be rewritten as: $$ \langle r_k - \tau_k T r_k, r_k - \tau_k Tr_k \rangle = \langle r_k, r_k \rangle - \tau_k \langle r_k, Tr_k \rangle - \tau_k \langle Tr_k, r_k \rangle + \tau^2_k \langle T r_k, T r_k \rangle = \langle r_k, r_k \rangle - 2 \tau_k \langle r_k, Tr_k \rangle + \tau^2_k \langle T r_k, T r_k \rangle $$ So, we have to minimize the following function over $\tau_k$ : $$ f(\tau_k) = \langle r_k, r_k \rangle - 2 \tau_k \langle r_k, Tr_k \rangle + \tau^2_k \langle T r_k, T r_k \rangle $$ It can be easily done using derivative of $f(\tau_k)$ over $\tau_k$: $$ f'(\tau_k) = -2 \langle r_k, T r_k \rangle + 2 \tau_k \langle T r_k, Tr_k \rangle = 0 $$ Hence on every step we have to take the following $\tau_k$: $$ \tau_k = \frac{\langle r_k, T r_k \rangle}{\langle Tr_k, Tr_k\rangle} $$
iterative
that outputs accuracy –– a numpy array of relative errors $\big\{\frac{\|x_{k+1} - x\|_2}{\|x\|_2}\big\}$ after num_iter iterations using $\tau_k$ from the previous task. Note: The only loop you are allowed to use here is a loop for $k$.
In [1]:
# INPUT: N - int (positive), alpha - float (positive), num_iter - integer (positive),
# y - np.array (shape: Nx1, convolved signal), s - np.array (shape: Nx1, original signal)
# OUTPUT: rel_error - np.array size (num_iter x 1)
def iterative(N, num_iter, y, s, alpha):
T = gen_toeplitz(N, alpha)
x_k = np.zeros(s.shape)
rel_error = []
for k in range(num_iter):
r_k = T @ x_k - y
tau_k = np.dot(r_k, T @ r_k) / np.dot(T @ r_k, T @ r_k)
x_k = x_k - tau_k * (T @ x_k - y)
rel_error.append(np.linalg.norm(x_k - s)/np.linalg.norm(s))
return rel_error
num_iter=1000, x=s[::20] and do a convergence plot for $\alpha = \frac{1}{2}$ and $\alpha = \frac{1}{5}$.
In [566]:
num_iter = 1000
x = audio[::5]
N = len(x)
x_conv_alpha1 = convolution(x, N, 0.5)
x_conv_alpha2 = convolution(x, N, 0.2)
conv_err_alpha1 = iterative(len(x), num_iter, x_conv_alpha1, x, 0.5)
conv_err_alpha2 = iterative(len(x), num_iter, x_conv_alpha2, x, 0.2)
plt.figure(figsize=(15, 7))
plt.plot(conv_err_alpha1, label=r"$\alpha = 0.5$")
plt.plot(conv_err_alpha2, label=r"$\alpha = 0.2$")
plt.legend(loc="upper right")
plt.ylabel("Relative error")
plt.xlabel("Iteration number")
plt.title(r'Convergence plot for $\alpha = 0.1$ and $\alpha = 0.2$')
plt.show()
x=s[::20], num_iter=1000 and $\alpha=\frac{1}{5}$. Explain what happens with the convergence if you add small random noise of amplitude $10^{-3}\max(x)$ to $y$. The answer to this question should be an explanation supported by plots and/or tables.
In [571]:
num_iter = 1000
#Before we have already cut every forth signal
x = audio[::5]
N = len(x)
noise = np.random.random(x.shape) * 0.001 * np.max(x)
x_conv_alpha1 = convolution(x, N, 0.2)
x_conv_noise = x_conv_alpha1 + noise
conv_err_alpha1 = iterative(len(x), num_iter, x_conv_alpha1, x, 0.2)
conv_err_noise = iterative(len(x), num_iter, x_conv_noise, x, 0.2)
plt.figure(figsize=(15, 7))
plt.plot(conv_err_alpha1, label="Deconvolution")
plt.plot(conv_err_noise, label="Deconvolution + noise")
plt.legend(loc="upper left")
plt.ylabel("Relative error")
plt.xlabel("Iteration number")
plt.title(r'Convergence plot for deconvolution w/ and w/o noise ($\alpha = 0.2$)')
plt.show()
1.
2.
3.
(3 pts) Differentiate with respect to $A$ the function $$ f(A) = \mathrm{sin}(x^\top A B C D x), $$ where $x$ is a vector and $A, B, C, D$ are square matrices.
(7 pts) Differentiate with respect to $y, A, X$ the function $$f(y, A, X) = \mathrm{tr}(\mathrm{diag}(y) A X),$$ where $y \in \mathbb{R}^n$ and $A, X \in \mathbb{R}^{n \times n}$. Here
1.
2.
3.
y_i, & \text{if}\ i = j \\
0, & \text{otherwise}
In [85]:
def naive_multiplication(A, B):
"""
Implement naive matrix multiplication with explicit for cycles
Parameters: Matrices A, B
Returns: Matrix C = AB
"""
sizeA = A.shape
sizeB = B.shape
C = np.zeros((sizeA[0], sizeB[1]))
for i in range(sizeA[0]):
for j in range(sizeB[1]):
for k in range(sizeA[1]):
C[i,j] += A[i,k] * B[k,j]
return C
2. (7 pts) Implement the Strassen algorithm.
In [308]:
def strassen(A, B):
"""
Implement Strassen algorithm for matrix multiplication
Parameters: Matrices A, B
Returns: Matrix C = AB
"""
n = A.shape[0]
if (n == 1):
return A[0,0] * B[0,0]
# Here we suppose that input matrices have the same size and ONLY square
if ((n & (n - 1)) != 0):
padA = int(2 ** np.ceil(np.log2(n)))
A = np.pad(A, ((0, padA - n),(0, padA - n)), 'constant')
B = np.pad(B, ((0, padA - n),(0, padA - n)), 'constant')
n = padA
n2 = int(n/2)
A11, A12, A22, A21 = A[:n2:, :n2:], A[:n2:, n2::], A[n2::, n2::], A[n2::, :n2:]
B11, B12, B22, B21 = B[:n2:, :n2:], B[:n2:, n2::], B[n2::, n2::], B[n2::, :n2:]
m1 = strassen(A11 + A22, B11 + B22)
m2 = strassen(A21 + A22, B11)
m3 = strassen(A11, B12 - B22)
m4 = strassen(A22, B21 - B11)
m5 = strassen(A11 + A12, B22)
m6 = strassen(A21 - A11, B11 + B12)
m7 = strassen(A12 - A22, B21 + B22)
if (isinstance(m1, np.ndarray)):
C = np.bmat([[m1 + m4 - m5 + m7, m3 + m5], [m2 + m4, m1 - m2 + m3 + m6]])
else:
C = np.matrix([[m1 + m4 - m5 + m7, m3 + m5], [m2 + m4, m1 - m2 + m3 + m6]])
return C
3. (5 pts) Compare three approaches: naive multiplication, Strassen algorithm and standard NumPy function.
Provide a plot in log-scale of dependence between the matrix size and the runtime of multiplication. You will have three lines, do not forget to add legend, axis labels and other attributes (see our requirements)
Consider the matrix size in the range of 100 to 700 with step 100, e.g. $n=100, 200,\ldots, 700$.
Justify the results theoretically (e.g., use the known formulas for total multiplication complexity of naive and Strassen algorithms).
In [307]:
import time
matrix_size = [x for x in range(100, 800, 100)]
naive_impl = []
strassen_impl = []
numpy_impl = []
for n in matrix_size:
print(n)
A = np.random.random((n, n))
B = np.random.random((n, n))
start = time.time()
naive_multiplication(A,B)
end = time.time()
naive_impl.append(end - start)
start = time.time()
strassen(A,B)
end = time.time()
strassen_impl.append(end - start)
start = time.time()
np.matmul(A,B)
end = time.time()
numpy_impl.append(end - start)
In [393]:
plt.figure(figsize=(15,7))
plt.loglog(matrix_size, naive_impl, label='Naive algorithm')
plt.loglog(matrix_size, strassen_impl, label='Strassen algorithm')
plt.loglog(matrix_size, numpy_impl, label='Numpy algorithm')
plt.ylabel(r"$t/t_0$")
plt.xlabel("Square matrix size, items")
plt.title("Different multiplication algorithm performance")
plt.legend(loc='upper left')
plt.show()
Regardless of theoretical expectations, Strassen algorithm works slower than naive implementation and there is some reason behind it. First of all, it is hard to achieve some accelerations using python directly because it is not fast and causes big overhead to process every line of code (as an example, numpy efficient multiplication function falls back to C/C++ function which calculates result). Second, since on every step a big amount of memory are allocated there is some extra time consumed. And finally it is because of recursion, which is also quite heavy procedure in sense of memory allocation. Interestingly, one might observe some step on the plot. Steps are reasonable because all matrices are padded to the closest power of 2 and only then multiplied.
1. (2 pts) Compute the singular values of some predownloaded image (via the code provided below) and plot them. Do not forget to use logarithmic scale.
In [360]:
%matplotlib inline
import matplotlib.pyplot as plt
from PIL import Image, ImageDraw
import requests
import numpy as np
url = 'https://people.sc.fsu.edu/~jburkardt/data/png/lena.png'
face_raw = Image.open(requests.get(url, stream=True).raw)
face = np.array(face_raw).astype(np.uint8)
plt.imshow(face_raw)
plt.xticks(())
plt.yticks(())
plt.title('Original Picture')
plt.show()
plt.figure(figsize=(15,7))
u, s, vh = np.linalg.svd(face)
plt.semilogy(s/s[0])
plt.xlabel("Number of sing. value")
plt.ylabel(r"$\sigma_i/\sigma_0$")
plt.title("Distribution of singular value for three color channels")
plt.show()
2. (3 pts) Complete a function compress, that performs SVD and truncates it (using $k$ singular values/vectors). See the prototype below.
Note, that in colourful case you have to split your image to channels and work with matrices corresponding to different channels separately.
Plot approximate reconstructed image $M_\varepsilon$ of your favorite image such that $rank(M_\varepsilon) = 5, 20, 50$ using plt.subplots.
In [485]:
def compress(image, k):
"""
Perform svd decomposition and truncate it (using k singular values/vectors)
Parameters:
image (np.array): input image (probably, colourful)
k (int): approximation rank
--------
Returns:
reconst_matrix (np.array): reconstructed matrix (tensor in colourful case)
s (np.array): array of singular values
"""
u, s, v = np.linalg.svd(image)
diag_s = s.transpose() * np.identity(len(s))
u_shr = u[:,:k]
v_shr = v[:k,:]
reconst_matrix = u_shr @ diag_s[:k,:k] @ v_shr
return reconst_matrix, s[:k]
In [486]:
plt.figure(figsize=(30,15))
for i,k in enumerate([5, 20, 50]):
plt.subplot(1, 3, i + 1)
img, sv = compress(face, k)
plt.xticks(())
plt.yticks(())
plt.imshow(img, cmap='gray')
plt.title('Number of singular values = {}'.format(k), fontsize=26)
3. (3 pts) Plot the following two figures for your favorite picture
In [477]:
rank = [x for x in range(1, 200, 2)]
rel_error = []
compr_rate = []
for r in rank:
img, sv = compress(face, r)
rel_error.append(np.linalg.norm(img - face, 'fro')/np.linalg.norm(face, 'fro'))
compr_rate.append((face.shape[0] * r + face.shape[1] * r + r) / (face.shape[0] * face.shape[1]))
plt.figure(figsize=(14,7))
plt.xlabel("Rank of approximation")
plt.ylabel("Relative error")
plt.title("Relative error of approximation")
plt.plot(rank, rel_error)
plt.show()
plt.figure(figsize=(14,7))
plt.xlabel("Rank of approximation")
plt.ylabel("Compression rate")
plt.title("Compression rate over approximation rank")
plt.plot(rank, compr_rate)
plt.show()
4. (2 pts) Consider the following two pictures. Compute their approximations (with the same rank, or relative error). What do you see? Explain results.
In [517]:
url1 = 'http://sk.ru/resized-image.ashx/__size/550x0/__key/communityserver-blogs-components-weblogfiles/00-00-00-60-11/skoltech1.jpg'
url2 = 'http://www.simpsoncrazy.com/content/characters/poster/bottom-right.jpg'
image_raw1 = Image.open(requests.get(url1, stream=True).raw)
image_raw2 = Image.open(requests.get(url2, stream=True).raw)
image1 = np.array(image_raw1).astype(np.uint8)
image2 = np.array(image_raw2).astype(np.uint8)
plt.figure(figsize=(18, 6))
plt.subplot(1,2,1)
plt.imshow(image_raw1)
plt.title('One Picture')
plt.xticks(())
plt.yticks(())
plt.subplot(1,2,2)
plt.imshow(image_raw2)
plt.title('Another Picture')
plt.xticks(())
plt.yticks(())
plt.show()
In [531]:
# Your code is here
rank = 50
dim = image1.ndim
image_zip1 = np.zeros(image1.shape)
image_zip2 = np.zeros(image2.shape)
for i in range(dim):
image_zip1[:,:,i], sing_val1 = compress(image1[:,:,i], rank)
image_zip2[:,:,i], sing_val2 = compress(image2[:,:,i], rank)
plt.figure(figsize=(18, 6))
plt.subplot(1,2,1)
plt.imshow(Image.fromarray(image_zip1.astype(np.uint8), 'RGB'))
plt.title('One Picture')
plt.xticks(())
plt.yticks(())
plt.subplot(1,2,2)
plt.imshow(Image.fromarray(image_zip2.astype(np.uint8), 'RGB'))
plt.title('Another Picture')
plt.xticks(())
plt.yticks(())
plt.show()
To figure out why the second image is much more distorted than the second one (basically, in the first case most of objects are well recognizable compared to the second case) let's plot their singular values.
In [521]:
_, sing_val1, _ = np.linalg.svd(image1[:,:,0])
_, sing_val2, _ = np.linalg.svd(image2[:,:,0])
plt.figure(figsize=(18, 6))
plt.subplot(1,2,1)
plt.semilogy(sing_val1/sing_val1[0])
plt.ylabel(r"$\sigma_i/\sigma_0$")
plt.xlabel("Number of sing. value")
plt.title('One Picture, sing. values')
plt.subplot(1,2,2)
plt.semilogy(sing_val2/sing_val2[0])
plt.ylabel(r"$\sigma_i/\sigma_0$")
plt.xlabel("Number of sing. value")
plt.title('Another Picture, sing. values')
plt.show()
As we know the more singular values we leave the better quality of image will be, but at the same time, the greater singular values we have the more ones contribute to the final image quality (it is obvious from definition of SVD). Keeping that in mind, we can observe that in the first image we have up to 60 singular values from $\sigma_0$ to $0.01\sigma_0$, so the contribution of the rest is less important. But in the second case we have up to 150 singular values within the same range. Hence, cutting the same number of singular values from both images is wrong on its head in terms of similarity of approximated images. Taking singular values within the same relative range will be more fair, so let's check it.
In [532]:
# Your code is here
dim = image1.ndim
image_zip1 = np.zeros(image1.shape)
image_zip2 = np.zeros(image2.shape)
for i in range(dim):
image_zip1[:,:,i], sing_val1 = compress(image1[:,:,i], 60)
image_zip2[:,:,i], sing_val2 = compress(image2[:,:,i], 150)
plt.figure(figsize=(18, 6))
plt.subplot(1,2,1)
plt.imshow(Image.fromarray(image_zip1.astype(np.uint8), 'RGB'))
plt.title('One Picture')
plt.xticks(())
plt.yticks(())
plt.subplot(1,2,2)
plt.imshow(Image.fromarray(image_zip2.astype(np.uint8), 'RGB'))
plt.title('Another Picture')
plt.xticks(())
plt.yticks(())
plt.show()
So both images are still well-recognizable.
The norm is called absolute if $\|x\|=\| \lvert x \lvert \|$ holds for any vector $x$, where $x=(x_1,\dots,x_n)^T$ and $\lvert x \lvert = (\lvert x_1 \lvert,\dots, \lvert x_n \lvert)^T$. Give an example of a norm which is not absolute.
Write a function ranks_HOSVD(A, eps)
that calculates Tucker ranks of a d-dimensional tensor $A$ using High-Order SVD (HOSVD) algorithm, where eps is the relative accuracy in the Frobenius norm between the approximated and the initial tensors. Details can be found here on Figure 4.3.
def ranks_HOSVD(A, eps):
return r #r should be a tuple of ranks r = (r1, r2, ..., rd)
In [ ]: