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.
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}} &amp;~ F(\mathbf{y}) + G(\mathbf{x}) \\
\mathrm{s.t.} &amp; ~ \mathbf{Ax} = \mathbf{y}
\end{align}
where
\begin{align}
F(\mathbf{y}) &amp;= ||\mathbf{y} - \mathbf{b}||_1 \\
G(\mathbf{x}) &amp;= 0 \\
\mathbf{A} &amp;= [\mathbf{A_1}^T, \mathbf{A_2}^T, ..., \mathbf{A_n}^T]^T \\
\mathbf{b} &amp;= [\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}^*) &amp;= \mathrm{sup}_{\mathbf{y}} \left[\left&lt;\mathbf{y},\mathbf{y}^*\right&gt; - F(\mathbf{y})\right] \\
F^*(\mathbf{y}^*) &amp;= \left&lt;\mathbf{b},\mathbf{y^*}\right&gt; + \mathcal{I}_{\mathbf{y^*}&lt;\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}^*) &amp;= \arg\min_{\mathbf{y}'}\left\{\frac{||\mathbf{y}'-\mathbf{y}^*||^2}{2\sigma} + F^*(\mathbf{y}^*)\right\} \\
\mathrm{prox}_{F^*\sigma}(\mathbf{y}^*) &amp;= \max(\mathbf{y}^* + \sigma \mathbf{b}, \mathbf{b}) \\
\mathrm{prox}_{G\tau}(\mathbf{x}) &amp;= \mathbf{x} \\
\end{align}
Thus, the algorithm in this case reads
\begin{align}
\mathbf{y}^{n+1} &amp;= \max(\mathbf{y}^n + \sigma \mathbf{A} \mathbf{w}^n + \sigma \mathbf{b}, \mathbf{b}) \\
\mathbf{x}^{n+1} &amp;= \mathbf{x}^n - \tau \mathbf{A}^T \mathbf{y}^{n+1} \\
\mathbf{w}^{n+1} &amp;= \mathbf{x}^{n+1} + \theta (\mathbf{x}^{n+1} - \mathbf{x}^n),
\end{align}
with $\theta = 1$, and $\sigma \tau ||A||^2 &lt; 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$.