In [5]:
%matplotlib inline
import scipy as sc
import numpy as np
#import pandas as pd
#import matplotlib as mpl
import matplotlib.pylab as plt
In [7]:
A = 1.5*np.mat(np.random.randn(2,2))
fig = plt.figure(figsize=(6,6))
p = 2
q = 2
X = norm_ball(p=q)
#X = np.matrix(2*np.random.rand(2,50)-1)
Y = A*X
plt.plot(Y[0,:].T, Y[1,:].T,'b-')
#plt.plot(2*Y[0,:].T, 2*Y[1,:].T,'b-')
#plt.plot(Y[0,:].T, Y[1,:].T,'b.')
for a in np.arange(1,6,1):
X = norm_ball(p=p)
Y = a*X
#print(X)
plt.plot(Y[0,:].T, Y[1,:].T,'k:',alpha=0.5)
#ax.set_xlim((-2,2))
#ax.set_ylim((-2,2))
#plt.plot(X[0,:].tolist())
ax = fig.gca()
ax.set_xlim((-3,3))
ax.set_ylim((-3,3))
ax.axis('equal')
plt.show()
Given a $p \geq 1$, a $p$-norm ball is the set of points $$ B_p = \{ x : \| x \|_p \leq 1 \} $$ Design an algorithm to sample uniform points in a given norm-ball $B_p$, where $x \in \mathbb{R}^m$.
Let $C = AB$. Hence, $c_{i,j} = a_i^\top b_j$ where $a_i^\top$ is the $i$'th row of $A$ and $b_j$ is the $j$'th column of $B$
\begin{eqnarray} \| A B \|_F^2 & = & \sum_{i} \sum_j |c_{i,j} |^2 \\ & = & \sum_{i} \sum_j |a^\top_{i} b_j |^2 \\ & \leq & \sum_{i} \sum_j \|a_i\|^2 \|b_j\|^2 \\ & = & \left(\sum_{i} \|a_i\|^2 \right) \left(\sum_j \|b_j\|^2 \right)\\ & = & \| A \|_F^2 \| B \|_F^2 \\ \end{eqnarray}Hence $$ \| A B \|_F \leq \| A \|_F \| B \|_F $$
QR for polynomials
In [158]:
x = np.mat(np.arange(-1.,1.,0.01)).T
N = len(x)
degree = 10
#A = np.hstack((np.power(x,0), np.power(x,1), np.power(x,2)))
B = np.hstack((np.power(x,i) for i in range(degree+1)))
B = B[:,np.random.permutation(degree+1)]
#B = np.hstack((np.power(x,i) for i in range(degree+1)))
#B = np.random.randn(N,degree+1)
A, R, = np.linalg.qr(B)
# Append an extra basis for outliers
K = A.shape[1]
plt.figure(figsize=(5,20))
#plt.show()
idx = range(0,N,20)
idx.append(N-1)
for i in range(K):
plt.subplot(K,1,i+1)
plt.stem(x[idx],A[idx,i])
plt.plot(x,A[:,i])
plt.gcf().gca().set_xlim([-1.2, 1.2])
plt.gcf().gca().set_ylim([-0.3, 0.3])
plt.gcf().gca().axis('off')
plt.show()
A reflector is like a projection onto some hyperplane but we go twice further. $$ H = I - 2P $$ The projection is onto the hyperplane orthogonal to $v$ $$ P = I - \frac{v}{\|v\|} \frac{v^\top}{\|v\|} = I - \frac{v v^\top}{v^\top v} $$ where, $v$ is the difference vector $ v = \|x\|e_1 - x $.
$$ H = I - 2I + 2qq' = -I + 2qq' = -(I - 2qq') $$If $P$ is an orthogonal projector, the reflector $I-2P$ is an orthonormal transformation as we have $P^\top P = P^2 = P$ for an orthogonal projector. $(I-2P)^\top (I-2P) = I - 2P^\top -2P + 4P^\top P = I$
In [161]:
import numpy as np
import numpy.linalg as la
np.set_printoptions(precision=3, suppress=True)
def house(x):
'''Householder Reflection'''
m = len(x)
e = np.mat(np.zeros((m,1)))
e[0] = 1
v = la.norm(x)*e - x
z = v/la.norm(v)
H = np.eye(m) - 2*z*z.T
return z,H
x = np.mat('[1;2;3;4]')
z,H = house(x)
print(x)
print(H*x)
print(x - 2*z*(z.T*x))
In [164]:
np.set_printoptions(precision=3, suppress=True)
m = 5
n = 7
A = np.mat(np.random.randn(m, n))
R = A.copy()
Q = np.eye(m)
print('A')
print(A)
for i in range(np.min((n,m))-1):
z,H = house(R[i:,i])
R[i:,i:] = H*R[i:,i:]
Q[:,i:] = Q[:,i:]*H.T
#print('q =')
#print(q)
#print('Q =')
#print(Q)
print('R =')
print(R)
print('--'*4)
print(Q*R)
print('--'*4)
print(A)
We have explicitely constructed $H$ (or $P$) but this is not needed. It is sufficient to store the Householder vectors $v$ to construct/implement transformations by $H$.
In [166]:
def random_house(m, k=None):
'''Generate a random Householder Reflection'''
if k is None:
k = m
e = np.mat(np.zeros((k,1)))
e[0] = 1
x = np.mat(np.random.randn(k,1))
v = la.norm(x)*e - x
q = v/la.norm(v)
q = np.vstack((np.zeros((m-k,1)), q))
H = np.eye(m) - 2*q*q.T
return q,H
q, H = random_house(5)
print(q)
np.linalg.matrix_rank(H)
Out[166]:
In [185]:
N = 5
A = [random_house(N,k)[1] for k in range(N,1,-1)]
qq = np.hstack([random_house(N,k)[0] for k in range(N,1,-1)])
qq
Out[185]:
In [155]:
A[:,0].T*A[:,1]
Out[155]: