Lecture 04 CVX Demonstration

Ridge vs LASSO on Overfitting

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import eval_legendre
import cvxpy as cvx
np.set_printoptions(precision=2, suppress=True)

# Setup the problem
N = 20
x = np.linspace(-1,1,N)
a = np.array([1, 0.5, 0.5, 1.5, 1])
y = a[0]*eval_legendre(0,x) + a[1]*eval_legendre(1,x) + \
    a[2]*eval_legendre(2,x) + a[3]*eval_legendre(3,x) + \
    a[4]*eval_legendre(4,x) + 0.25*np.random.randn(N)

# Solve the regression problem
d = 20
theta_true = np.zeros(d)
theta_true[0:5] = a
X = np.zeros((N, d))
for p in range(d):
  X[:,p] = eval_legendre(p,x)

lambd = 0
theta       = cvx.Variable(d)
objective   = cvx.Minimize( cvx.sum_squares(X*theta-y) )
prob        = cvx.Problem(objective)
prob.solve()
theta_vanilla = theta.value

lambd = 0.5
theta       = cvx.Variable(d)
objective   = cvx.Minimize( cvx.sum_squares(X*theta-y) + lambd*cvx.sum_squares(theta) )
prob        = cvx.Problem(objective)
prob.solve()
theta_ridge = theta.value

lambd = 1
theta       = cvx.Variable(d)
objective   = cvx.Minimize( cvx.sum_squares(X*theta-y) + lambd*cvx.norm1(theta) )
prob        = cvx.Problem(objective)
prob.solve()
theta_LASSO = theta.value


# Plot the curves
t = np.linspace(-1, 1, 500)
Xhat = np.zeros((500,d))
for p in range(d):
  Xhat[:,p] = eval_legendre(p,t)

yhat_true = np.dot(Xhat, theta_true)
yhat_vanilla = np.dot(Xhat, theta_vanilla)
yhat_ridge = np.dot(Xhat, theta_ridge)
yhat_LASSO = np.dot(Xhat, theta_LASSO)
In [ ]:
# Regression Curves
fig = plt.figure(figsize = (12,8))
plt.figure(1) 
plt.subplot(221)
plt.plot(x, y, 'o')
plt.plot(t, yhat_true, linewidth=4, c='red')
plt.subplot(222)
plt.plot(x, y, 'o')
plt.plot(t, yhat_vanilla, linewidth=4, c='blue')
plt.subplot(223)
plt.plot(x, y, 'o')
plt.plot(t, yhat_ridge, linewidth=4, c='black')
plt.subplot(224)
plt.plot(x, y, 'o')
plt.plot(t, yhat_LASSO, linewidth=4, c='green')
Out[ ]:
[<matplotlib.lines.Line2D at 0x7fb11cc83390>]
In [ ]:
fig = plt.figure(figsize = (12,8))
plt.figure(1) 
plt.subplot(221)
plt.stem(theta_true,use_line_collection=True,linefmt='r-', markerfmt='ro',basefmt='r-')
plt.subplot(222)
plt.stem(theta_vanilla,use_line_collection=True,linefmt='b-', markerfmt='bo',basefmt='b-')
plt.subplot(223)
plt.stem(theta_ridge,use_line_collection=True,linefmt='k-', markerfmt='ko',basefmt='k-')
plt.subplot(224)
plt.stem(theta_LASSO,use_line_collection=True,linefmt='g-', markerfmt='go',basefmt='g-')
Out[ ]:
<StemContainer object of 3 artists>

Demo of CVX: Example 1

$$ \text{argmin}_{\mathbf{\theta}} \;\; \|\mathbf{X}\mathbf{\theta} - \mathbf{y}\|^2 \;\; \text{subject to} \;\; \|\mathbf{\theta}\|_1 \le \alpha $$$$ \text{argmin}_{\mathbf{\theta}} \;\; \|\mathbf{X}\mathbf{\theta} - \mathbf{y}\|^2 + \lambda \|\mathbf{\theta}\|_1 $$
In [ ]:
# LASSO linear regression
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import eval_legendre
import cvxpy as cvx
np.set_printoptions(precision=2, suppress=True)

N = 50
x = np.sort(2*np.random.rand(N)-1)
a = np.array([-0.001, 0.01, 0.55, 1.5, 1.2])
y = a[0]*eval_legendre(0,x) + a[1]*eval_legendre(1,x) + \
    a[2]*eval_legendre(2,x) + a[3]*eval_legendre(3,x) + \
    a[4]*eval_legendre(4,x) + 0.2*np.random.randn(N)
ysum = sum(np.abs(a))

P = 30
X = np.zeros((N, P))
for p in range(P):
  X[:,p] = eval_legendre(p,x)

# ===== Call CVX ==========
theta = cvx.Variable(P)
objective   = cvx.Minimize( cvx.sum_squares(X*theta-y) )
constraints = [cvx.norm1(theta) <= ysum]
prob        = cvx.Problem(objective, constraints)
prob.solve()
beta = theta.value
# =========================

t    = np.linspace(-1, 1, 500)
Xhat = np.zeros((500,P))
for p in range(P):
  Xhat[:,p] = eval_legendre(p,t)
yhat = np.dot(Xhat, beta)  
plt.plot(x,y,'o',markersize=12)
plt.plot(t,yhat, linewidth=4)
plt.axis((-1, 1, -3, 3));
plt.show()
print(beta)

[-0.03  0.    0.54  1.57  1.05 -0.    0.    0.02  0.    0.    0.    0.
 -0.01 -0.   -0.   -0.   -0.   -0.    0.   -0.    0.   -0.   -0.04 -0.
 -0.    0.    0.    0.    0.   -0.  ]

Demo of CVX: Example 2

$$ \text{argmin}_{\mathbf{\theta}} \log\left(\sum_{n=1}^N \exp(\mathbf{x}_n^T\mathbf{\theta} - y_n) \right) + \lambda \|\mathbf{\theta}\|^2 $$
In [3]:
import cvxpy as cvx
import numpy as np
import matplotlib.pyplot as plt

N = 100
d = 2
X = np.random.randn(N, d)
y = np.random.randn(N)
lambd = 0.01


# Setup the grid for visualization
xgrid = np.linspace(-1, 1, 200)
ygrid = np.linspace(-1, 1, 200)
theta0, theta1 = np.meshgrid(xgrid, ygrid)
cost = np.zeros((len(xgrid),len(ygrid)))
for i in range(len(xgrid)):
  for j in range(len(ygrid)):
    theta_tmp = np.hstack((theta0[i,j],theta1[i,j]))
    cost[i,j] = np.log(np.sum(np.exp(np.dot(X, theta_tmp)-y))) + \
              lambd*np.sum(theta_tmp**2)

# ===== Call CVX ==========
theta = cvx.Variable(d)
objective = cvx.Minimize(cvx.log_sum_exp(X*theta-y) + lambd*cvx.sum_squares(theta))
constraints = []
prob = cvx.Problem(objective, constraints)
optimal_objective_value = prob.solve()
thetahat = theta.value
# =========================
print("status:", prob.status)
print(optimal_objective_value)
print(theta.value)


# Plotting
fig = plt.figure(figsize = (6,6))
plt.contour(xgrid, ygrid, cost, 20, cmap='RdGy')
plt.scatter(thetahat[0], thetahat[1],c='blue')  

status: optimal
5.306236216397515
[-0.01 -0.34]
Out[3]:
<matplotlib.collections.PathCollection at 0x7f4776b12898>

Demo of CVX: Example 3

$$ \text{argmin}_{x} \; \| x - x_0\|^2 \;\; \text{subject to} \;\; Ax = y $$
In [4]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import eval_legendre
import cvxpy as cvx
np.set_printoptions(precision=2, suppress=True)

m = 3
n = 2*m
A = np.random.randn(m,n)
xstar = np.random.randn(n)
y = np.dot(A, xstar)
x0 = np.random.randn(n)

# ===== Call CVX ==========
x = cvx.Variable(n)
objective   = cvx.Minimize( cvx.sum_squares(x-x0) )
constraints = [A*x == y]
prob        = cvx.Problem(objective, constraints)
prob.solve()
theta = x.value
# =========================

# ===== Direct inversion ======
AAt = np.dot(A, A.T)
e   = y-np.dot(A,x0)
theta1 = x0 + np.dot( A.T, np.dot( np.linalg.pinv(AAt), e ) )
# =============================

print(theta, theta1)

[ 0.38 -1.22  0.13 -0.23  1.36 -0.67] [ 0.38 -1.22  0.13 -0.23  1.36 -0.67]
In [ ]:
%%shell
jupyter nbconvert --to html /content/ECE595_lecture04_cvx.ipynb

[NbConvertApp] Converting notebook /content/ECE595_lecture04_cvx.ipynb to html
[NbConvertApp] Writing 483362 bytes to /content/ECE595_lecture04_cvx.html
Out[ ]: