How to solve a semi-definite programming with l1-norm objective?

I am currently working on a problem related to a semi-definite programming (SDP) with l1-norm objective. After two-week survey in the literature, I still can't find a possible way to solve the following problem:

min_A   -tr(A*C) + \lambda*||A-A0||
 s.t.   A>=0, tr(A)<=\gamma

where C and A0 are given nxn constant matrices, and \lambda is a coefficient term that balances the contribution. The goal is to find an nxn matrix A such that the objective function is minimized such that A is POSITIVE SEMI-DEFINITE and its trace is less than some threshold \gamma.

I wonder if there is any paper or solver that has handled this kind of problem. Thanks in advance if anyone could provide any clues.

Answer 1

[John D'Errico]

This is an optimization problem, but I would argue it has been posed poorly. If you want to constrain the matrix A to be positive semi-definite, then DON'T solve for the matrix A. Instead, solve for a Cholesky factor of A. So your unknowns will be a triangular matrix.

For example, suppose we want to solve a problem for the 3x3 matrix A, such that A is P.S.D.? I'll pick a semi-random problem. Then I'll solve it using a problem based formulation.

help optimproblem

 OPTIMPROBLEM Create an optimization problem
 
    PROB = OPTIMPROBLEM creates an optimization problem with default
    properties.
 
    PROB = OPTIMPROBLEM(NAME, VALUE, ...) creates an optimization problem 
    with specified values of optional parameters:
 
      'ObjectiveSense'   Sense of optimization: 'maximize', 
                         'minimize' (default), 'max', 'min', or a struct
                         with each field containing one of those strings
      'Objective'        An OptimizationExpression array or a struct with 
                         scalar OptimizationExpressions as fields
      'Constraints'      An OptimizationConstraint array or a struct with 
                         OptimizationConstraint arrays as fields
      'Description'      Problem description: character array or string.
    
    The following simple example illustrates how an optimization problem is
    created and solved
 
    % Create an optimization problem
    prob = optimproblem('Description', 'Simple Example');
   
    % Create a 2-by-1 optimization variable
    x = optimvar('x', 2, 1, 'LowerBound', 0);
 
    % Define the objective
    prob.Objective = -5*x(1) - x(2);
 
    % Define the constraints
    prob.Constraints.MyConstraint1 = x(1) + x(2) <= 5;
    prob.Constraints.MyConstraint2 = 2*x(1) + 0.5*x(2) <= 8;
 
    % Solve the problem
    sol = solve(prob);
 
    See also OPTIMVAR, OPTIMCONSTR, OPTIMEXPR

    Documentation for optimproblem
       doc optimproblem

% arbitrary numbers for lambda and the maximum trace of A.
trmax = 3;
lambda = 1; 
C = rand(3);
A0 = rand(3); A0 = A0 + A0' + randn(size(A0))/100

A0 = 3×3

1.8589 1.4917 1.2471 1.4817 1.3345 0.3594 1.2544 0.3707 0.3411

eig(A0)

ans = 3×1

3.5595 -0.5115 0.4865

So A0 is not in fact truly PSD. But it is not terribly far off.

Now formulate the problem, combined with all constraints. The unknown will be a vector of

X = optimvar('X',[1,6])

X = 
  1×6 OptimizationVariable array with properties:

  Array-wide properties:
          Name: 'X'
          Type: 'continuous'
    IndexNames: {{}  {}}

  Elementwise properties:
    LowerBound: [-Inf -Inf -Inf -Inf -Inf -Inf]
    UpperBound: [Inf Inf Inf Inf Inf Inf]

  See variables with show.
  See bounds with showbounds.

prob = optimproblem('Description','PSD matrix problem');

% Sadly, squareform does not work here.
L = reshape([X(1),X(2),X(3),0,X(4),X(5),0,0,X(6)],[3,3]);
A = L*L';

prob.ObjectiveSense = 'minimize';
prob.Objective = -trace(A*C) + lambda*norm(A - A0,1); % 1-norm used
prob.Constraints.Aconstraint = -A <= 0;
prob.Constraints.tracelimit = trace(A) <= trmax;

X0.X = rand(1,6);
Xsol = solve(prob,X0)

Solving problem using fmincon.

Feasible point with lower objective function value found.


Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

Xsol = struct with fields:
    X: [1.2594 0.9781 0.6749 4.2153e-05 -9.8384e-06 -2.9467e-05]

How well did it do?

L = tril(ones(3));
L(find(L)) = Xsol.X;

The matrix L is lower triangular.

L

L = 3×3

1.2594 0 0 0.9781 0.0000 0 0.6749 -0.0000 -0.0000

And so Asol was:

Asol = L*L'

Asol = 3×3

1.5862 1.2318 0.8499 1.2318 0.9566 0.6601 0.8499 0.6601 0.4554

trace(Asol)

ans = 2.9982

So the trace of A lies at or near the upper limit for the trace.

Is Asol positive semi-definite? It must be so, because it was formed from a Cholesky product.

format long g
eig(Asol)

ans = 3×1

4.68376370308797e-10 1.74264820884218e-09 2.99823203228449

Sadly, this answer has been posted far too late to help @Vincent. But perhaps others may find it useful. I could have done it as easily without the problem formulation for an optimization, so just as a call to fmincon, but since I'm posting an answer in 2023, I might as well make my answer be one that uses the current capabilities of MATLAB.