---
abstract: 'We provide a simple and efficient algorithm for the projection operator for weighted $\ell_1$-norm regularization subject to a sum constraint, together with an elementary proof. The implementation of the proposed algorithm can be downloaded from the author’s homepage.'
author:
- |
    Weiran Wang\
    Toyota Technological Institute at Chicago\
    `weiranwang@ttic.edu`
date: 'March 1, 2015'
title: 'An $\calO(n\log n)$ projection operator for weighted $\ell_1$-norm regularization with sum constraint'
---

The problem
===========

In this report, we consider the following optimization problem: $$\begin{aligned}
 \label{e:proj}
\min_{\x} & \quad \frac{1}{2} \norm{\x-\y}^2 + \sum_{i=1}^n d_i \abs{x_i}, \\
\text{s.t.} & \quad \x^\top \1=1, \nonumber % \label{e:lincons}\end{aligned}$$ where $\y=[y_1,\dots,y_n]^\top \in\bbR^n$, $d_i\ge 0$, $i=1,\dots,n$, and $\1$ is the $n$-dimensional vector consisting of all $1$’s. This is a quadratic program and the objective function is strictly convex (even though it is non-smooth), so there is a unique solution which we denote by $\x=[x_1,\dots,x_n]^\top$ with a slight abuse of notation.

Notice if $d_1=\dots=d_n$ and the constraint were absent, the problem has a closed form solution known as the soft-shrinkage operator (see, e.g., [@BeckTeboul09a]), which is widely used for solving $\ell_1$-regularized problem in learning sparse representations. But our problem is more involved due to the constraint that couples all dimensions of $\x$. Nonetheless, we give an efficient algorithm with time complexity $\calO(n\log n)$ for this problem using only the KKT theorem.

Our motivation for also comes from sparse coding. @Yu_09a propose the local coordinate coding (LCC) algorithm for learning sparse representations induced by locality. Given a data sample $\uu\in \bbR^n$ and a set of landmark points $\{\vv_j\}_{j=1}^C$ where $\vv_j \in \bbR^n$, $j=1,\dots,C$, the LCC algorithm reconstructs $\uu$ from the landmark points while enforcing the faraway landmark points to contribute less than nearby landmark points (or to have smaller reconstruction coefficients). Let the reconstruction coefficient of $\vv_j$ be $w_j$, $j=1,\dots,C$. Then the optimization problem for these coefficients in LCC is $$\begin{aligned}
\label{e:lcc}
\min_{\w} & \quad \norm{\uu-\sum_{j=1}^C w_j \vv_j}^2 + \lambda \sum_{j=1}^C {\norm{\uu-\vv_j}^2 \abs{w_j} } \\
\text{s.t.} & \quad \sum_{j=1}^C w_j=1, \nonumber\end{aligned}$$ where $\lambda>0$ is some trade-off parameter. The constraint in ensures that the representation is translation invariant. There are different ways of solving this problem, e.g., @ElhamifVidal11a have a similar optimization problem which they solve with Alternating Direction Method of Multipliers [@Boyd_11a]. One simple way of solving is to use the gradient proximal algorithm and its Nesterov’s acceleration scheme (see [@BeckTeboul09a] and the reference therein), where one iteratively takes a short gradient step for the smooth quadratic term and projects the new estimate with the weighted $\ell_1$ regularization term subject to the sum constraint, where the projection operator solves exactly .

The solution
============

We solve the problem using only the KKT theorem [@NocedalWright06a], which states the necessary and sufficient condition[^1] satisfied by the solution $\x$. The Lagrangian of is $$\calL(\x,\alpha)= \frac{1}{2} \norm{\x-\y}^2 +  \sum_{i=1}^n d_i \abs{x_i}+\alpha (\x^\top \1-1),$$ where $\alpha$ is the Lagrange multipliers associated with the constraint. And the KKT system of this problem is

\[e:kkt\] $$\begin{aligned}
x_i-y_i+ d_i+\alpha=0, &\qquad \text{if}\quad x_i>0,\\
x_i-y_i- d_i+\alpha=0, &\qquad \text{if}\quad x_i<0,\\
- d_i \le -y_i+\alpha \le  d_i, &\qquad \text{if}\quad x_i=0, \label{e:abs}\\
\sum_{i=1}^n x_i=1,\end{aligned}$$

where we have used the fact that the sub-differential of $\abs{x}$ is $[-1,1]$ at $x=0$ to obtain .

Denote $y_i^- = y_i- d_i$, $y_i^+ = y_i+ d_i$, $i=1,\dots,n$, which can be computed beforehand. We can then rewrite in terms of $\alpha$:

\[e:alpha\] $$\begin{aligned}
\alpha < y_i^- \; \Longleftrightarrow x_i>0, \label{e:alpha1}\\
\alpha > y_i^+ \; \Longleftrightarrow x_i<0, \label{e:alpha2}\\
y_i^- \le \alpha\le y_i^+ \; \Longleftrightarrow x_i=0, \label{e:alpha3}\\
% \sum_{i:\; \alpha < y_i^-} (y_i^- -\alpha)&\; + \sum_{i:\; \alpha > y_i^+} (y_i^+ -\alpha) =1. \label{e:alpha4}
\sum_{i:\; x_i>0} (y_i^- -\alpha) + \sum_{i:\; x_i<0} (y_i^+ -\alpha) =1. \label{e:alpha4}\end{aligned}$$

Obviously, the Lagrange multiplier $\alpha$ is the key to our problem. Once the value of $\alpha$ is determined, we can easily obtain the optimal solution by setting

\[e:sol\] $$\begin{aligned}
x_i &= y_i^- -\alpha \qquad   \text{if}\;\; y_i^-  > \alpha, \\
x_i &= y_i^+ -\alpha \qquad  \text{if}\;\; y_i^+  < \alpha, \\
x_i &= 0  \qquad\qquad\ \   \text{otherwise}.\end{aligned}$$

We can sort all dimensions of $y_i^-$ and $y_i^+$ together (a total of $2N$ scalars) into an ascending $z$-sequence: $$% z_0= -\infty < z_1 \le z_2 \le \cdots \le z_{2N}\le z_{2N+1}=\infty.
z_1 \le z_2 \le \cdots \le z_{2N}.$$ An important observation is that the $z$-sequence partitions the real axis into $4N+1$ disjoint sets, each being either a single point set $\{z_j\}$, $j=1,\dots,2N$ or an open interval of the form $(-\infty,z_1)$, $(z_j,z_{j+1})$, $j=1,\dots,2N-1$, or $(z_{2N},\infty)$ and the Lagrange multiplier $\alpha$ for the solution must lie in one of them.

We then test each of the $4N+1$ sets as follows. Assuming that $\alpha$ lies in one set, we can use – to conjecture the positive, negative, and zero dimensions of a possible solution $\hat{\x}$. After that, we use to compute a hypothesized value $\hat{\alpha}$ for the Lagrange multiplier, i.e., $$\begin{aligned}
\alpha=\frac{\sum\limits_{i:\; \hat{x}_i>0} y_i^- + \sum\limits_{i:\; \hat{x}_i<0} y_i^+ -1}{\sum\limits_{i:\; \hat{x}_i>0} 1 + \sum\limits_{i:\; \hat{x}_i<0} 1 }.\end{aligned}$$ If the computed $\hat{\alpha}$ indeed lies in the assumed set (a point or an open interval), we have a KKT point and thus the solution.

Since the problem is strictly convex and there exists a unique global optimum, this procedure will find the exact solution with no more than $4N+1$ tests. We can do this efficiently by sorting $y_i^-$ and $y_i^+$ separately ($\calO(n\log n)$ operations) and gradually merging the two sorted sequences (an $\calO(n)$ operation). Therefore the total cost of our procedure for solving is of order $\calO(n\log n)$.

Algorithm \[alg:proj\] gives the detailed pseudocode for solving , whose MATLAB and `C++` implementation can be downloaded at `https://eng.ucmerced.edu/people/wwang5`.

$\y \in\bbR^n$ and $\dd=[d_1,\dots,d_n]$ where $d_i\ge 0$, $i=1,\dots,n$. Sort $\y- \dd$ into $\y^-$: $y_1^- \le y_2^- \le \dots \le y_n^-$.  And sort $\y+ \dd$ into $\y^+$: $y_1^+ \le y_2^+ \le \dots \le y_n^+$. $i\leftarrow 1$,  $j\leftarrow 1$ $s_1 \leftarrow \sum_{i=1}^n y_i^-$,  $s_2 \leftarrow 0$,  $t\leftarrow n$ $\alpha\leftarrow (s1+s2)/t$; **return**

$k\leftarrow i$ $s_1\leftarrow s_1-y_k^-$, $t\leftarrow t-1$, $k\leftarrow k+1$ $\alpha\leftarrow y_i^-$; **return** $left\leftarrow y_i^-$, $i\leftarrow k$ $\alpha\leftarrow y_j^+$; **return** $left\leftarrow  y_j^+$ $s_2\leftarrow s_2+y_j^+$, $t\leftarrow t+1$, $j\leftarrow j+1$ $k\leftarrow i$ $s_1\leftarrow s_1-y_k^-$, $t\leftarrow t-1$, $k\leftarrow k+1$ $\alpha\leftarrow y_i^-$; **return** $left\leftarrow  y_i^-$, $i\leftarrow k$ $s_2\leftarrow s_2+y_j^+$, $t\leftarrow t+1$, $j\leftarrow j+1$ $right\leftarrow y_i^-$ $right\leftarrow y_j^+$ $\alpha\leftarrow (s_1+s_2-1)/t$; **return** $\alpha$ is the Lagrange multiplier of the problem , use to obtain $\x$.

[5]{} \[1\][\#1]{} \[1\][`#1`]{} urlstyle \[1\][doi: \#1]{}

A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. *SIAM J. Imaging Sciences*, 20 (1):0 183–202, 2009.

S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. *Foundations and Trends in Machine Learning*, 30 (1):0 1–122, 2011.

E. Elhamifar and R. Vidal. Sparse manifold clustering and embedding. In J. Shawe-Taylor, R. S. Zemel, P. Bartlett, F. Pereira, and K. Q. Weinberger, editors, *Advances in Neural Information Processing Systems (NIPS)*, volume 24, pages 55–63. MIT Press, Cambridge, MA, 2011.

J. Nocedal and S. J. Wright. *Numerical Optimization*. Springer Series in Operations Research and Financial Engineering. Springer-Verlag, New York, second edition, 2006.

K. Yu, T. Zhang, and Y. Gong. Nonlinear learning using local coordinate coding. In Y. Bengio, D. Schuurmans, J. Lafferty, C. K. I. Williams, and A. Culotta, editors, *Advances in Neural Information Processing Systems (NIPS)*, volume 22. MIT Press, Cambridge, MA, 2009.

[^1]: Strictly speaking, our objective is convex and non-smooth, so the condition is that the zero vector $\0$ lies in the sub-differential at the solution $\x$.
