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)+ 0.5 * np.random.randn(m)


# 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_ = .1

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 0x3004886d0>]

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.1

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

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

Converged in 44 iterations.

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

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

# Subsample by multiplexing
y_test = A_test.dot(x0)+ 0.5 * np.random.randn(m)

# Range of lambda values
lambda_values = np.logspace(0,-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_test @ beta
    error = np.mean((y_test - predictions) ** 2)  # Mean Squared Error
    errors.append(error)

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

# 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)

False

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= 164.32041829954989

False

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))

import numpy as np
import matplotlib.pyplot as plt
from skimage import io

# Load the combined image
image = io.imread('lena.jpg', as_gray=True)

# Get image dimensions
height, width = image.shape

# Split into left and right halves
midpoint = width // 2
lena_noisy = image[:, :midpoint]
lena_clean = image[:, midpoint:]

# Display both images
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.title("Lena Noisy")
plt.imshow(lena_noisy, cmap='gray')
plt.axis('off')

plt.subplot(1, 2, 2)
plt.title("Lena Clean")
plt.imshow(lena_clean, cmap='gray')
plt.axis('off')

plt.tight_layout()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from skimage import io, color
from skimage.restoration import denoise_tv_chambolle

# Apply total variation denoising (isotropic TV)
lambda_tv = 0.1  # regularization strength
img_denoised = denoise_tv_chambolle(lena_noisy, weight=lambda_tv)

# Plot
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.title("Original (Noisy)")
plt.imshow(lena_noisy, cmap='gray')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.title("Denoised (TV)")
plt.imshow(img_denoised, cmap='gray')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.title("Ground Truth")
plt.imshow(lena_clean, cmap='gray')
plt.axis('off')

plt.tight_layout()
plt.show()

import cvxpy as cp

Y = lena_noisy  # input image
n, m = Y.shape
print(n,m)
X = cp.Variable((n, m))

tv_penalty = cp.tv(X)
# Denoising objective
lambda_tv = 0.1
objective = cp.Minimize(0.5 * cp.sum_squares(X - Y) + lambda_tv * tv_penalty)

# Solve
problem = cp.Problem(objective)
problem.solve()

# Result
img_denoised = X.value

351 343

# Plot
plt.figure(figsize=(12, 4))
plt.subplot(1, 3, 1)
plt.title("Original (Noisy)")
plt.imshow(lena_noisy, cmap='gray')
plt.axis('off')

plt.subplot(1, 3, 2)
plt.title("Denoised (TV)")
plt.imshow(img_denoised, cmap='gray')
plt.axis('off')

plt.subplot(1, 3, 3)
plt.title("Ground Truth")
plt.imshow(lena_clean, cmap='gray')
plt.axis('off')

plt.tight_layout()
plt.show()

STAT 207: Advanced Optimization Topics¶

Linear Programming¶

Examples¶

Lasso Problem¶

Quantile regression¶

$\ell_1$ Regression¶

Dantzig selector¶

Quadratic Programming¶

Examples¶