How to Formulate a Constraint Which Ensures All Variables Have the Same Sign
The solution from the blog you linked goes as following (Coordinating Variable Signs by Paul Rubin, Web Archive):


  Someone asked me today (or yesterday, depending on whose time zone you
  go by) how to force a group of variables in an optimization model to
  take the same sign (all nonpositive or all nonnegative). Assuming that
  all the variables are bounded, you just need one new binary variable
  and a few constraints.
  
  Assume that the variables in question are $ {x}_{1}, {x}_{2}, \ldots, {x}_{n} $ and that they are all bounded, say $ \forall \; 1 \leq i \leq n, \; {l}_{i} \leq {x}_{i} \leq {u}_{i} $.   If we are going to
  allow variables to be either positive or negative, then clearly we
  need $ {l}_{i} &lt; 0 &lt; {u}_{i} $.   We introduce a new binary variable $ y \in \left\{ 0, 1 \right\} $ and for each $ i $ the constraints
  become:
  
  $$ \forall \; 1 \leq i \leq n, \; {l}_{i} \left( 1 - y \right) \leq {x}_{i} \leq {u}_{i} y $$
  
  If $ y = 0 $ every original variable must be between its lower bound
  and 0, nonpositive. If $ y = 1 $ every original variable must be
  between 0 and its upper bound, nonnegative.
  
  Note that trying to enforce strictly positive or strictly negative
  rather than nonnegative or nonpositive is problematic, since
  optimization models abhor strict inequalities. The only work around I
  know is to change "strictly positive" to "greater than or equal to $
&gt; \epsilon $" for some strictly positive $ \epsilon $, which creates
  holes in the domains of the variables (making values between 0 and $ \epsilon $ infeasible).


The problem with above approach is that it will make a Convex Optimization Problem into Non Convex Optimization Problem because of the multiplication between variables.

For many solvers I guess it will make no difference. Yet if you use Convex Optimization solver or even Disciplined Convex Programming (DCP) based solver (Like CVX) it won't work.

What I suggest doing, in this case, is take advantage that you have only 2 cases. One where all variables are Non Negative and the other when all of them are Non Positive.

Then all you need is to solve the problem with 2 different sets of constraints (Forming 2 different problems). When the constraints are:


$ \forall \; 1 \leq i \leq n, \; {l}_{i} \leq {x}_{i} \leq 0 $
$ \forall \; 1 \leq i \leq n, \; 0 \leq {x}_{i} \leq {u}_{i} $


Then chose the solution with the lowest objective value.  

The main advantage is while the solution using Binary Variable (Which probably will do the same) require special MILP / MIQP solver the approach I suggested, for Linear Least Squares problem only require Non Negative Least Squares solver which is easier to find (Usually faster).

Example with Code

Let's say we want to solve the following problem:

$$\begin{aligned}
\arg \min_{x} \quad &amp; \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \\
\text{subject to} \quad &amp; \operatorname{sign} \left( {x}_{i} \right) = \operatorname{sign} \left( {x}_{j} \right) &amp;&amp; \forall i, j \in \mathcal{S}
\end{aligned}$$

Where $ \mathcal{S} = \left\{ i \mid {x}_{i} \neq 0 \right\} $, namely the set of all indices for which the element of $ x $ isn't zero.
So basically we say all elements should have the same sign or be zero.

Let's say we have a building block for solving Non Negative Least Squares. For instance one could use MATLAB's lsqnonneg().
To solve the problem above one could solve it one time for $ A $ and the other for $ - A $. Then take the argument which minimize the objective (If we take the argument of $ - A $ we need to negate it).

So the MATLAB function would look like:



function [ vX ] = SolveLsSameSign( mA, vB )

[vX1, resNorm1] = lsqnonneg(mA, vB);
[vX2, resNorm2] = lsqnonneg(-mA, vB);

if(resNorm1 &lt; resNorm2)
    vX = vX1;
else
    vX = -vX2;
end

end


I created a simple script to compare results to CVX using the Mixed Integer Quadratic Programming (MIQP) which is given by:



cvx_solver('Mosek');
% cvx_solver('Gurobi');

cvx_begin()
    % cvx_precision('best');
    variable vX(numCols, 1);
    variable varY(1) binary;
    minimize( norm(mA * vX - vB, 2) );
    subject to
        vX &gt;= (1 - varY) * lowerBound;
        vX &lt;= varY * upperBound;
cvx_end


Pay attention that it requires MOSEK or Gurobi to be installed in order to solve this MIQP problem.

The full MATLAB code is available on my StackExchange Signal Processing Q52099 GitHub Repository.