Problem set 3 (90 pts)

Important note: the template for your solution filename is Name_Surname_PS3.ipynb

For this problem set we do not run the bot, so try to debug your solutions with your own simple tests

Problem 1 (20 pts)

  • (5 pts) Prove that $\mathrm{vec}(AXB) = (B^\top \otimes A)\, \mathrm{vec}(X)$ if $\mathrm{vec}(X)$ is a columnwise reshape of a matrix into a long vector. What does it change if the reshape is rowwise?

Note:

  1. To make a columnwise reshape in Python one should use np.reshape(X, order='f'), where the string 'f' stands for the Fortran ordering.
  2. If $\mathrm{vec}(X)$ is a rowwise reshape, $$\mathrm{vec}(AXB)=(A \otimes B^\top) \mathrm{vec}(X).$$

  • (2 pts) What is the complexity of a naive computation of $(A \otimes B) x$? Show how it can be reduced.

  • (3 pts) Let matrices $A$ and $B$ have eigendecompositions $A = S_A\Lambda_A S_A^{-1}$ and $B = S_B\Lambda_B S^{-1}_B$. Find eigenvectors and eigenvalues of the matrix $A\otimes I + I \otimes B$.

  • (10 pts) Let $A = \mathrm{diag}\left(\frac{1}{1000},\frac{2}{1000},\dots \frac{999}{1000}, 1, 1000 \right)$. Estimate analytically the number of iterations required to solve linear system with $A$ with the relative accuracy $10^{-4}$ using
    • Richardson iteration with the optimal choice of parameter (use $2$-norm)
    • Chebyshev iteration (use $2$-norm)
    • Conjugate gradient method (use $A$-norm).
In [1]:
# Your solution is here

Problem 2 (40 pts)

Spectral graph partitioning and inverse iteration

Given connected graph $G$ and its corresponding graph Laplacian matrix $L = D - A$ with eigenvalues $0=\lambda_1, \lambda_2, ..., \lambda_n$, where $D$ is its degree matrix and $A$ is its adjacency matrix, Fiedler vector is an eignevector correspondng to the second smallest eigenvalue $\lambda_2$ of $L$. Fiedler vector can be used for graph partitioning: positive values correspond to the one part of a graph and negative values to another.

Inverse power method (15 pts)

To find the Fiedler vector we will use the inverse iteration with adaptive shifts (Rayleigh quotient iteration).

  • (5 pts) Write down the orthoprojection matrix on the space orthogonal to the eigenvector of $L$, corresponding to the eigenvalue $0$ and prove (analytically) that it is indeed an orthoprojection.

  • (5 pts) Implement the spectral partitioning as the function partition:

In [3]:
# INPUT:
# A - adjacency matrix (scipy.sparse.csr_matrix)
# num_iter_fix - number of iterations with fixed shift (int)
# shift - (float number)
# num_iter_adapt - number of iterations with adaptive shift (int) -- Rayleigh quotient iteration steps
# x0 - initial guess (1D numpy.ndarray)
# OUTPUT:
# x - normalized Fiedler vector (1D numpy.ndarray)
# eigs - eigenvalue estimations at each step (1D numpy.ndarray)
# eps - relative tolerance (float)
def partition(A, shift, num_iter_fix, num_iter_adapt, x0, eps):
    x = x0
    eigs = np.array([0])
    return x, eigs

Algorithm must halt before num_iter_fix + num_iter_adapt iterations if the following condition is satisfied $$ \boxed{\|\lambda_k - \lambda_{k-1}\|_2 / \|\lambda_k\|_2 \leq \varepsilon} \text{ at some step } k.$$

Do not forget to use the orthogonal projection from above in the iterative process to get the correct eigenvector. It is also a good idea to use shift=0 before the adaptive stragy is used. This, however, is not possible since the matrix $L$ is singular, and sparse decompositions in scipy does not work in this case. Therefore, we first use a very small shift instead.

  • (3 pts) Generate a random lollipop_graph using networkx library and find its partition. Draw this graph with vertices colored according to the partition.

  • (2 pts) Start the method with a random initial guess x0, set num_iter_fix=0 and comment why the method can converge to a wrong eigenvalue.

Spectral graph properties (15 pts)

  • (5 pts) Prove that multiplicity of the eigenvalue $0$ in the spectrum of the graphs Laplacian is the number of its connected components.
  • (10 pts) The second-smallest eigenvalue of $L(G)$, $\lambda_2(L(G))$, is often called the algebraic connectivity of the graph $G$. A basic intuition behind the use of this term is that a graph with a higher algebraic connectivity typically has more edges, and can therefore be thought of as being “more connected”.
    To check this statement, create few graphs with equal number of vertices using networkx, one of them should be $C_{30}$ - simple cyclic graph, and one of them should be $K_{30}$ - complete graph. (You also can change the number of vertices if it makes sense for your experiments, but do not make it trivially small).
    • Find the algebraic connectivity for the each graph using inverse iteration.
    • Plot the dependency $\lambda_2(G_i)$ on $|E_i|$.
    • Draw a partition for a chosen graph from the generated set.
    • Comment on the results.

Image bipartition (10 pts)

Let us deal here with a graph constructed from a binarized image. Consider the rule, that graph vertices are only pixels with $1$, and each vertex can have no more than $8$ connected vertices (pixel neighbours), $\textit{i.e}$ graph degree is limited by 8.

  • (3 pts) Find an image with minimal size equal to $(256, 256)$ and binarize it such that graph built on black pixels has exactly $1$ connected component.
  • (5 pts) Write a function that constructs sparse adjacency matrix from the binarized image, taking into account the rule from above.
  • (2 pts) Find the partition of the resulting graph and draw the image in accordance with partition.
In [2]:
# Your solution is here

Problem 3 (30 pts)

Say hi to the drone

You received a radar-made air scan data of a terrorist hideout made from a heavy-class surveillance drone. Unfortunately, it was made with an old-fashioned radar, so the picture is convolved with the diffractive pattern. You need to deconvolve the picture to recover the building plan.

In [3]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hankel2
radiointel = np.load('../files/radiointel.npy')
plt.subplot(1,2,1)
plt.imshow( np.abs(radiointel) )
plt.title('Intensity')
plt.subplot(1,2,2)
plt.imshow( np.angle(radiointel), cmap='bwr' )
plt.title('Phase')
plt.show()

In this problem you asked to use using FFT-matvec and make the convolution operator for the picture of the size $N\times N$, where $N=300$ with the following kernel (2D Helmholtz scattering): $$ T_{\overline{i_1 j_1}, \overline{i_2 j_2} } \equiv eG_{i_1-j_1,i_2-j_2} = \frac{-1j}{4} H^{(2)}_0 \left( k_0 \cdot \Delta r_{ \overline{i_1 j_1}, \overline{i_2 j_2} } \right), \quad i_1,j_1, i_2, j_2 = 0,\dots, N-1 $$

except when both $i_1=i_2$ and $j_1 = j_2$.

In that case set $$T_{i_1=i_2, j_1=j_2} = 0$$.

Here $1j$ is the imaginary unit, $H^{(2)}_0(x)$ - (complex-valued) Hankel function of the second kind of the order 0. See 'scipy.special.hankel2'.

$$ \Delta r_{ \overline{i_1 j_1}, \overline{i_2 j_2} } = h \sqrt{ (i_1-i_2)^2 + (j_1-j_2)^2 } $$ $$ h = \frac{1}{N-1}$$ $$k_0 = 50.0$$

Tasks:

  1. (5 pts) Create the complex-valued kernel $eG$ ($2N-1 \times 2N-1$)-sized matrix according with the instructions above. Note that at the point where $\Delta r=0$ value of $eG$ should be manually zet to zero. Store in the variable eG. Plot the eG.real of it with plt.imshow

  2. (5 pts) Write function Gx that calculates matvec of $T$ by a given vector $x$. Make sure all calculations and arrays are in dtype=np.complex64. Hint: matvex with a delta function in pl

  3. (3 pts) What is the complexity of one matvec?

  4. (2 pts) Use scipy.sparse.linalg.LinearOperator to create an object that has attribute .dot() (this object will be further used in the iterative process). Note that .dot() input and output must be 1D vectors, so do not forget to use reshape.
  5. (15 pts) Write a function that takes an appropriate Krylov method(s) and solves linear system $Gx=b$ to deconvolve radiointel. The result should be binary mask array (real, integer, of 0s and 1s) of the plane of the building. Make sure it converged sufficiently and you did the post-processing properly. Plot the result as an image.

Note: You can use standart fft and ifft from e.g. numpy.fft

1. Kernel (5 pts)

In [ ]:
k0 = #
N = #

def make_eG(k0, N):
    # INPUT:  
    # k0 #dtype = float
    # N #dtype = int
    
    # OUTPUT:
    # np.array, shape = (2N-1, 2N-1), dtype = np.complex64
    return eG

eG = make_eG(k0=k0, N=N)

plt.imshow(eG.real)

2. Matvec (5 pts)

In [ ]:
def Gx(x, eG):
    # input:  
    # x, np.array, shape=(N**2, ), dtype = np.complex64
    # eG, np.array, shape=(2N-1, 2N-1), dtype = np.complex64
    # output:
    # matvec, np.array, shape = (N**2, ), dtype = np.complex64

    return matvec

3. Complexity (3 pts)

Big-O complexity of one matvec operation is ... It can be shown...

4. LinearOperator (2 pts)

In [ ]:
L_Gx =

5. Reconstruction (15pts)

In [ ]:
def normalize(mask): #proper normalization to binary mask
    mask = np.clip(mask, a_min=0, a_max=1)
    mask = np.round(mask)
    mask = np.asarray(mask, dtype=int)
    return mask

errs=[]
def callback(err): #callback function to store the history of convergence
    global errs
    errs.append(err)
    return 

mask = #some_solver(, , , callback = callback)

plt.figure()
plt.imshow( normalize(mask) , cmap='binary')
plt.title('Restored building plan')
plt.colorbar()

plt.figure()
plt.semilogy(errs)
plt.title('Convergence')