examples.ipynb

Kernel: Python 3 (Anaconda)

Examples

import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt

Data-fitting

Given a signal $f \in \mathbb{R}^m$ we want to fit a model $$f_i = \sum_{j=1}^n u_j\exp(-\beta(x_i - x_j)^2).$$ Typically, $m>n$ in such applications as we are looking for a simple explanation of the signal.

m = 51
n = 10

# nodes
x = np.linspace(0,1,m)
xi = np.linspace(0,1,n)

# sample
noise = np.random.randn(m)
f = np.cos(15*x) + .5*noise

# matrix
K = np.zeros((m,n))
for k in range(n):
    K[:,k] = np.exp(-100*(x - xi[k])**2)
    

# SVD
U, S, V = np.linalg.svd(K)

plt.subplot(131)
plt.plot(S)
plt.title('singular values')

plt.subplot(132)
plt.plot(U[:,5])
plt.title('left singular vector')

plt.subplot(133)
plt.plot(V[0,:])
plt.title('right singular vector')

Text(0.5,1,'right singular vector')

# project f on range of K

plt.subplot(121)
plt.plot(x,f,'*',x,U[:,0:n]@U[:,0:n].T@f,'o',x,np.cos(15*x))
plt.legend(('f','projected f','truth'))

plt.subplot(122)
plt.plot(x,noise,'*',x,U[:,0:n]@U[:,0:n].T@noise,'o')
plt.legend(('noise','projected noise'))

<matplotlib.legend.Legend at 0x7f10cef952b0>

# solve
u = np.linalg.lstsq(K,f)[0]

plt.plot(x,f,'*')
plt.plot(x,K@u,x,np.cos(15*x))
plt.legend(('f','fit','truth'))
plt.show

<function matplotlib.pyplot.show>

Convolution

# set up matrix
n = 100
t = np.linspace(0,1,n)
K = la.toeplitz(np.exp(-100*t**2)/(10))

# data
x = np.linspace(-1,1,n)
u = np.exp(-(10*x)**50)
f = K@u
noise = 0.1*np.random.randn(n)

# plot
plt.plot(x,u,x,f)

[<matplotlib.lines.Line2D at 0x7f10ceec62b0>,
 <matplotlib.lines.Line2D at 0x7f10ceec6438>]

# SVD
U, S, V = np.linalg.svd(K)

plt.subplot(131)
plt.semilogy(S,'*')
plt.title('singular values')

plt.subplot(132)
plt.plot(U[:,10])
plt.title('left singular vectors')

plt.subplot(133)
plt.plot(V[0,:])
plt.title('right singular vectors')

Text(0.5,1,'right singular vectors')

# noise
plt.semilogy(S,'*')
plt.semilogy(U.T@(f+noise),'*')
plt.semilogy((U.T@(f+noise))/S,'o')
plt.legend(('sigma','u_i^Tf','picard'))

<matplotlib.legend.Legend at 0x7f10ced9d780>

# regularized inverse
plt.semilogy(S)
plt.semilogy((S**2 + 1e-6)/S)
plt.semilogy(S + (1e-3/S)**20)
plt.axis([0,99,1e-16,1e16])

[0, 99, 1e-16, 1e+16]

# noise
plt.semilogy(S,'*')
plt.semilogy(U.T@(f+noise),'*')
plt.semilogy((U.T@(f+noise))*S/(S**2 + 1e-6),'o')
plt.legend(('sigma','u_i^Tf','picard'))

<matplotlib.legend.Legend at 0x7f10ceaf4c88>

# regularized inverse
k = 30
alpha = 1e-7
u1 = (V[0:k,:][e-mail](1/S[0:k])@U[:,0:k].T)@f
u2 = (V[0:n,:][e-mail](S/(S**2 + alpha))@U[:,0:n].T)@f
plt.plot(x,u1,x,u2,x,u)

[<matplotlib.lines.Line2D at 0x7f10cea0a1d0>,
 <matplotlib.lines.Line2D at 0x7f10cea0a358>,
 <matplotlib.lines.Line2D at 0x7f10cea0ab70>]

# noisy data
k = 5
alpha = 1e-1
u1 = (V[0:k,:][e-mail](1/S[0:k])@U[:,0:k].T)@(f+noise)
u2 = (V[0:n,:][e-mail](S/(S**2 + alpha))@U[:,0:n].T)@(f+noise)

plt.plot(x,u1,x,u2,x,u)

[<matplotlib.lines.Line2D at 0x7f10cea33128>,
 <matplotlib.lines.Line2D at 0x7f10cea332b0>,
 <matplotlib.lines.Line2D at 0x7f10cea33ac8>]

Tomography

n = 100
m = 10
x = np.linspace(0,1,n)
zi = [10, 12, 21, 26, 35, 47, 63, 78, 81,99]

# matrix
K = np.zeros((m,n))
for k in range(m):
    K[k,0:zi[k]] = 1

# data
u = x + .5*x*np.sin(5*x)
noise = 0.1*np.random.randn(m)
f = K@u

# SVD
U, S, V = np.linalg.svd(K)

plt.subplot(131)
plt.semilogy(S,'*')
plt.title('singular values')

plt.subplot(132)
plt.plot(U[:,0])
plt.title('left singular vectors')

plt.subplot(133)
plt.plot(V[0,:])
plt.title('right singular vectors')

Text(0.5,1,'right singular vectors')

# project ground truth on row-space
plt.plot(x,u,x,u-V[m:,:].T@V[m:,:]@u)

[<matplotlib.lines.Line2D at 0x7fc2f360e160>,
 <matplotlib.lines.Line2D at 0x7fc2f360e2e8>]

# noise
plt.semilogy(S,'*')
plt.semilogy(U.T@(f+noise),'*')

[<matplotlib.lines.Line2D at 0x7fc2f358f160>]

# minimum norm solution
u1 = (V[0:m,:][e-mail](1/S[0:m])@U[:,0:m].T)@(f+noise)

plt.plot(x,u,x,u1)

[<matplotlib.lines.Line2D at 0x7fc2f34aa6d8>,
 <matplotlib.lines.Line2D at 0x7fc2f34aa860>]

# many solutions!
u2 = u1 + V[m:,:][e-mail](n-m)
plt.subplot(121)
plt.plot(x,u1,x,u2)
plt.subplot(122)
plt.plot(zi,K@u1,zi,K@u2)

[<matplotlib.lines.Line2D at 0x7fc2f33163c8>,
 <matplotlib.lines.Line2D at 0x7fc2f32b94e0>]