Robust PCA¶
Solve the robust PCA problem taken from Xinyang Yi, et al. “Fast algorithms for robust PCA via gradient descent.” Advances in neural information processing systems. 2016.
Problem Description¶
\[\min_{M,S}||M||_{\text{nuc}}+\lambda||S||_1,\]
\[\text{s.t. }Y=M+S,\]
where \(M,S\in R^{d_1,d_2}\) are matrix form optimization variables, \(Y\in R^{d_1,d_2}\) is a given matrix, and \(||\cdot||_{\text{nuc}}\) denotes the nuclear norm.
Modules Importing¶
Import all necessary modules and add PyGRANSO src folder to system path.
[1]:
import time
import torch
from pygranso.pygranso import pygranso
from pygranso.pygransoStruct import pygransoStruct
Data Initialization¶
Specify torch device, and generate data.
Use GPU for this problem. If no cuda device available, please set device = torch.device(‘cpu’)
[2]:
device = torch.device('cuda')
d1 = 3
d2 = 4
torch.manual_seed(1)
eta = .05
# All the user-provided data (vector/matrix/tensor) must be in torch tensor format.
# As PyTorch tensor is single precision by default, one must explicitly set `dtype=torch.double`.
# Also, please make sure the device of provided torch tensor is the same as opts.torch_device.
Y = torch.randn(d1,d2).to(device=device, dtype=torch.double)
Function Set-Up¶
Encode the optimization variables, and objective and constraint functions.
Note: please strictly follow the format of comb_fn, which will be used in the PyGRANSO main algortihm.
[3]:
# variables and corresponding dimensions.
var_in = {"M": [d1,d2],"S": [d1,d2]}
def user_fn(X_struct,Y):
M = X_struct.M
S = X_struct.S
# objective function
f = torch.norm(M, p = 'nuc') + eta * torch.norm(S, p = 1)
# inequality constraint, matrix form
ci = None
# equality constraint
ce = pygransoStruct()
ce.c1 = M + S - Y
return [f,ci,ce]
comb_fn = lambda X_struct : user_fn(X_struct,Y)
User Options¶
Specify user-defined options for PyGRANSO
[4]:
opts = pygransoStruct()
opts.torch_device = device
opts.print_frequency = 10
opts.x0 = .2 * torch.ones((2*d1*d2,1)).to(device=device, dtype=torch.double)
opts.opt_tol = 1e-6
Main Algorithm¶
[5]:
start = time.time()
soln = pygranso(var_spec = var_in,combined_fn = comb_fn,user_opts = opts)
end = time.time()
print("Total Wall Time: {}s".format(end - start))
╔═════ QP SOLVER NOTICE ════════════════════════════════════════════════════════════════════════╗
║ PyGRANSO requires a quadratic program (QP) solver that has a quadprog-compatible interface, ║
║ the default is osqp. Users may provide their own wrapper for the QP solver. ║
║ To disable this notice, set opts.quadprog_info_msg = False ║
╚═══════════════════════════════════════════════════════════════════════════════════════════════╝
═════════════════════════════════════════════════════════════════════════════════════════════════════════════════╗
PyGRANSO: A PyTorch-enabled port of GRANSO with auto-differentiation ║
Version 1.2.0 ║
Licensed under the AGPLv3, Copyright (C) 2021-2022 Tim Mitchell and Buyun Liang ║
═════════════════════════════════════════════════════════════════════════════════════════════════════════════════╣
Problem specifications: ║
# of variables : 24 ║
# of inequality constraints : 0 ║
# of equality constraints : 12 ║
═════╦═══════════════════════════╦════════════════╦═════════════════╦═══════════════════════╦════════════════════╣
║ <--- Penalty Function --> ║ ║ Total Violation ║ <--- Line Search ---> ║ <- Stationarity -> ║
Iter ║ Mu │ Value ║ Objective ║ Ineq │ Eq ║ SD │ Evals │ t ║ Grads │ Value ║
═════╬═══════════════════════════╬════════════════╬═════════════════╬═══════════════════════╬════════════════════╣
0 ║ 1.000000 │ 9.26721133015 ║ 0.81282032303 ║ - │ 1.922768 ║ - │ 1 │ 0.000000 ║ 1 │ 2.502220 ║
10 ║ 0.071790 │ 0.03336327939 ║ 0.46473565836 ║ - │ 1.67e-15 ║ S │ 1 │ 1.000000 ║ 1 │ 11.88437 ║
20 ║ 0.011973 │ 0.00404030852 ║ 0.33746530772 ║ - │ 1.33e-15 ║ S │ 1 │ 1.000000 ║ 1 │ 0.002421 ║
30 ║ 0.004175 │ 0.00139375036 ║ 0.33386777336 ║ - │ 7.99e-15 ║ S │ 2 │ 0.500000 ║ 1 │ 112.4975 ║
40 ║ 0.001997 │ 6.6642624e-04 ║ 0.33376748890 ║ - │ 2.20e-14 ║ S │ 4 │ 0.125000 ║ 4 │ 3.07e-05 ║
═════╩═══════════════════════════╩════════════════╩═════════════════╩═══════════════════════╩════════════════════╣
F = final iterate, B = Best (to tolerance), MF = Most Feasible ║
Optimization results: ║
═════╦═══════════════════════════╦════════════════╦═════════════════╦═══════════════════════╦════════════════════╣
F ║ │ ║ 0.33375284367 ║ - │ 6.48e-14 ║ │ │ ║ │ ║
B ║ │ ║ 0.33375284367 ║ - │ 6.48e-14 ║ │ │ ║ │ ║
MF ║ │ ║ 1.30399559824 ║ - │ 1.11e-16 ║ │ │ ║ │ ║
═════╩═══════════════════════════╩════════════════╩═════════════════╩═══════════════════════╩════════════════════╣
Iterations: 46 ║
Function evaluations: 104 ║
PyGRANSO termination code: 0 --- converged to stationarity and feasibility tolerances. ║
═════════════════════════════════════════════════════════════════════════════════════════════════════════════════╝
Total Wall Time: 2.862501621246338s