Max's pricing problem. #136
Description
This issue is to help the great Max Schaller with his pricing problem.
@maxschaller I used claude to try and run your full pricing problem. It changed a few things: it first changes the variable to be C_T (C transpose directly) and then does some fancy linear algebra trick to initialize the problem well. Please let me know if these changes make any sense.
import cvxpy as cp
import numpy as np
import time
# simulate data
n, N = 100, 190
np.random.seed(0)
Psim = 1 + np.random.rand(n, N)
p_nom_sim = np.mean(Psim, axis=1, keepdims=True)
Pi = np.log(Psim / p_nom_sim)
Pitilde = np.vstack((Pi, np.ones((1, Pi.shape[1]))))
diagsim = np.diag(np.random.uniform(-5.0, -1.0, size=n))
tall = np.clip(np.random.randn(n, 10), -1, 1)
lowranksim = tall @ tall.T
Esim = diagsim + 0.3 * lowranksim
D = np.exp(Esim @ Pi) # assumes that all dnom are 1
rank = 10
lam = 0.1
# variables
Etilde = cp.Variable((n, n + 1))
B = cp.Variable((n, rank))
# Reformulation: define C_T directly instead of C.T to avoid matrix transpose
C_T = cp.Variable((rank, n)) # This is effectively C.T
s = cp.Variable(n)
# Initialize variables to help convergence
# Use the simulated values as a starting point
Etilde_init = np.hstack([Esim, np.zeros((n, 1))])
Etilde.value = Etilde_init
# SVD of Esim to initialize B and C_T
U, S, Vt = np.linalg.svd(Esim)
B.value = U[:, :rank] @ np.diag(np.sqrt(S[:rank]))
C_T.value = np.diag(np.sqrt(S[:rank])) @ Vt[:rank, :]
# Initialize s from diagonal of Esim minus the low-rank approximation
s.value = np.diag(Esim) - np.diag(B.value @ C_T.value)
# objective and constraint
# Note: For the regularization term, we use C_T instead of C
# sum_squares(C) = sum_squares(C_T) since Frobenius norm is transpose-invariant
f = cp.sum(cp.multiply(D, Etilde @ Pitilde) - cp.exp(Etilde @ Pitilde)) / N
obj = cp.Maximize(f - (lam / 2) * (cp.sum_squares(B) + cp.sum_squares(C_T)))
constr = [Etilde[:, :-1] == B @ C_T + cp.diag(s)]
# solve problem
prob = cp.Problem(obj, constr)
start_time = time.time()
prob.solve(nlp=True, verbose=True)
end_time = time.time()
print(f"\n{'='*50}")
print(f"Total solve time: {end_time - start_time:.2f} seconds")
print(f"Optimal value: {prob.value}")
print(f"Status: {prob.status}")The problem runs successfully in 6 iterations after this new initialization.
I will post the full logs in the comments.