import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set(font_scale=2)
plt.style.use('seaborn-whitegrid')
from scipy import sparse
from PIL import Image
img = Image.open('ssn.png')
xmat = (255 - np.asarray(img).max(axis=2))/255
print(xmat.shape)
print(xmat.min(),xmat.max())
plt.imshow(xmat);
x
¶x = xmat.flatten()
x.shape
More about this blur matrix on a later MP...
imat, jmat = np.meshgrid(np.arange(xmat.shape[0]), np.arange(xmat.shape[1]), indexing='ij')
ivec = np.atleast_2d(imat.flatten())
jvec = np.atleast_2d(jmat.flatten())
A = np.fmax(0, 1 - np.sqrt((ivec.T - ivec)**2 + (jvec.T - jvec)**2)/5)
A /= A.sum(axis=1)
b = A @ x
b2D=b.reshape(xmat.shape)
plt.imshow(b2D)
import scipy.linalg as sla
P, L, U = sla.lu(A)
If $Ax = P L U x = b$, then there are two steps:
y = sla.solve_triangular(L, np.dot(P.T, b), lower=True)
x_solve = sla.solve_triangular(U, y)
plt.imshow(x_solve.reshape(xmat.shape))
np.linalg.solve
?¶x_solve_2 = la.solve(A,b)
plt.imshow(x_solve_2.reshape(xmat.shape))
Suppose you have many social security numbers to un-blur. You should factorize your blur matrix only once and then perform several triangular solves.
Let's time:
%timeit sla.solve(A, b)
%timeit P, L, U = sla.lu(A)
%timeit sla.solve_triangular(U, y)
plt.figure()
plt.subplot(131)
plt.spy(A); plt.xticks([]); plt.yticks([]);
plt.subplot(132)
plt.spy(L); plt.xticks([]); plt.yticks([]);
plt.subplot(133)
plt.spy(U); plt.xticks([]); plt.yticks([]);
A_csr = sparse.csr_matrix(A)
%timeit sparse.linalg.spsolve(A_csr,b)
x_solve_3 = sparse.linalg.spsolve(A_csr,b)
plt.imshow(x_solve_3.reshape(xmat.shape))
b_noisy = b + 1e-2 * np.random.rand(b.size)
y = sla.solve_triangular(L, np.dot(P.T, b_noisy), lower=True)
x_solve = sla.solve_triangular(U, y)
plt.imshow(x_solve.reshape(xmat.shape))