# Model fitting#

## Data model#

In this example we discuss model fitting and show an example with CVXPY. We are given data $$(x_i,y_i)\in \mathcal X \times \mathcal Y$$, $$i=1, \ldots, m$$:

• For $$\mathcal X= {\bf R}^n$$, $$x$$ is a feature vector.

• For $$\mathcal Y= {\bf R}$$, $$y$$ is a (real) outcome or label.

• For $$\mathcal Y= \{-1,1\}$$, $$y$$ is a (boolean) outcome.

Our goal is to find a model or predictor $$\psi: \mathcal X \to \mathcal Y$$ so that $$\psi(x)\approx y$$ for data $$(x,y)$$ that we haven’t seen:

• For $$\mathcal Y ={\bf R}$$, $$\psi$$ is a regression model.

• For $$\mathcal Y =\{-1,1\}$$, $$\psi$$ is a classifier.

We choose $$\psi$$ based on observed data and prior knowledge.

## Loss minimization model#

Let our data model be parametrized by $$\theta\in {\bf R}^n$$. We define a loss function $$L: \mathcal X \times \mathcal Y \times {\bf R}^n \to {\bf R}$$ where $$L(x_i,y_i,\theta)$$ is the loss (miss-fit) for the data point $$(x_i,y_i)$$, using the model parameter $$\theta$$.

We choose $$\theta$$ to minimize the total loss $$\sum_i L(x_i,y_i,\theta)$$. Our model is then $$\psi(x) = {\rm argmin}_y L(x,y,\theta)$$.

## Model fitting via regularized loss minimization#

An important concept in model fitting is regularization functions $$r:{\bf R}^n \to {\bf R} \cup \{\infty\}$$. The function $$r(\theta)$$ measures model complexity, enforces constraints, or represents a prior.

With regularization, we choose $$\theta$$ by minimizing the regularized loss

$(1/m) \sum_i L(x_i,y_i,\theta) + r(\theta).$

For many useful cases, this is a convex problem. Our model again is $$\psi(x) = {\rm argmin}_y L(x,y,\theta)$$.

## Example#

In the following code we do an example of model fitting with CVXPY. We are given (boolean) features $$z\in \{0,1\}^{10}$$ and (boolean) outcomes $$y\in \{-1,1\}$$. We generate a new feature vector $$x \in \{0,1\}^{55}$$ which contains all products $$z_iz_j$$ (co-occurence of pairs of original features).

To fit our model, we use logistic loss, or $$L(x,y,\theta) = \log (1+ \exp(-y\theta^T x))$$, and an $$\ell_1$$ regularizer $$r(\theta) = \|\theta\|_1$$. We train on $$m=200$$ examples and test on $$100$$ examples. We plot the train and test error as we vary $$\lambda$$.

# Generate data for logistic model fitting problem.
import numpy as np

# Construct Z given X.
def pairs(Z):
m, n = Z.shape
k = n * (n + 1) // 2
X = np.zeros((m, k))
count = 0
for i in range(n):
for j in range(i, n):
X[:, count] = Z[:, i] * Z[:, j]
count += 1
return X

np.random.seed(1)
n = 10
k = n * (n + 1) // 2
m = 200
TEST = 100
sigma = 1.9
DENSITY = 1.0
theta_true = np.random.randn(n, 1)
idxs = np.random.choice(range(n), int((1 - DENSITY) * n), replace=False)
for idx in idxs:
theta_true[idx] = 0

Z = np.random.binomial(1, 0.5, size=(m, n))
Y = np.sign(Z.dot(theta_true) + np.random.normal(0, sigma, size=(m, 1)))
X = pairs(Z)
X = np.hstack([X, np.ones((m, 1))])
Z_test = np.random.binomial(1, 0.5, size=(TEST, n))
Y_test = np.sign(Z_test.dot(theta_true) + np.random.normal(0, sigma, size=(TEST, 1)))
X_test = pairs(Z_test)
X_test = np.hstack([X_test, np.ones((TEST, 1))])

# Form model fitting problem with logistic loss and L1 regularization.
import cvxpy as cp

theta = cp.Variable((k + 1, 1))
lambd = cp.Parameter(nonneg=True)
loss = cp.sum(
cp.log_sum_exp(cp.hstack([np.zeros((m, 1)), -cp.multiply(Y, X @ theta)]), axis=1)
)
reg = cp.norm(theta[:k], 1)
prob = cp.Problem(cp.Minimize(loss / m + lambd * reg))

# Compute a trade-off curve and record train and test error.
TRIALS = 100
train_error = np.zeros(TRIALS)
test_error = np.zeros(TRIALS)
lambda_vals = np.logspace(-4, 0, TRIALS)
for i in range(TRIALS):
lambd.value = lambda_vals[i]
prob.solve(solver=cp.SCS)
train_error[i] = (
np.sign(Z.dot(theta_true)) != np.sign(X.dot(theta.value))
).sum() / m
test_error[i] = (
np.sign(Z_test.dot(theta_true)) != np.sign(X_test.dot(theta.value))
).sum() / TEST

# Plot the train and test error over the trade-off curve.
import matplotlib.pyplot as plt

%config InlineBackend.figure_format = 'svg'
plt.plot(lambda_vals, train_error, label="Train error")
plt.plot(lambda_vals, test_error, label="Test error")
plt.xscale("log")
plt.legend(loc="upper left")
plt.xlabel(r"$\lambda$", fontsize=16)
plt.show()


Below we plot $$|\theta_{k}|$$, $$k=1,\ldots,55$$, for the $$\lambda$$ that minimized the test error. Each $$|\theta_{k}|$$ is placed at position $$(i,j)$$ where $$z_iz_j = x_k$$. Notice that many $$\theta_{k}$$ are $$0$$, as we would expect with $$\ell_1$$ regularization.

# Solve model fitting problem with the lambda that minimizes test error.
idx = np.argmin(test_error)
lambd.value = lambda_vals[idx]
prob.solve(solver=cp.SCS)

# Plot the absolute value of the entries in theta corresponding to each feature.
P = np.zeros((n, n))
count = 0
for i in range(n):
for j in range(i, n):
P[i, j] = np.abs(theta.value[count])
count += 1
row_labels = range(1, n + 1)
column_labels = range(1, n + 1)

fig, ax = plt.subplots()
heatmap = ax.pcolor(P, cmap=plt.cm.Blues)

# put the major ticks at the middle of each cell
ax.set_xticks(np.arange(P.shape[1]) + 0.5, minor=False)
ax.set_yticks(np.arange(P.shape[0]) + 0.5, minor=False)

# want a more natural, table-like display
ax.invert_yaxis()
ax.xaxis.tick_top()

ax.set_xticklabels(column_labels, minor=False)
ax.set_yticklabels(row_labels, minor=False)

plt.xlabel(r"$z_i$", fontsize=16)
ax.xaxis.set_label_position("top")
plt.ylabel(r"$z_j$", fontsize=16)
plt.show()

/tmp/ipykernel_1888/4118099202.py:11: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
P[i, j] = np.abs(theta.value[count])