TFOCS: Projecting on box-affine constraint

In TFOCS, how can we code a projection to a set like this:
$${ x \in \mathbb{R}^n: \alpha^\top x =0, \ x_i \in [0,1]}$$
Would a block structure that separately models the two sets work? If yes, is there a direct way to project to the first (affine equality) set?

Well, we’ve had one crazy person (I’m joking) use CVX to implement their custom project. It’s crazy, of course, because it is an inefficient way to compute a projection! The whole point is to build problems up using combinations of efficient operators and projectors. But it was a very novel projection, so it made a certain kind of sense.

It might be instructive to do it here as well. What you need to be able to do is solve this optimization problem for fixed \bar{x}:

$$\begin{array}{ll}\text{minimize} & |x-\bar{x}|_2 \ \text{subject to} & \alpha^T x = 0 \ & 0\preceq x \preceq \vec{1}\end{array}$$

In CVX, this model looks like:

cvx_begin
    variable x(n)
    minimize(norm(x-xbar))
    alpha' * x == 0
    0 <= x <= 1
cvx_end

Of course, in the end, you would want to come up with a more efficient way to compute this projection, but at least you have something to check it against. It will likely be easier to minimize \frac{1}{2}\|x-\bar{x}\|_2^2 in a custom algorithm.

Thanks :slight_smile: However, I’m looking for TFOCS code. For instance, using two h functions in a block structure, where one is proj_box and the other some sort of weighted proj_simplex. Is this possible? Or it’s better to implement a new home-baked h?

No, even the scalar \alpha^T x is the combination of a linear operator and a standard projection. So this is tfocs_SCD material.

I see. However, in this case dualization seems overkill, since projection to the box-affine set can be directly encoded by a new proj operator. I wonder if TFOCS can somehow figure out when it needs to dualize vs directly compute the projection… Does this make sense?

CVX does intelligently decide whether to solve the primal or dual. But TFOCS includes no such intelligence at all, and likely never will. It is a toolbox for building custom solvers for custom models. I think your best bet is implementing a custom projector.

It’s not obvious, but yes, we can do this, and we can do it rather efficiently (order N log N, where N is the length of the vector). Here’s some barebones code that will do it, and I’ll put this into the next TFOCS release. Let me know if you have feedback. You should test this with the cvx model that Michael suggested to make sure it is working.

The trick is that if we knew the correct Lagrange multiplier for the affine constraint, then the problem decouples and each index is easily solved. To find the right value for the affine constraint, we have to solve a 1 dimensional equation, but this equation is piecewise constant with 2N “turning points”: in between these turning points, we can disregard the box constraints (or “pin” active constraints to the constraint value), and the equation becomes linear and we can solve it exactly.

I don’t know if this has been done before for your exact model, but the trick has been used (and rediscovered) in many settings.

Below is the code to compute the projection (i.e., the proximity operator) of y onto the sets you described in your post.

function x = projBoxAffine( y, a )
n           = length(y);
projBox     = @(x) max( 0, min( 1, x ) );
% Turning points for constraints = 0
T1 = y./a;
% Turning points for constraints = 1
T2 = (y-1)./a;
T = sort(union(T1,T2));
lwrBound = 1;
uprBound = 2*n;
for i = 1:ceil(log2(2*n))
  indx    = round( (lwrBound+uprBound)/2 );
  beta    = T(indx);

  % Our trial solution
  x       = projBox( y - beta*a );
  % Refine beta on the support (ignore constraints): see if it satisifies constraints
  S       = find( x > 0 & x < 1 );
  betaEst = (a(S)'*y(S) + sum(a(x==1)) )/(a(S)'*a(S));
  % Check if bestEst is in the admissible range
  % e.g. if betaEst > beta, is it less than T(indx+1)? and vice-versa
  if betaEst > beta
      if indx == 2*n || betaEst < T(indx+1)
        break;
      else
        lwrBound = indx + 1; % we need to increase beta
      end
  else
      if indx == 1 || betaEst > T(indx-1)
        break;
      else
        uprBound = indx - 1; % we need to decrease beta
      end
  end
end
x       = projBox( y - betaEst*a );
S       = find( x > 0 & x < 1 );
betaEst = (a(S)'*y(S) + sum(a(x == 1)) )/(a(S)'*a(S));
x       = projBox( y - betaEst*a );