In [1]:
import torch
import matplotlib.pyplot as plt
from torch.autograd import Variable
import os
import numpy as np
In [2]:
def equalizer(U):
""" equalizes the input matrix such that norms of the columns are equal
Args:
U: input matrix to be equalized.
Returns:
sol: the equalizer, such that U @ sol is equalized
"""
r = np.size(U,1)
G = U.T @ U
sol = np.identity(r)
for i in range(r):
l, u = np.linalg.eigh(G)
vi = np.mean(u,1)*np.sqrt(r-i)
vi = vi.reshape(len(vi),1)
if i<r-1:
tmp, _, _ = np.linalg.svd(vi)
else:
tmp = np.array([1])
Ri = np.identity(r)
Ri[i:,i:] = tmp
sol = sol @ Ri
G = tmp.T @ G @ tmp
if i<r-1:
G = G[1:,1:]
return sol;
In [3]:
def dropout_solution(M, r, lam):
""" solves the matrix dropout problem.
Args:
M: input matrix to be factorized
r: the factorization size
lam: the regularization parameter. lam = p/(1-p), where p is the dropout rate
Returns:
sol: sol = Ustar @ Vstar.T
Ustar, Vstar: a pair of factors which is globally optimum
"""
U, s, V = np.linalg.svd(M, full_matrices=False)
idx = s>0
s, U, V = s[idx], U[:,idx], V[idx,:] # compact SVD
rho_lim = min(len(s),r)
s, U, V = s[range(rho_lim)], U[:,range(rho_lim)], V[range(rho_lim),:].T
for i in reversed(range(1,rho_lim+1)):
kappa=np.mean(s[range(i)])
cutoff=(lam*i*kappa)/(r+lam*i)
if s[i-1] > cutoff:
s = s - cutoff
s[s<=0] = 0
s = np.append(s,np.zeros(r-rho_lim))
rootS = np.diag(np.sqrt(s))
U =np.hstack((U, np.zeros([np.size(U,0),r-rho_lim])))
V =np.hstack((V, np.zeros([np.size(V,0),r-rho_lim])))
eq=equalizer(np.matmul(U,rootS))
Ustar=U @ rootS @ eq
Vstar=V @ rootS @ eq
sol=Ustar @ Vstar.T
return (sol, Ustar, Vstar);
In [4]:
# N is sample size; D_in is input dimension;
# H is the number of hidden nodes; D_out is output dimension.
N, D_in, H, D_out = 1000, 20, 10, 15
torch.manual_seed(0)
# Create training points from standard Gaussian distribution
x = Variable(torch.randn(N, D_in))
# Create the transformation matrix M that generates labels y = M @ x.
um = np.random.randn(D_in, H)
vm = np.random.randn(H, D_out)
M = Variable(torch.Tensor(um @ vm), requires_grad=False)
# Regression: Create observed variables as y = M @ x
y = Variable(torch.mm(x,M).data, requires_grad=False)
# dropout rate and the corresponding regularizer
thetas = np.array([1/2,2/3,10/11])
lt = len(thetas)
# initialize performance measures
num_itr = 20000
obj = torch.Tensor(num_itr,lt).zero_()
loss_obj = torch.Tensor(num_itr,lt).zero_()
reg = torch.Tensor(num_itr,lt).zero_()
opt = torch.Tensor(lt).zero_()
# Lets initialize the weights to be highly non-equalized
scaling = torch.from_numpy(np.diag(2**np.arange(H))).float()
W1 = torch.mm(scaling,torch.rand(H, D_in))
W1 = torch.norm(M.data)**.5 * W1 / torch.norm(W1)
W2 = torch.mm(torch.rand(D_out, H), scaling)
W2 = torch.norm(M.data)**.5 * W2 / torch.norm(W2)
for theta_idx in range(lt):
theta = thetas[theta_idx]
lam = (1-theta)/theta
sol, U, V = dropout_solution(M.data, H, lam)
uu = np.diag(U.T@U)
vv = np.diag(V.T@V)
opt[theta_idx] = np.linalg.norm(M.data.numpy()-U@V.T)**2 + lam * np.sum(uu*vv)
model = torch.nn.Sequential(
torch.nn.Linear(D_in, H, bias=False),
torch.nn.Dropout(p=1-theta),
torch.nn.Linear(H, D_out, bias=False),
)
model[0].weight.data = W1.clone()
model[2].weight.data = W2.clone()
# MSE loss
loss_fn = torch.nn.MSELoss(size_average=False)
learning_rate = 1e-7
for t in range(num_itr):
y_pred = model(x)
loss = loss_fn(y_pred, y)
a=model[0].weight.data # H * D_in
b=model[2].weight.data # D_out * H
aa=torch.mm(a,torch.t(a)).diag()
bb=torch.mm(torch.t(b),b).diag()
obj[t,theta_idx] = torch.norm(M.data-torch.mm(torch.t(a),torch.t(b)))**2 + lam * torch.sum(aa*bb)
_, sm, _ = np.linalg.svd((b @ a).numpy(),full_matrices=False)
nuc_norm = sm.sum()
reg[t,theta_idx] = torch.std(aa*bb/(nuc_norm**2))
if t%2000 == 0:
print('itr: {}, truth: {}, loss: {}, obj: {}, var: {}'.format(t, opt[theta_idx], loss.data[0]*1/N, obj[t,theta_idx], reg[t,theta_idx]))
model.zero_grad()
loss.backward()
for param in model.parameters():
param.data -= learning_rate * param.grad.data
In [5]:
path = "../plots"
try:
if not os.path.exists(path):
os.makedirs(path)
except OSError:
print ("Creation of the directory %s failed" % path)
In [6]:
plt.close()
plt.figure()
for theta_idx in range(lt):
theta=thetas[theta_idx]
lam=(1-theta)/theta
plt.plot(range(num_itr),reg[:,theta_idx].numpy(),label=r"$\lambda$ = %.1f"%lam)
plt.legend()
plt.xlabel('Iterations')
plt.ylabel('Var of importance scores')
plt.title('Equalization')
plt.yscale('linear')
plt.xscale('log')
plt.grid()
plt.savefig("%s/nn_dout_r=%d.pdf"%(path,H))
plt.show()
In [7]:
plt.close()
for theta_idx in range(lt):
theta=thetas[theta_idx]
lam=(1-theta)/theta
plt.figure(theta_idx)
plt.plot(range(num_itr),obj[:,theta_idx].numpy(),label=r"Dropout, $\lambda$ = %.1f" %lam)
plt.plot(range(num_itr),opt[theta_idx]*np.ones(num_itr),label='Truth')
plt.legend()
plt.xlabel('Iterations')
plt.ylabel('Objective ')
plt.title('Convergence in objective')
plt.yscale('log')
plt.xscale('log')
plt.grid()
plt.savefig('%s/nn_dout_r=%d_theta=%.1f.pdf'%(path,H,thetas[theta_idx]))
plt.show()
In [ ]: