# Algorithm to solve $L_1$ optimization of $\sum_i ||\mathbf{A_i x} - \mathbf{b_i}||_1$

Is there is an efficient algorithm to solve the following optimization:

$$\mathbf{x}^* = \arg\min_\mathbf{x}\sum_i ||\mathbf{A_i x} - \mathbf{b_i}||_1$$

for given $$\mathbf{b_i}, \mathbf{A_i}\ \forall i$$, where $$|| \mathbf{y} ||_1$$ is the $$L_1$$ norm (i.e. $$\sum_j |y_j|$$) and $$\mathbf{A_i}$$ are big matrices which we don't have their explicit representations, only access to operations $$\mathbf{A_i y}$$ and $$\mathbf{A_i^T z}$$?

I can solve the problem for small matrices $$\mathbf{A_i}$$ using linear programming using the same trick from here, but if the matrices are large, there is no obvious way (at least for me) to solve it using linear programming.

Any ideas and suggestions are welcome.

• In other words, you are looking for an algorithm which is more efficient than linear programming. Dec 29 '18 at 11:35
• If you stack your matrices $\mathbf{A_i}$ and your vectors $\mathbf{b_i}$ together, can’t you reduce your problem to the case $n = 1$ (i.e., there is only one matrix/vector pair)? Dec 29 '18 at 11:36
• I'm looking for algorithms that are more memory-efficient than linear programming. My problem is that the matrices are very large and it's almost impossible to store in memory. I only have access to the operations $A\cdot$ and $A^T\cdot$ Dec 29 '18 at 17:07

I have found the answer for this problem. One way is, as suggested in the comment, to stack the matrix $$\mathbf{A_i}$$ and $$\mathbf{b_i}$$ together and apply first-order convex optimization algorithm, such as Primal-Dual Hybrid Gradient (PDHG), inexact ADMM, or their variants.

For PDHG, the problem to solve is

\begin{align} \mathrm{min}_{\mathbf{x},\mathbf{y}} &~ F(\mathbf{y}) + G(\mathbf{x}) \\ \mathrm{s.t.} & ~ \mathbf{Ax} = \mathbf{y} \end{align}

where

\begin{align} F(\mathbf{y}) &= ||\mathbf{y} - \mathbf{b}||_1 \\ G(\mathbf{x}) &= 0 \\ \mathbf{A} &= [\mathbf{A_1}^T, \mathbf{A_2}^T, ..., \mathbf{A_n}^T]^T \\ \mathbf{b} &= [\mathbf{b_1}^T, \mathbf{b_2}^T, ..., \mathbf{b_n}^T]^T. \\ \end{align}

Thus, the convex conjugate of $$F(\mathbf{y})$$ is

\begin{align} F^*(\mathbf{y}^*) &= \mathrm{sup}_{\mathbf{y}} \left[\left<\mathbf{y},\mathbf{y}^*\right> - F(\mathbf{y})\right] \\ F^*(\mathbf{y}^*) &= \left<\mathbf{b},\mathbf{y^*}\right> + \mathcal{I}_{\mathbf{y^*}<\mathbf{b}} \end{align} where $$\mathcal{I}_\mathcal{C}$$ is an indicator function with value $$0$$ if the condition $$\mathcal{C}$$ fulfilled, otherwise $$\infty$$.

As the algorithm uses proximal operator, we need to get the proximal operator for $$F$$ and $$G$$,

\begin{align} \mathrm{prox}_{F^*\sigma}(\mathbf{y}^*) &= \arg\min_{\mathbf{y}'}\left\{\frac{||\mathbf{y}'-\mathbf{y}^*||^2}{2\sigma} + F^*(\mathbf{y}^*)\right\} \\ \mathrm{prox}_{F^*\sigma}(\mathbf{y}^*) &= \max(\mathbf{y}^* + \sigma \mathbf{b}, \mathbf{b}) \\ \mathrm{prox}_{G\tau}(\mathbf{x}) &= \mathbf{x} \\ \end{align}

Thus, the algorithm in this case reads

\begin{align} \mathbf{y}^{n+1} &= \max(\mathbf{y}^n + \sigma \mathbf{A} \mathbf{w}^n + \sigma \mathbf{b}, \mathbf{b}) \\ \mathbf{x}^{n+1} &= \mathbf{x}^n - \tau \mathbf{A}^T \mathbf{y}^{n+1} \\ \mathbf{w}^{n+1} &= \mathbf{x}^{n+1} + \theta (\mathbf{x}^{n+1} - \mathbf{x}^n), \end{align}

with $$\theta = 1$$, and $$\sigma \tau ||A||^2 < 1$$.

All the operations do not need the direct access to the element of the matrix $$\mathbf{A}$$ and only require the operation of $$\mathbf{A} \cdot$$ and $$\mathbf{A}^T \cdot$$.

One possible approach would be to use stochastic gradient descent. This will be memory-efficient, as you desired. In particular, in each iteration we randomly pick an index $$i$$, then update $$x$$ based on the gradient of the function $$f_i(x) = \| A_i x - b_i \|_1$$, and repeat many times. Each iteration only requires access to a single matrix $$A_i$$ at a time.