Note:
np.reshape(X, order='f'), where the string 'f' stands for the Fortran ordering. But we might get better analytical results if one tries to derive $k$ from the following formula: $$ \tau_k = \frac{\|x_k - x_* \|_A^2}{\| x_* \|^2_A} \leq \underset{deg(q) \leq k, q(0) = 1}{\mbox{inf}} \Big( \underset{i = 1...n}{\mbox{max}} q(\lambda_i)^2 \Big) $$ On each step there will be polynomial $q$ such that deg$q$ = n and all roots are eigenvalues of $A$, so $\tau_k$ = 0. Such polynomial will be available from $n_{th}$ step, so we need a number of iterations equal to the number of eigenvalues $A$ has or 1001, so $\boxed{k \geq 1001}$
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.
To find the Fiedler vector we will use the inverse iteration with adaptive shifts (Rayleigh quotient iteration).
partition:
In [7]:
import numpy as np
import scipy as sp
from scipy import sparse
from scipy.sparse import linalg
import networkx as nx
from networkx.linalg.algebraicconnectivity import fiedler_vector
import matplotlib.pyplot as plt
In [248]:
# 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 operator_P(x):
return x - np.ones_like(x) * np.sum(x) / x.shape[0]
def is_small(eigs, eps):
if (abs(eigs[-2] - eigs[-1]) / eigs[-1] <= eps):
return True
def partition(A, shift, num_iter_fix, num_iter_adapt, x0, eps):
x_k = x0 / np.linalg.norm(x0)
eigs = np.array([0])
L, I = nx.laplacian_matrix(nx.from_scipy_sparse_matrix(A)), sp.sparse.identity(x_k.shape[0])
for _ in range(num_iter_fix):
x_k = operator_P(x_k)
x_k = sp.sparse.linalg.spsolve(L - shift * I, x_k)
x_k /= np.linalg.norm(x_k, 2)
eigs = np.append(eigs, np.dot(x_k, L.dot(x_k)))
if is_small(eigs, eps):
return x_k, eigs
for _ in range(num_iter_adapt):
x_k = operator_P(x_k)
x_k = sp.sparse.linalg.spsolve(L - np.dot(x_k, L.dot(x_k)) * I, x_k)
x_k /= np.linalg.norm(x_k, 2)
eigs = np.append(eigs, np.dot(x_k, L.dot(x_k)))
if is_small(eigs, eps):
return x_k, eigs
return x_k, 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.
lollipop_graph using networkx library and find its partition. Draw this graph with vertices colored according to the partition.
In [249]:
m, n = (10, 20)
G = nx.lollipop_graph(m, n)
A = nx.adjacency_matrix(G)
x0 = np.random.random((A.shape[0],)).astype(np.float64)
eigs_builtin = np.sort(np.linalg.eigh(nx.laplacian_matrix(G).todense())[0])
x, eigs = partition(A, 0.01, 10, 100, x0, 0.00001)
print("Are second smallest the same?", np.allclose(eigs_builtin[1], eigs[-1]))
In [250]:
pos = nx.spring_layout(G)
plt.figure(figsize=(14,7))
nx.draw_networkx(G, pos, node_color=np.sign(x))
x0, set num_iter_fix=0 and comment why the method can converge to a wrong eigenvalue.
In [251]:
x0 = np.random.random((A.shape[0],)).astype(np.float32)
eigs_builtin = np.sort(np.linalg.eigh(nx.laplacian_matrix(G).todense())[0])
x, eigs = partition(A, 1e-3, 0, 5, x0, 1e-5)
print("Are second smallest the same?", np.allclose(eigs_builtin[1], eigs[-2]))
print(eigs_builtin[1], eigs[-1])
Fixed shift is important because w/o one it starts converging to the number equal to initial shift, which can be seen in code explicitly $x_0^T L x_0$
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).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.
In [265]:
from PIL import Image
import requests
url = "https://findicons.com/files/icons/2773/pictonic_free/256/lin_ubuntu.png"
img = np.array(Image.open(requests.get(url, stream=True).raw).convert('RGBA')).astype(np.uint8)[:,:,3]
img = 1.0 * (img < 100)
plt.imshow(img)
Out[265]:
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 [170]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import hankel2
radiointel = np.load('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$$(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
(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 pts) What is the complexity of one matvec?
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.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
In [234]:
k0 = 50.0
N = 300
def make_eG(k0, N):
# INPUT:
# k0 #dtype = float
# N #dtype = int
# OUTPUT:
# np.array, shape = (2N-1, 2N-1), dtype = np.complex64
# eG = np.zeros((2*N - 1, 2*N - 1), dtype=np.complex64)
row = np.square(np.linspace(-N + 1, N - 1, 2*N - 1))
column = np.square(np.linspace(-N + 1, N - 1, 2*N - 1))
row_m, col_m = np.meshgrid(row, column, indexing = 'ij')
eG = -1j * hankel2(0, k0 / (N - 1) * np.sqrt(row_m + col_m)) / 4
eG[N - 1, N - 1] = 0
return np.roll(eG.astype(np.complex64), shift=N-1, axis=(0,1))
eG = make_eG(k0=k0, N=N)
plt.imshow(eG.real)
Out[234]:
In [235]:
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
N = int(np.sqrt(x.shape[0]))
x_ready = np.pad(x.reshape((N, N)), ((0, N - 1), (0, N - 1)), 'constant')
matvec = np.fft.ifft2(np.fft.fft2(eG) * np.fft.fft2(x_ready))[0:N,0:N]
return matvec.reshape(N**2, 1).astype(np.complex64)
Big-O complexity of one matvec operation is $\mathcal{O}(N^2\log{}N)$
In [236]:
from scipy.sparse import linalg
L_Gx = linalg.LinearOperator((N ** 2, N ** 2), matvec=lambda x, eG=eG: Gx(x, eG))
In [238]:
from scipy.sparse.linalg import gmres
maxiter = 200
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, _ = gmres(L_Gx, radiointel.reshape(N**2,1), maxiter=maxiter, callback = callback)
plt.figure(figsize=(7,7))
plt.imshow(normalize(mask.real).reshape((N, N)), cmap='binary')
plt.title('Restored building plan', size=16)
plt.colorbar()
plt.figure(figsize=(14,7))
plt.semilogy(errs)
plt.grid()
plt.title('Convergence', size=16)
Out[238]:
In [ ]: