Exercises

import numpy as np
from scipy.optimize import minimize

def f(x):
    return (x[0] - 3)**2 + (x[1] + 1)**2

def grad_f(x):
    return np.array([2*(x[0] - 3), 2*(x[1] + 1)])

res = minimize(f, [0, 0], method='BFGS', jac=grad_f)
print(f"x* = {res.x}")  # [3.0, -1.0]
print(f"f(x*) = {res.fun:.2e}")  # ~0
print(f"Converged: {res.success}")

import numpy as np
from scipy.optimize import brentq

f = lambda x: np.exp(x) - 3*x
x_star = brentq(f, 0, 2)
print(f"Root: {x_star:.10f}")
print(f"f(x*) = {f(x_star):.2e}")

import numpy as np

def soft_threshold(v, lam):
    return np.sign(v) * np.maximum(np.abs(v) - lam, 0)

v = np.array([-3, -1, 0.5, 0, 2, -0.3, 4])
result = soft_threshold(v, 1.0)
print(f"Result: {result}")
# [-2, 0, 0, 0, 1, 0, 3]
print(f"Zeros: {np.sum(result == 0)}")  # 4

from scipy.optimize import linprog

c = [-3, -5]  # negate for maximization
A_ub = [[1, 1], [2, 1]]
b_ub = [10, 14]
bounds = [(0, None), (0, None)]

res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds)
print(f"x={res.x[0]:.1f}, y={res.x[1]:.1f}")
print(f"Max profit = {-res.fun:.1f}")

import numpy as np
from scipy.optimize import fixed_point

x_star = fixed_point(np.cos, 0.5)
print(f"Fixed point: {x_star:.10f}")
print(f"cos(x*) = {np.cos(x_star):.10f}")
print(f"|x* - cos(x*)| = {abs(x_star - np.cos(x_star)):.2e}")

import numpy as np
from scipy.optimize import minimize

def rosen(x):
    return (1-x[0])**2 + 100*(x[1]-x[0]**2)**2

def rosen_grad(x):
    return np.array([
        -2*(1-x[0]) - 400*x[0]*(x[1]-x[0]**2),
        200*(x[1]-x[0]**2)
    ])

def rosen_hess(x):
    return np.array([
        [2 - 400*(x[1]-3*x[0]**2), -400*x[0]],
        [-400*x[0], 200]
    ])

x0 = np.array([-1.5, 2.0])
for method in ['Nelder-Mead', 'CG', 'BFGS', 'Newton-CG']:
    kw = {}
    if method in ('CG', 'BFGS', 'Newton-CG'):
        kw['jac'] = rosen_grad
    if method == 'Newton-CG':
        kw['hess'] = rosen_hess
    res = minimize(rosen, x0, method=method, **kw)
    print(f"{method:12s}: nit={res.nit:4d}, nfev={res.nfev:4d}, "
          f"f*={res.fun:.2e}")

import numpy as np

N = 12
gains = 10**((20 - 3*np.arange(N))/10)
P, sigma2 = 20.0, 1.0
noise_floor = sigma2 / gains
idx = np.argsort(noise_floor)
nf_sorted = noise_floor[idx]

for K in range(N, 0, -1):
    mu = (P + np.sum(nf_sorted[:K])) / K
    p = np.maximum(mu - nf_sorted, 0)
    if p[K-1] > 0:
        break

powers = np.zeros(N)
powers[idx] = p
print(f"Zero-power channels: {np.sum(powers == 0)}")

import numpy as np
import cvxpy as cp

np.random.seed(42)
m, n, k = 50, 200, 10
A = np.random.randn(m, n)
x_true = np.zeros(n); x_true[:k] = np.random.randn(k)
b = A @ x_true + 0.1*np.random.randn(m)
lam = 0.5

# CVXPY
x_cv = cp.Variable(n)
prob = cp.Problem(cp.Minimize(0.5*cp.sum_squares(A@x_cv-b) + lam*cp.norm(x_cv,1)))
prob.solve()

# ISTA
L = np.linalg.norm(A.T @ A, 2)
x_is = np.zeros(n)
for _ in range(1000):
    x_is = np.sign(x_is - A.T@(A@x_is-b)/L) * np.maximum(
        np.abs(x_is - A.T@(A@x_is-b)/L) - lam/L, 0)

print(f"||x_cvxpy - x_ista|| = {np.linalg.norm(x_cv.value - x_is):.4f}")

from scipy.optimize import fsolve
import numpy as np

def sys(v):
    x, y = v
    return [x**2 + y**2 - 25, x**2 - y - 3]

def jac(v):
    x, y = v
    return [[2*x, 2*y], [2*x, -1]]

starts = [(3,3), (-3,3), (3,-3), (-3,-3)]
solutions = set()
for s in starts:
    sol, info, ier, msg = fsolve(sys, s, fprime=jac, full_output=True)
    if ier == 1:
        solutions.add(tuple(np.round(sol, 8)))

for s in sorted(solutions):
    print(f"({s[0]:.4f}, {s[1]:.4f})")

import numpy as np

def proj_simplex(v):
    n = len(v)
    u = np.sort(v)[::-1]
    cssv = np.cumsum(u) - 1
    rho = np.max(np.where(u > cssv / np.arange(1, n+1)))
    theta = cssv[rho] / (rho + 1)
    return np.maximum(v - theta, 0)

v = np.array([1.5, -0.5, 0.3, 0.8, -0.1])
w = proj_simplex(v)
print(f"Projection: {w}")
print(f"Sum: {np.sum(w):.10f}")
print(f"All >= 0: {np.all(w >= -1e-10)}")

import numpy as np
from scipy.optimize import minimize, LinearConstraint, Bounds

f = lambda x: x[0]**2 + x[1]**2
grad_f = lambda x: 2*x
hess_f = lambda x: 2*np.eye(2)

lc = LinearConstraint([[1, 1]], [1], [1])
bounds = Bounds([0.1, 0.1], [np.inf, np.inf])

res = minimize(f, [0.5, 0.5], method='trust-constr',
               jac=grad_f, hess=hess_f,
               constraints=lc, bounds=bounds)
print(f"x* = {res.x}")
print(f"f* = {res.fun:.6f}")

ISTA converges as $O(1/k)$ while FISTA converges as $O(1/k^2)$. On the log-scale plot, ISTA should show slope $-1$ and FISTA slope $-2$ (asymptotically).

The Goemans-Williamson algorithm guarantees at least $0.878$ of the optimal cut in expectation. With 1000 random hyperplanes, you should find a good cut.

Anderson acceleration typically reduces iteration count by 3-10x compared to plain iteration. The residual norms should decrease much faster.

Small $\lambda$ preserves detail but retains noise; large $\lambda$ over-smooths. The optimal $\lambda$ trades off between the two. TV denoising preserves sharp edges, unlike Gaussian smoothing.

For large $n$, Hessian-vector products cost $O(mn)$ while forming and applying the Hessian costs $O(mn + n^2)$. The hessp approach should be faster when $n$ is large.

ADMM decomposes the problem into easy subproblems. The $x$-update is a least-squares solve; the $z$-update is soft-thresholding. Convergence may require 100-500 iterations but each is cheap.

After fixing the phase of the desired signal constraint, this becomes an SOCP. CVXPY handles complex variables natively. The resulting beamformer should show a beam toward $\theta_0$ and nulls toward interferers.