import numpy as np

def soft_thresholding(rho, lambda_):
    if rho < - lambda_:
        return (rho + lambda_)
    elif rho > lambda_:
        return (rho - lambda_)
    else:
        return 0

def lasso_coordinate_descent(X, y, lambda_, 
                             num_iters=100, tol=1e-4,
                            verbose = False):
    m, n = X.shape
    beta = np.zeros(n)
    beta_prev = np.zeros(n)
    
    for iteration in range(num_iters):
        for j in range(n):
            X_j = X[:, j]
            residual = y - X @ beta + beta[j] * X_j  # partial residual
            rho = np.dot(X_j, residual)
            beta[j] = soft_thresholding(rho, lambda_)  # update rule
        
        # Check for convergence
        if np.linalg.norm(beta - beta_prev, ord=2) < tol:
            if verbose: 
                print(f"Converged in {iteration + 1} iterations.")
            break
        
        beta_prev = beta.copy()
    
    return beta

import numpy as np
import matplotlib.pyplot as plt

# random seed
np.random.seed(24)

# Size of signal
n = 256

# Sparsity (# nonzeros) in the signal
s = 10

# Number of samples (undersample by a factor of 8)
m = 64

# Generate and display the signal
x0 = np.zeros(n)
nonzero_indices = np.random.choice(np.arange(n), s)
x0[nonzero_indices] = np.random.randn(s)

# Generate the random sampling matrix
A = np.random.randn(m, n) / m

# Subsample by multiplexing
y = A.dot(x0)

# Plot the true signal
plt.figure()
plt.title("True Signal x0")
plt.xlabel("Index")
plt.ylabel("X0")
plt.plot(np.arange(1, n+1), x0)
plt.show()

# Example usage
lambda_ = 0.01

beta = lasso_coordinate_descent(A, y, lambda_,
                               verbose = True)

plt.plot(np.arange(1, n+1), beta)

Converged in 3 iterations.

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

# Range of lambda values
lambda_values = np.logspace(-1,-2, 100)

# Store the coefficients for each lambda
coefficients = []
errors = []

for lambda_ in lambda_values:
    beta = lasso_coordinate_descent(A, y, lambda_)
    coefficients.append(beta)
    predictions = A @ beta
    error = np.mean((y - predictions) ** 2)  # Mean Squared Error
    errors.append(error)

coefficients = np.array(coefficients)
errors = np.array(errors)

def lasso_greedy_coordinate_descent(X, y, lambda_, num_iters=100, tol=1e-4,
                            verbose = False):
    m, n = X.shape
    beta = np.zeros(n)
    beta_prev = np.zeros(n)
    
    for iteration in range(num_iters):
        max_decrease = 0
        for j in range(n):
            X_j = X[:, j]
            residual = y - X @ beta + beta[j] * X_j  # partial residual
            rho = np.dot(X_j, residual)
            new_beta_j = soft_thresholding(rho, lambda_)
            
            # Calculate the decrease in the objective function
            decrease = abs(beta[j] - new_beta_j)
            if decrease > max_decrease:
                max_decrease = decrease
                max_index = j
                best_beta_j = new_beta_j
        
        # Update the coordinate with the largest decrease
        beta[max_index] = best_beta_j
        
        # Check for convergence
        if np.linalg.norm(beta - beta_prev, ord=2) < tol:
            if verbose:
                print(f"Converged in {iteration + 1} iterations.")
            break
        
        beta_prev = beta.copy()
    
    return beta

# Example usage
lambda_ = 0.01

beta = lasso_greedy_coordinate_descent(A, y, lambda_,
                               verbose = True)

plt.plot(np.arange(1, n+1), beta)

Converged in 23 iterations.

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

# Plotting the coefficient paths and prediction errors side by side
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Coefficient paths plot
for i in range(coefficients.shape[1]):
    ax1.plot(lambda_values, coefficients[:, i])
ax1.set_xscale('log')
ax1.invert_xaxis()  # Invert the x-axis
ax1.set_xlabel('Lambda')
ax1.set_ylabel('Coefficients')
ax1.set_title('Lasso Paths')

# Prediction errors plot
ax2.plot(lambda_values, errors)
ax2.set_xscale('log')
ax2.invert_xaxis()  # Invert the x-axis
ax2.set_xlabel('Lambda')
ax2.set_ylabel('Mean Squared Error')
ax2.set_title('Prediction Error for Different Lambda Values')

plt.tight_layout()
plt.show()

from scipy.optimize import linprog

## Method 1

vF = np.ones(2 * n)

mAeq = np.hstack((A, -A))
vBeq = y

vLowerBound = np.zeros(2 * n)
vUpperBound = np.inf * np.ones(2 * n)

res = linprog(vF, A_eq=mAeq, b_eq=vBeq, bounds=list(zip(vLowerBound, vUpperBound)))

vX = res.x[:n] - res.x[n:]

np.allclose(x0, vX)

True

import cvxpy as cp
import numpy as np

x = cp.Variable(n)

# Create the optimization problem
objective = cp.Minimize(cp.norm(x, 1))
constraints = [A @ x == y]
problem = cp.Problem(objective, constraints)

# Solve the problem
problem.solve()

# Retrieve the solution
x_sol = x.value

print("obj val=", problem.solve())

np.allclose(x0, x_sol)

obj val= 11.788696421666058

True

from IPython.display import Image, display

# Adjust the file path as necessary
file_path = "qp.jpg"

# Display the image with 80% width
display(Image(filename=file_path, width=600))

# Adjust the file path as necessary
file_path = "ridge.jpg"

# Display the image with 80% width
display(Image(filename=file_path, width=600))

# Adjust the file path as necessary
file_path = "lasso.jpg"

# Display the image with 80% width
display(Image(filename=file_path, width=600))

# Adjust the file path as necessary
file_path = "lena.jpg"

# Display the image with 80% width
display(Image(filename=file_path, width=600))

STAT 207: Advanced Optimization Topics¶

Linear Programming¶

Examples¶

Lasso Problem¶

Quantile regression¶

$\ell_1$ Regression¶

Dantzig selector¶

Quadratic Programming¶

Examples¶