How to add L1 norm as a constraint in PCA

回答済み

Alvaro Mendez Civieta

2020年12月11日 11:12

• 
• 
• 

I am trying to solve the PCA problem adding an extra \(L_1\) constraint into it. Given a data matrix \(X\), my objective function to solve is 

\[max_w{w^tX^tXw} \quad s.t. \quad \|w\|_2=1, \quad \|w\|_1 \leq t\]

I have already worked out how to solve the objective function + the first constraint but I have no idea on how to include the \(L_1\) constraint.

 

So far this is what I have,

import gurobipy as gp
from gurobipy import GRB
import numpy as np
from sklearn.datasets import load_boston
from itertools import product

# Load data
boston = load_boston()
x = boston.data
x = x[:, [0, 2, 4, 5, 6, 7, 10, 11, 12]]  # select non-categorical variables
x_center = x - np.mean(x, 0)
n, p = x_center.shape
Quad = np.dot(x_center.T, x_center)

# Define model
pca_gurobi = gp.Model()
weights = pca_gurobi.addVars(p, lb=-GRB.INFINITY, name="weights")  # pca loading
obj_pca = (1/n)*sum(Quad[i, j] * weights[i] * weights[j] for i, j in product(range(p), repeat=2))
pca_gurobi.addConstr(sum(weights[i]**2 for i in range(p)) == 1)

# Add L1 constraint
# lasso_threshold = 5.0
# First try (Error due to sum over GenExpr)
# pca_gurobi.addConstr(sum(gp.abs_(weights[i]) for i in range(p)) <= lasso_threshold)

# Second try: Error
# abs_weights = pca_gurobi.addVars(p, name="abs_weights")
# for v, absv in zip(weights, abs_weights):
#     pca_gurobi.addConstr(gp.abs_(v)==absv)
# pca_gurobi.addConstr(sum(abs_w for abs_w in abs_weights.tolist()) <= lasso_threshold)

# Solve model
pca_gurobi.setObjective(obj_pca, GRB.MAXIMIZE)
pca_gurobi.params.NonConvex = 2  # Let gurobi know that constraint is non convex
pca_gurobi.optimize()
weights_pca_gurobi = np.array([weights[i].X for i in range(p)])

My second approach on adding the L1 norm was based on this post here:

https://support.gurobi.com/hc/en-us/community/posts/360073709712-Norm-1-of-a-MVar-as-Constraint?input_string=How%20to%20add%20L1%20norm%20as%20a%20constraint 

And I saw that the answer was based on using MVar instead of normal variables, so I converted my problem to use MVar, but it somehow resulted into an infeasible problem. This is the code for my second approach:

spca_gurobi2 = gp.Model()
s_weights2 = spca_gurobi2.addMVar(p, lb=-GRB.INFINITY, name="sparse_weights")  # pca loading

lasso_threshold = 5.0

obj_spca2 = 0
for i, j in product(range(p), repeat=2):
    list_w = s_weights2.tolist()
    obj_spca2 += Quad[i, j] * list_w[i] * list_w[j]
obj_spca2 = (1/n)*obj_spca2

spca_gurobi2.addConstr(sum(w**2 for w in s_weights2.tolist()) == 1)
# First try
abs_s_weights2 = spca_gurobi2.addMVar(p, name="abs_sparse_weights")
for v, absv in zip(s_weights2.tolist(), abs_s_weights2.tolist()):
    spca_gurobi2.addConstr(absv == gp.abs_(v))
spca_gurobi2.addConstr(sum(abs_w for abs_w in abs_s_weights2.tolist()) <= lasso_threshold)

spca_gurobi2.setObjective(obj_spca2, GRB.MAXIMIZE)
spca_gurobi2.params.NonConvex = 2  # Let gurobi know that constraint is non convex

spca_gurobi2.optimize()
weights_spca_gurobi2 = np.array([s_weights2[i].X for i in range(p)])

I am pretty new to Gurobi so I would really appreciate any advice here.

0

• 正式なコメント

  

  Simranjit Kaur

  • Gurobi Staff

  2025年06月23日 03:03 編集済み

  This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?.
• 

  Matthias Miltenberger

  • Gurobi Staff

  2020年12月15日 12:36

  Hi Alvaro,

  You need to write the L1 constraints like this:

  for v, absv in zip(weights, abs_weights):
      pca_gurobi.addConstr(abs_weights[absv] == gp.abs_(weights[v]))

  Please note that looping over the variable dicts will generate index objects as looping variables, so you need to use those to specify the actual decision variables in the respective dicts weights and abs_weights.

  The error message here is quite misleading because the input data for the addConstr call doesn't contain a variable.

  I hope that helps.

  Cheers,
  Matthias

  0

投稿コメントは受け付けていません。